Study of final-state interactions of protons in neutrino-nucleus scattering with INCL and NuWro cascade models

The modeling of neutrino-nucleus interactions constitutes a challenging source of systematic uncertainty for the extraction of precise values of neutrino oscillation parameters in long-baseline accelerator neutrino experiments. To improve such modeling and minimize the corresponding uncertainties, a new generation of detectors is being developed, which aim to measure the complete final state of particles resulting from neutrino interactions. In order to fully benefit from the improved detector capabilities, precise simulations of the nuclear effects on the final-state nucleons are needed. This article presents the study of the in-medium propagation of knocked-out protons, i.e., final-state interactions (FSI), comparing the NuWro and INCL cascade models. The INCL model is used here for the first time to predict exclusive final states of neutrino interactions. This study of INCL in the framework of neutrino interactions features various novelties, including the production of nuclear clusters (e.g., deuterons, $\alpha$ particles) in the final state. The paper includes a complete characterization of the final state after FSI, comparisons to available measurements of single transverse variables, and an assessment of the observability of nuclear clusters.

The modeling of neutrino-nucleus interactions constitutes a challenging source of systematic uncertainty for the extraction of precise values of neutrino oscillation parameters in long-baseline accelerator neutrino experiments. To improve such modeling and minimize the corresponding uncertainties, a new generation of detectors is being developed, which aim to measure the complete final state of particles resulting from neutrino interactions. In order to fully benefit from the improved detector capabilities, precise simulations of the nuclear effects on the final-state nucleons are needed. This article presents the study of the in-medium propagation of knocked-out protons, i.e., final-state interactions (FSI), comparing the NuWro and INCL cascade models. The INCL model is used here for the first time to predict exclusive final states in measured neutrino interaction cross sections. This study of INCL in the framework of neutrino interactions features various novelties, including the production of nuclear clusters (e.g., deuterons, α particles) in the final state. The paper includes a complete characterization of the final state after FSI, comparisons to available measurements of single transverse variables, and an assessment of the observability of nuclear clusters.

I. INTRODUCTION
The neutrino oscillation physics enters the precision era: notably, the NOVA [1] and T2K [2] long-baseline experiments are featuring measurements of the neutrino mixing angle θ 23 and of the largest mass splitting in the atmospheric sector at few percent precision. Sensitivity studies combining future T2K and NOVA data together with measurements at reactor experiments, show the potential of more than 3σ significance for possible chargeparity (CP) violation hints and mass ordering determination [3,4]. The next-generation experiments DUNE [5] and HyperKamiokande [6] are aiming to establish the mass ordering and possibly discover charge-parity violation with 5σ significance, as well as to measure the value of the CP-parameterizing phase (δ CP ) with a precision better than 15 degrees. Such results will be enabled by unprecedented statistics of produced and detected neutrinos, requiring an exceptionally robust and precise control of systematic uncertainties.
The largest and most complex systematic uncertainty in present neutrino long-baseline experiments stems from the modeling of neutrino-nucleus interactions. The so called "near detectors", placed near the neutrino source before any standard neutrino oscillation can occur, are designed to characterize the neutrino flux and to measure the neutrino-nucleus cross section in order to tune the interaction models and minimize the corresponding uncertainties. In order to cope with the increasing needs in precision, a new generation of near detectors is being developed based on the concept of a precise and exclusive reconstruction of all the final-state particles produced in neutrino-nucleus interactions. This concept is at the core of the design of the upgrade of T2K near detector ND280 [7]. A new active target [8] will allow 3-dimensional reconstruction capabilities for low momentum tracks. Such a detector also enable the reconstruction of neutron kinetic energy by measuring time-of-flight between the neutrino interaction vertex and the first neutron secondary interaction [9]. The sensitivity of this type of detector for the most relevant systematic uncertainties in neutrino-nucleus interactions below 1 GeV has been recently documented in Ref. [10]. Notably, the exclusive reconstruction of the hadronic part of the final state of neutrino-nucleus interactions permits a more precise reconstruction of the neutrino energy on an event-byevent basis, yet it poses new challenges on the modeling of such interactions. Neutrino-nucleus interaction models have historically been developed and tuned to describe inclusive processes, thus they characterize the cross section as a function of the outgoing lepton kinematics. A new effort is ongoing to expand their predictivity to the kinematics of the hadrons for an exclusive description of the final state [11,12].
Most of the available neutrino-interaction models and Monte Carlo simulations describe the initial nuclear state and the fundamental interaction separately from the re-interactions of the final-state hadrons with the nucleus. Typically, final-state interactions (FSI) of resulting hadrons are simulated with a cascade mechanism. It is worth mentioning that many nuclear models, originally developed to describe electron-nucleus scattering, include the effect of an outgoing hadrons distortion in a given nuclear potential while solving the lepton-nucleus interaction [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30]. Extensive feasibility studies of implementing such models into generators are nowadays performed in the community, e.g., in Ref. [30]. While such an improvement is certainly needed to cope with the precision of the next generation of long-baseline experiments, we focus on FSI modeling with currently available tools given the needs of running experiments, notably the analysis of the data from the upgraded ND280. This paper presents the study of FSI with two different nuclear models implemented in the IntraNuclear Cascade Liege (INCL) code [31][32][33][34] and in NuWro [35]. Fundamental neutrino-nucleus interactions without pion production are studied, where pions in the final state can only be produced through nucleon FSI. This work can be expanded in the future to antineutrino interactions (focusing on neutrons FSI) and to (anti)neutrino interactions with pion production (focusing on pion FSI). The study presented here and its possible further developments aim at evaluating in a detailed way the possible uncertainties inherent to the FSI simulation in neutrino-nucleus scattering. This study is also a step forward into a more precise and complete implementation of FSI effects in Monte Carlo simulations of neutrino-nucleus interactions.
An important novelty in this paper is a discussion of the production of nuclear clusters (like α particles and deuterons) which is modeled in the cascade mechanism of INCL, for the first time, in neutrino-nucleus interactions. Such predictions may significantly impact the analysis and interpretation of the experimental data, notably for the identification of low momentum particles and the measurement of energy deposits around the neutrino vertex (also known as vertex activity).
The NuWro and INCL nuclear models are described in Section II. The procedure to simulate fundamental neutrino-nucleus interaction events and, subsequently, simulate the intra-nuclear cascade with different models in INCL or NuWro, is described in Section III. Section IV introduces the studied variables and the analysis strategy. Section V presents the results from the simulations of the different nuclear models and their comparison. We study different variables to characterize and quantify the impact of FSI on the kinematics of the leading proton, with particular focus on Single Transverse Variables (STV) introduced in Ref. [36] and on the visible energy deposited around the vertex. Section VI re-ports a reasoned comparison of the studied simulations with available measurements of STV. Finally the main conclusions are reported in Section VII.

A. NuWro
Since 2005, the theoretical group of the University of Wroc law has developed NuWro as a comprehensive Monte Carlo lepton-nucleus event generator [37], optimizing it for use in accelerator-based neutrino oscillation experiments, i.e., the few-GeV energy region. Depending on the energy transferred from the interacting leptonic probe to the hadronic system, NuWro provides quasielastic [38], hyperon production [39], single-pion production and more inelastic channels (DIS) [40] for scattering off free nucleons. After including complex nuclear targets, additional channels such as two-body processes [41], coherent pion production [42], and neutrino scattering off atomic electrons [43] are included. The framework utilizes various nuclear models to provide predictions for the dynamics of target nucleons (e.g., global or local Fermi gas, spectral functions [13,15], or a momentumdependent nuclear potential [38]). Finally, FSI are simulated by an intranuclear cascade which propagates the outgoing nucleons [44] and produced pions [35] through the residual nucleus. In the context of this work, the technical aspects of modeling quasielastic neutrino-nucleus scattering and the subsequent FSI are emphasized. More details on the used NuWro version (19.02.2) could be found in Refs. [44,45].
In the quasielastic interaction channel, nuclear modeling enters utilizing the plane-wave impulse approximation that factorizes the one nucleon knock-out processes into the interaction on a single off-shell nucleon convoluted with a particular hole spectral function, i.e., the probability of leaving the residual system with specific excitation energy and recoil. The approach relies on the calculation by O. Benhar et al. [13] that takes into account the electron scattering input to the single-particle wave functions and adding the correlated part evaluated within the local-density approximation. Additionally, in this model, the prescription by A. Ankowski et al. [46] is applied to go beyond the factorized picture and account for the effects of distorting the final nucleon wave function by an optical potential. Alternatively, the target nucleons can be treated as constituing the ideal Fermi gas, parametrized through nuclear density or its average value and referred to as local (LFG) or relativistic Fermi gas (RFG), respectively. Finally, the primary interaction vertex is constrained by the conserved vector current (CVC) and partially conserved axial current (PCAC) hypotheses. The vector form factors are provided by the BBBA05 parametrization [47], while the axial form factor has a dipole shape with g A = 1.267 and the axial mass parameter M A = 1.03 GeV/c 2 , according to the discussion in Ref. [48].
Modeling final-state interactions is a challenging manybody problem that bears a tension between numerical efficiency and accuracy of nuclear calculations. NuWro solution is based on seminal papers by N. Metropolis et al. [49,50], which describe an algorithm of the space-like cascade model, and applied up-to-date physics ingredients. In this approach, mean free paths are attributed to the particles propagated in straight lines with steps of ∆x through a continuous medium. Such Monte Carlo sampling uses the standard non-interaction probability formula where λ = (ρσ) −1 is the mean free path calculated locally, expressed in nuclear density ρ and an effective interaction cross section σ. The maximal step of ∆x = 0.2 fm is sufficient to grasp the structure of commonly used density profiles. By default, the nucleons constituting the nuclear medium originate from the LFG model, and therefore, meet its Pauli blocking rules (applied on an event-by-event basis). The cascade terminates when all the moving hadrons leave the nucleus or do not have enough kinetic energy and are stuck in nuclear potential (with the separation energy of 7 MeV). The remnant nucleus is in an excited state, and its deexcitation is not modeled. In Ref. [44], the nucleon part of the NuWro cascade has been exhaustively tested, aiming to reproduce the nuclear transparency in exclusive (e, e p) scattering experiments. The essence of this model lies in nucleon-nucleon cross sections, which replicate the PDG dataset [51], the fraction of single-pion production adjusted to follow the fits of Ref. [52], and the center-of-momentum frame angular distributions of Ref. [53]. Additionally, the cross sections are modified with in-medium corrections [54,55] and twonucleon correlation effects [44]. Finally, the pion-nucleon interaction dynamics is taken from the model of L.L. Salcedo et al. [56]. This aspect, together with the formation zone effect for the inelastic scattering channels, has been presented and compared to data in Ref. [35].
The standard version for the projectile employs the impact parameter formalism taking into account the Coulomb distortion. A 'working sphere', where all events take place, is defined with radius R max : where R 0 and a are radius and diffuseness of the target nucleus density respectively. For carbon, R max = 5.7 fm. If the trajectory of the projectile, defined by a randomly chosen asymptotic impact parameter and modified by the Coulomb field of the target, intersects the working sphere, the projectile enters the nucleus and can interact with the nucleons or not (a "transparent" event). The simulation of neutrino interactions in INCL is not yet available. In this study, the neutrino-nucleon interaction inside the nucleus is enforced by using the neutrino vertex from the NuWro code. Then, the produced particles from the neutrino-nucleon interaction may undergo final-state interactions and initiate a cascade. The exact matching procedure to interface NuWro and INCL will be described in Sec. III. The overall normalization is taken from the total cross section calculated by NuWro. The INCL nuclear model is essentially classical, with some additional ingredients to mimic quantum effects. Each nucleon in the nucleus has its position and momentum and moves freely in a square potential well. The radius of the potential well depends on the kinetic energy of the nucleon. Nucleon momenta are distributed uniformly in a sphere whose radius is defined by the maximal Fermi momentum. Position and momentum are correlated in INCL. In such a classical picture, the particle is allowed to move in a sphere with the maximum radius fixed by its momentum, and so its position is sampled in this sphere. This picture has been refined to take into account the quantum properties of the wave functions. Based on a HFB (Hartree-Fock-Bogoliubov) formalism, the correlation is still valid but less strict, i.e. the nucleon has a non-zero probability for going beyond the maximum radius. Further details can be found in Ref. [33].
The cascade reaction can be described as an avalanche of independent binary collisions. Different types of events inside the cascade could occur: collisions, decays and, at the surface, reflections or transmissions. The particles travel along straight lines and the different possible fates are calculated: • two particles reach the minimal distance to interact, • a particle hits the border of the potential well then it deflects back or leaves the nucleus, • a particle decays (e.g. a ∆ resonance or ω meson).
The fate with the shortest time is chosen.
Nucleons that do not participate in the cascade are defined as spectators. In the beginning, all nucleons are spectators except for the projectile. A spectator could become a participant while interacting with the projectile or another participant. To prevent spontaneous nucleon emission (nuclear boiling), the spectator nucleons of the target cannot interact between themselves. The projectile feels the nuclear potential of the target. A participant nucleon can become a spectator again if its energy decreases below a threshold (details are in Ref. [31]). While such an option has been included to improve the agreement with some sets of data [31], the modeling of very low energy nucleons is still an open issue for cascade mechanisms. In the scope of this paper, this feature is disabled. Participants can be absorbed at the end of the cascade when the stopping time is reached and the nucleus is thermalized through equipartition of energy. In this case, the participant's energy is left in the nucleus, putting it in an excited state.
The cascade lasts until one of these conditions is satisfied: • the stopping time determined by the model is reached, • there are no more participants, • the mass number of the target is less than 4, • the projectile leaves without an interaction (transparent event).
Concerning a participant nucleon at the surface, it can be reflected or emitted. It leaves the nucleus if its energy is higher than the Fermi energy plus the value of its separation energy that is taken from mass tables based on experimental data [85]. However, a notable feature of INCL is that an outgoing nucleon (subject to FSI) with some probability could clusterize with other nucleons and leave the nucleus as a nuclear cluster (e.g., α particle).
The emission of nuclear clusters was extensively compared with representative experimental data. One can find more information about cluster emission in INCL in Ref. [32]. In addition to quantum effects embedded in the elementary cross sections used and the one previously quoted with the refined position-momentum correlation related to the quantum behavior of the wave function, Pauli-blocking is considered and checked. If the collision is blocked, the particle propagates until it reaches the next allowed interaction in the cascade and the products of that reaction are further propagated. There are two main Pauli blocking models implemented in INCL: the strict one, forbidding the interaction if the projectile momentum is below the Fermi momentum, and the statistical model that considers only nearby nucleons in the phase-space volume and acts upon the calculated occupation probability. The default option applied in INCL is a compromise between the strict Pauli blocking for the first collision and the statistical one for the subsequent ones. Motivation and details of this procedure can be found in Ref. [31,86]. In this paper, the first interaction is the neutrino interaction with the nucleus taken from NuWro. So for the first interaction, the Pauli blocking from NuWro is used, then for subsequent proton interactions statistical Pauli blocking is applied.

III. SIMULATION PROCEDURE
The study focuses on the comparison between NuWro and INCL of proton-induced FSI cascades in Charged-Current Quasi-Elastic (CCQE) neutrino interactions on Carbon using the T2K neutrino energy flux at the near detector from Ref. [87]. The initial state and the neutrino-nucleus interaction is simulated with NuWro with the Spectral Function (SF) approach. About 350000 CCQE events have been simulated with NuWro. In each event, the leading proton exiting the interaction of a neutrino on a neutron is then injected in the INCL nuclear model and FSI is simulated with INCL. The results are compared to the FSI cascade simulation done in NuWro.
In the INCL simulation, the neutron with the closest momentum-vector to the initial neutron simulated by NuWro for the neutrino interaction is selected. This procedure is independent of the exact values of neutron momenta in NuWro and INCL and is applied event-byevent for the whole sample. The chosen INCL neutron is then substituted with the outgoing particles (muon and proton) preserving their kinematic properties calculated from the neutrino-neutron interaction by NuWro. The Pauli blocking in the neutrino interaction is simulated by NuWro, while the Pauli Blocking along the proton FSI cascade is simulated differently by NuWro and INCL, as explained in Sec. II.
Although driven by the same experimental charge density distribution, NuWro and INCL differ in their approaches to model the nucleus that can result in small inconsistencies within our algorithm. The momentum distribution of nucleons in SF and INCL is quite different, as can be seen in Fig. 2: the high-momentum tail due to nucleon-nucleon correlations is missing in INCL (work is ongoing to include them, and we defer such results for a future paper). Another discrepancy stems from the lack of correlation between position and momentum in the spectral function approach, while having the Hartree-Foch-Bogoliubov formalism implemented in INCL, as presented in Fig. 3. Thus, the INCL position distribution sampled starting from the simulated NuWro interactions is slightly different than the intrinsic nucleon distribution in INCL, as shown in Fig. 4. In general, neutrons with high momentum in the SF tail tend to be attributed to the peripheral region of the nucleus, according to the INCL nuclear model shown in Fig. 3, while in the NuWro SF model, there is no intrinsic positionmomentum correlation for the target neutrons. Still, as shown in Fig. 4, the chosen INCL neutron tends to be at smaller radius with respect to the general neutron distribution in INCL: this is due to the fact that the momentum-vector (and not the momentum magnitude) is used to find the INCL neutron which best matches the NuWro neutron. It is worth mentioning that NuWro uses the local Fermi gas for nuclear model in the cascade itself. This model does include a momentum-position correlation, which is presented at the bottom of Fig. 3. One can see that there exist a contrary dependence relative to the INCL solution, which is an important point in comparing these two cascades regardless of the primary vertex simulation.
In order to check that the results on FSI characterization are robust against the approximations of the simulation procedure above and, more in general, against assumptions on the initial state nuclear model, similar studies are also performed using Relativistic Fermi Gas and using a simulation representing the INCL nuclear model in the initial state. This is discussed in Ap- pendix A.
In the factorized models considered here, the FSI interaction cannot change the fundamental neutrino-nucleus interaction cross section. Thus for all studies presented here the NuWro cross section from SF (or from RFG for the results in Appendix A) is considered.
The procedure to match NuWro and INCL models described above has been made possible by the high-level of modularity of Monte Carlo implementations. It also opens the road to possible further improvement of FSI modeling into such simulation programs.

IV. ANALYSIS STRATEGY
In this section, the analysis strategy to compare and characterize the FSI effects in INCL and NuWro models is described.
Various channels are possible, depending on the particles leaving the nucleus: their probability is quantified in the two models. The FSI effects on the leading proton are then characterized by comparing its kinematics before and after FSI. Additionally, the Single Transverse Variables (STV), as introduced in Ref. [36], are stud- ied and compared between the two models. Resolution effects or thresholds are not applied in the results of Sec. V A and some of the studied variables are not observable experimentally. Comparison to data is deferred to Section VI, where direct comparison to STV measurements from T2K [88] and MINERvA [89] are shown. In this case, acceptance cuts are applied to match the phase space covered by the two experiments.
The study focuses on the following STV: where p p T is the component of the proton momentum projected into the plane transverse to the neutrino direction (transverse component) and p µ T is the transverse component of the muon momentum. The illustration of the STV definition is presented in Fig. 5. The variable δ p T could be considered as the "missing transverse momentum" and, in absence of FSI in quasi-elastic events, it represents the Fermi motion of the initial nucleon.
The variable δα T is especially sensitive to the FSI of the leading proton. In transparent events (where the proton leaves the nucleus without FSI), the δα T distribution depends only on the stochastic Fermi motion of the initial neutron and therefore is expected to be uniform. FSI tends to decrease the proton transverse momentum, inducing larger values of missing transverse momentum and δα T .
The production of clusters is an important novelty brought by INCL: no other cascade model used within neutrino studies is able to simulate cluster production. Not only the kinematics of the simulated clusters but also their identification probability in scintillating detectors is studied. To simulate clusters' interactions inside a detector, a simple Geant4 model [64,90,91] 1 is built. A uniform, fully active block of hydrocarbon (with 1.06 g/cm 3 density corresponding to polystyrene) with dimensions big enough to contain a whole track (1900×1900×3000 cm 3 size) is simulated. The clusters are generated with momentum distribution as simulated by INCL FSI production.
A simplified analysis is performed to gauge the observability of clusters: an identification algorithm is developed based on the path traversed by clusters and the energy deposited in the detector, as shown in the simple schematic of Fig. 6.
The main idea of this algorithm is to try to identify the type of the particle based on its ionization curve. The measured deposited energy is used as an estimation of the particle momentum along its trajectory. The measured local deposits of energy (dE/dx) along the track are compared to the expected ones for different clusters {α, D, T, 3 He, p} at the estimated momenta. Finally the algorithm selects the cluster hypothesis that best matches the observed dE/dx along the track.

×
... Actual identification capabilities in an experiment depend highly on the detector geometry, granularity, and scintillator material. Here we consider a detector granularity of 1 cm 3 cubes, corresponding, for instance, to the geometry of the ND280 upgrade scintillating target su-perFGD [7]. We do not take into account any reflecting material or border effects.
Tracks are simulated as straight lines parallel to one of the cube faces and starting from the cube boundary. The last part of the track, which is shorter than 1 cm, is not used in the analysis since the exact length of the last step is unknown, given the simulated granularity. Thus we effectively remove the Bragg peak from this simplified identification algorithm. In more sophisticated analyses, the Bragg peak study could bring further information to help particle identification. Moreover, at the end of the track, a fraction of events undergoes inelastic interactions with the creation of secondary particles. We leave the investigation of such secondary interactions and their detector signature for future work: we do not include the energy of secondary particles in the analysis since we remove the last cube.
We consider only the energy deposited by the primary particle by ionization and we apply Birks quenching. Depending on the material, an additional overall suppression factor should be considered in energy to scintillation light conversion but is omitted in the present study.
Visible energy in a given simulation step is therefore calculated as where E dep step is the energy deposited in the step, k B is the Birks coefficient, and L step is the step length. Birks coefficient for protons is taken from Ref. [92], where k B = 0.0208 (cm/MeV); according to Ref. [93], the coefficients for protons and deuterons are assumed to be the same. The same Birks coefficient is also assumed for tritium. For α particles the results from Ref. [94] are used, where k B = 0.0085 (cm/MeV).
The total visible energy is evaluated as the sum of the visible energy in each cube is This visible energy is a proxy for the total kinetic energy of the particle. However, it is not a perfect estimator since the energy lost in the last cube, where inelastic events may happen, is discarded from the analysis. As an analogous estimator of the remaining kinetic energy of the cluster at each step along the track, the total 'residual' visible energy is estimated as The distribution of the visible energy deposited in each step E vis i , as a function of the left visible energy E vis,res i , is used to build the expectations for each of the five particle hypotheses {α, D, T, 3 He, p}. For each particle k, the mean of the simulated visible energy at a given step is estimated (E exp,k i ) for the corresponding value of remaining visible energy (E vis,res i ). The width of such distribution is also evaluated (σ k i ). Finally for each step of the simulated cluster track, the visible energy is compared to the expected one for each particle hypothesis: a χ 2 is built for each particle hypothesis and the hypothesis with lowest χ 2 is chosen to identify the particle.

V. RESULTS
In Tab. I the fraction of different final states produced by INCL and NuWro FSI simulation from CCQE interactions are reported. In the NuWro SF simulation, due to short-range correlations, the neutrino vertex contains two outgoing protons in 15% of events. In the INCL vertex, only the leading proton, defined as the proton with higher kinetic energy, is retained and triggers the cascade. We have also tested that complete removal of events with 2-protons at the neutrino vertex does not impact the conclusions on the characterization of the FSI cascade of the leading proton in NuWro and INCL. We consider only the leading proton and its secondary hadrons for the channel characterization in Tab. I. If the proton has the same energy before and after FSI, the event is characterized as "no FSI" (a transparent event).
INCL features an evident enhancement of events with absorption of the proton and only a muon in the final state: 7.7% of the total events in INCL against less than 0.1% in NuWro. Indeed, in the NuWro cascade, the vast majority of particles produced, and not Pauli blocked, leave the nucleus. The INCL nuclear model features larger probability of reabsorbing particles in the nucleus during FSI. This tendency is also confirmed by the smaller fraction of events with more than 1 proton in the final state (multinucleon production by FSI results into 25.6% of total events in NuWro and 9% in INCL).
We have also compared both INCL and NuWro cascades against existing transparency data. Here we define transparency as the percentage of events without FSI interaction on the leading proton exiting from the interaction vertex. Fig. 7 shows that both models are in good agreement with data while having a divergence in their prediction. Especially for the low momentum protons, INCL features smaller transparency. The difference between the two models is dominated by the impact of the short-range correlations on the FSI definition, as also shown in Refs. [44] and [65]. In the region of interest for The larger FSI strength in INCL suggests a larger dissipation of energy across the nucleus through interactions that tend to be more 'soft'. On the other hand, the nuclear model of INCL includes the probability to form clusters during the attempt of the nucleon to leave the nucleus, as explained in Sec. III. Indeed, the events with no proton in the final state are in the large majority due to charge exchange in NuWro (91% of neutron production) while in INCL the probability of cluster production and neutron production in events without protons in the final state is similar (around 30% each). The leading proton momentum before and after FSI is shown in Fig. 8. Events with the production of other particles, e.g., clusters, but no protons in the final state, are included only in the distribution of proton momentum before FSI. The proton momentum before FSI has by construction the same distribution for NuWro and INCL, as explained in Sec. III. In Fig. 8 the relative fraction of the different final states as a function of proton momentum before and after FSI is also shown. A larger diversity of channels is visible in the INCL case. As previously discussed, the smaller transparency of INCL induces a smaller number of events with at least one proton in the final state after FSI. In particular, all the events with full proton absorption are concentrated at low proton momentum before FSI, while the events with cluster production populate the distribution just below the main peak of initial proton momentum (around 350 MeV/c). Again, due to the smaller nuclear transparency of INCL, the events where the proton leaves the nucleus are much less often accompanied by additional nucleons. Thus in NuWro, the protons at very low momenta after FSI are almost always accompanied by other nucleons, while in INCL, they are mostly events with only one proton.
More direct FSI quantification on the leading proton is done by comparing the kinematics of the proton before and after FSI, as shown in Fig. 9. As expected, FSI decelerate the leading proton and induce smearing compared to the pre-FSI angular distribution that is peaked in the forward direction (but this effect is less evident in INCL). Interestingly, INCL tends, event-by-event, to decelerate protons more but still the distribution of momenta of all the protons leaving the nucleus after FSI tends to be larger in INCL, since the low momenta protons are absorbed in the nucleus.
In order to disentangle FSI effects from the physics of the initial nuclear state, STV are used. As explained in Sec. IV, the identified experimental observable most sensitive to FSI is δα T , while δp T is primarily sensitive to the initial nuclear state. This is confirmed by the results in Fig. 10, where such distributions are compared between NuWro and INCL. The shape of δp T is very similar between the two models since the initial nuclear momentum is taken from the same model. The most visible feature is the suppression of the large δα T values in INCL, corresponding to a suppression of low momentum protons. This conclusion is robust, independently of the nuclear model used for the fundamental interaction: we observe a similar behavior using RFG, as shown in Appendix A in Fig. 19. Part of the low momentum protons lost in INCL are emitted as clusters. Such clusters will not have the same detector signature of protons, as discussed in the following, and they do not fill the suppression at very large values of δα T , as can be seen in Fig. 10: indeed, there is still a sizable fraction of events, at low proton momenta, with full proton absorption and only a muon in the final state. There is no cluster production in NuWro.  SF and NuWro+INCL before and after FSI (top right, shape comparison). Distribution of difference in proton's cos(Θ) before and after FSI (∆ cos Θ = cos Θ af ter − cos Θ bef ore ) in NuWro SF and NuWro+INCL (bottom left, shape comparison). The "no FSI" channel is excluded. CCQE events with T2K neutrino energy flux are simulated.

B. Nuclear cluster production
In events with proton absorption multiple nucleons production in INCL, different particles can leave the nucleus, as shown in Fig. 11. As one can see, INCL features a notable production of deuterons, α-particles, 3 He isotopes and tritons. In this paper, we focus on studying production of these nuclear clusters. The production of clusters in nucleon scattering on nuclei is a well established evidence in nuclear physics, see for instance Ref. [95].
The production of clusters by FSI in neutrino interactions is a novelty of this study. While the presence of clusters can change the interpretation of various measurements in data (proton multiplicity, vertex activity, single transverse variables, etc.) the signature of such clusters is very detector-dependent, notably regarding granularity and threshold of measurable energy deposits. We provide here the fundamental information to characterize the observability and signature of such clusters in different detectors. In particular, the path length of cluster tracks, their mis-identification probabilities, evaluated with a simplified identification algorithm, and the energy released around the neutrino vertex are quantified. Figure 12 shows the length of the tracks produced by different particles as a function of their visible energy, as defined in Eq. 6. The sparse population of events at lower track lengths correspond to inelastic processes, which are included in the simulation but not in the characterization of visible energy, as explained in Sec. IV. In figure 13 we show the momentum distribution and the deposited ionization energy with Birks correction as a function of momentum. These distributions are the result of the simulation and the input to the algorithm described in Sec. IV. Table II shows the fraction of events where clusters travel more than 1 or 3 cm in the detector. While α and 3 He do not travel enough and will mostly contribute to vertex activity (energy deposited around the vertex), deuteron and tritium leave a visible track in a sizable fraction of events and, if searched for, can be quite cleanly identified. The primary source of misidentification are secondary interactions through inelastic processes, which typically happen at the end of the cluster track. In such inelastic interactions the cluster looses a large amount of energy and breaks the nucleus. The observed energy released by ionization along the cluster track, before the end of the track, does not include the energy released by secondary interactions and it is therefore smaller than what is expected for a cluster of that energy. As a consequence, particles could be misidentified with less ionizing ones. Therefore, tritium could be misidentified (∼10% of events) as proton or deuteron with similar probability, deuteron could be misidentified (∼20% of events) as protons, while protons cannot be misidentified as nuclear clusters. 3 He and α rarely travel more than 1 cm and when they do, they can be easily identified due to their very large energy deposit. The characterization of the inelastic interactions, the particles produced and their   The particles that are below tracking threshold still contribute to the total energy deposited around the vertex. To compute the contribution of clusters to the vertex activity, the visible energy deposited by ionization in a sphere of 1 (3) cm radius around the vertex is calculated, as shown in Fig. 14

VI. COMPARISON TO DATA
The prediction of NuWro and INCL are compared to STV measurements (δα T and δp T ) of T2K [88] and Minerνa [89] in Figs. 16 and 17. The contribution from non-QE interactions is simulated with NuWro and added on top of INCL predictions. Simulated events are selected to cope with the acceptance of the detectors. Such acceptance cuts, notably requiring large enough proton momentum to be reconstructed, tend to suppress the difference between the models. We look forward to measurements with lower tracking threshold, as expected for instance with ND280 upgrade [7] and the detectors of the SBN program [96]. The different shape in δα T is washed out, still the lower overall normalization of INCL, due to proton absorption, is visible and tends to slightly improve the comparison to data. The impact of nuclear clusters in the distributions is shown assuming perfect identification and same acceptance as protons. With present momentum threshold, nuclear clusters tend to distribute at high δα T , similarly to protons affected by FSI. Given the relatively high momentum threshold for proton tracking, the different proton momentum distributions after FSI in the two cascade models has a subleading effect.

VII. CONCLUSIONS
We have compared the simulation of the final-state interactions between the NuWro and INCL cascade models in CCQE events. We have characterized the produced particles and the change of kinematics of the leading proton induced by FSI. We have quantified the impact on observables like STV and vertex activity. We have also compared the results of the simulations with the existing experimental data provided by T2K and MINERνA  collaborations.
An essential novelty of this study is the simulation of cluster production by INCL in FSI of neutrino interactions. Such clusters can be detected as highly ionizing tracks or, if below the tracking threshold, can contribute to the energy deposited around the vertex (vertex activity). Nuclear clusters behave differently than protons from many points of view. Since nuclear clusters release more energy by ionization, the detector tracking threshold tends to require the clusters to have larger momentum than protons in order to be reconstructed. When identified as tracks, the particle momentum could be measured by its curvature, and the corresponding kinetic energy depends on the actual particle mass. Generally, the kinetic energy of a cluster leaving the nucleus is different from the kinetic energy of the initial proton producing the cluster. Moreover, for a complete recollection of all the energy of the particles leaving the nucleus, the correct reconstruction of their secondary interactions is needed. While we leave the detailed characterization of secondary interactions for future work, we have considered the fraction of energy lost in secondary interactions and their probability, which depends on the nature of the particles leaving the nucleus. Notably, emitted nuclear clusters have different probability of secondary interactions with respect to protons. For all these reasons, the accurate modeling and identification of particles in the final state, including nuclear clusters, are crucial to estimate correctly all the energy leaving the nucleus and the FSI corrections to it. Thus, the modeling of cluster production by FSI is relevant for correct interpretation of data to constrain nuclear models and it may affect the calorimetric method of total neutrino energy reconstruction and how it is shared between the leptonic and hadronic part of the final state. The precise neutrino energy reconstruction is a vital step for the new generation of oscillation experiments. The calorimetric method of the neutrino energy reconstruction, which takes into account all emitted particles after the neutrino interaction, allows for better neutrino energy resolution but also demands better modeling of the hadronic compound of the neutrino-nucleus interactions. Understanding the production and kinematical properties of nuclear clusters is crucial for improving the existing neutrino-nucleus interaction models and controlling the corresponding systematic uncertainties.
Regarding the effects of the final-state interactions on the leading proton in neutrino-nucleus interactions, the nuclear model of INCL features much smaller transparency than the NuWro cascade. Transparency is a measure that allows estimating the FSI strength in a given model. Therefore, such results originate in a characteristic combination of the nucleon-nucleon cross sections, Pauli blocking, and the specifics of the used nuclear model. In particular, INCL FSI simulation features a significant fraction of events without a proton in the final state, especially while propagating low momentum protons from the primary vertex. The INCL model also tends to re-absorb other particles produced during the FSI cascade: events with low momentum protons are mostly accompanied in NuWro with additional nucleons but not in INCL.
The accurate modeling of the rate of events without protons in the final state directly affects the correct interpretation of experimental data in long-baseline neutrino oscillation experiments. Regarding specific observables, the depletion of low momentum protons causes suppression of large values of δα T . Unfortunately, the proton tracking threshold in the detector induces a similar effect: as we show in the comparison to T2K and Minerva measurements, the available data have very limited sensitivity to the difference between the two FSI models due to the proton momentum threshold of present detectors. These results demonstrate the importance of measuring low momentum protons to characterize their FSI. Currently running and upcoming experiments are planning to perform such measurements by providing low tracking thresholds and precise calorimetry [97].
Summarizing, the study presented here is a progress along the road of a complete and precise estimation of uncertainties due to the modeling of final-state interactions effects of neutrino-scattering simulation. It is also a useful example of the modularity of Monte Carlo generators and a step forward toward improving the simulation of FSI in neutrino-nucleus interactions. The INCL FSI study presented in this paper relies on an approximated simulation where the initial state is taken from NuWro SF. As discussed in Sec. III, biases can be induced in the procedure of merging NuWro SF with the INCL nuclear model. The study is repeated using relativistic Fermi gas (RFG) in NuWro to model the initial state, in order to demonstrate that the conclusions on the characterization of FSI in INCL and the differences between NuWro and INCL FSI are robust against assumptions on the initial nuclear state. It is well known that RFG model is unable to reproduce the measured distributions of proton kinematics and STV [88,89]. RFG is here used only as a diagnostic tool to test the robustness of the FSI results. In this case, the momentum and position of the nucleon is provided by NuWro RFG and the INCL nuclear state is built around that nucleon.
We are also interested to test a more coherent simulation where both the initial state and the final-state nuclear effects are simulated with the INCL model. To this aim, the NuWro RFG simulation is reweighted to represent the INCL distribution of nucleons momenta and position shown in Fig. 18. The binding energy (for each nucleon position and momentum) and the total cross section are kept as in the original NuWro RFG model. In Tab. III, the fractions of different FSI channels are similar to those observed with SF simulation in Tab. I. The fractions change slightly due to the different momentum distribution of proton before FSI but similar trends are observable.
In Fig 19, δα T distributions are shown. RFG features an enhancement of small δα T values due to the large fraction of protons with large momentum, just below the Fermi limits. Still the characteristic suppression of such region induced by FSI in INCL is visible.