A quantum Monte Carlo based approach to intranuclear cascades

We propose a novel approach to intranuclear cascades which takes as input quantum MonteCarlo nuclear configurations and uses a semi-classical, impact-parameter based algorithm to modelthe propagation of protons and neutrons in the nuclear medium. We successfully compare oursimulations to available proton-carbon scattering data and nuclear-transparency measurements. Byanalyzing the dependence of the simulated observables upon the ingredients entering our intranuclearcascade algorithm, we provide a quantitative understanding of their impact. Particular emphasisis devoted to the role played by nuclear correlations, the Pauli exclusion principle, and interactionprobability distributions.


I. INTRODUCTION
The propagation of nucleons through the nuclear medium is an important aspect of nuclear reactions, heavy-ion collisions, and astrophysical environments. It is also crucial in the analysis of electron-nucleus scattering experiments (see e.g. Refs. [1][2][3][4][5][6][7][8][9][10]), neutrino oscillation experiments [11][12][13][14], and dark matter searches [15]. For example, scattering of neutrinos on nuclei can produce a number of outgoing hadrons. The multiplicity of the recoiling hadrons can indicate statistically if the incoming particle was more likely to be a neutrino or anti-neutrino [16]; their momenta are correlated with the original neutrino energy and direction, which can be used to gather additional information, for instance, on the leptonic CP phase from the sub-GeV atmospheric neutrino samples [17] or from searches for dark matter annihilation from the Sun [18,19].
The quantitative understanding of nucleons propagation in the nuclear medium would in principle require a fully quantum-mechanical description of the hadronic final state. Due to its tremendous difficulty, this problem has been tackled by introducing different approximations. The seminal papers by Serber and Metropolis et. al. [20][21][22][23] laid the foundations for the use of Monte Carlo techniques in semi-classical intranuclear cascades (INC) that assume classical propagation between consecutive scatterings. The latter are modeled using free-space elementary cross sections whose final state is modified to account for the Pauli principle. The results obtained by these first implementations of INC agree at least qualitatively with experimental data on nuclear transparencies and the frequency and angular distribution of emitted fast protons [21].
Experiments in heavy-ion reactions have spurred theoretical efforts to describe the dynamical evolution of nucleus-nucleus collisions using transport methods [24][25][26][27][28] based on the Kadanoff-Baym equations [29,30]. The impossibility of solving these equations exactly on real-world computers requires the introduction of approximations, such as the gradient expansion leading to the Boltzmann-Uehling-Uhlenback equations (BUU). For instance, the Giessen Boltzmann-Uehling-Uhlenback (GiBUU) model is based on this truncated set of semiclassical kinetic equations, which describe the dynamics of the hadronic system explicitly in phase space and in time [31,32]. It can be difficult to estimate systematic effects coming from the truncation of the transport equations. Although initially developed to simulate heavy-ion collisions, GiBUU has been extended to the description of lepton and photon scattering on nuclei [33,34].
Over the years, several studies have been devoted to improving the accuracy and extending the predictive power of INC models, with a focus on the analysis of nuclear spallation processes, where a hadronic probe with energy from a few tens of MeV to a few hundred MeV strikes a nucleus [35][36][37][38]. In contrast to the standard Glauber approach [39], where the so-called frozen approximation is utilized, INC simulations explicitly account for the motion of the background particles or scattering centers [40]. The validity of the semi-classical propagation assumed in the INC has yet to be fully assessed. However, genuine quantum-mechanical effects can either be safely neglected because they are expected to play a minor role-as for coherent scattering in nuclear transparency calculations-or they can be effectively parametrized by means of trajectory deflections and nuclear collective excitations, as in the analysis of nuclear spallation processes [41].
Important examples of state-of-the-art INC codes used in hadron-and lepton-nucleus scattering analyses are the Liège INC [42], NEUT [43], nuance [44], PEANUT (used within FLUKA) [45,46], NuWro [47], and GENIE [48] programs. Despite some differences in technical aspects and degree of sophistication, all the above INC models use as input elementary cross sections and mean-field properties of nuclei, such as single-nucleon densities. A notable exception is NuWro, which has recently been ex-tended to incorporate effective nucleon-nucleon correlations [49].
In this work, we propose a novel cascade model that employs nuclear configurations obtained from quantum Monte Carlo (QMC) calculations, which retains all correlation effects. The probability that the propagating nucleon scatters with a background particle is modeled using two different weight functions that depend on the impact parameter of the two nucleons and their cross section. The nucleon-nucleon cross sections used in the simulation are taken to be the same as in vacuum; in-medium effects are partially accounted for by imposing the Pauli exclusion principle and an effective nuclear binding.
Our semi-classical approach allows for the calculation of exclusive quantities in nuclear scattering, such as the fully differential phase space and the number of outgoing nucleons. To test the validity of our model, we compare our simulations with available data on proton-nucleus scattering cross sections and nuclear transparencies.
On the technical side, the INC code developed here is publicly available at https://github.com/jxi24/ IntranuclearCascade. The front-end is written in Python 3, while computationally intensive code is written in C++, using pybind11 to generate the interface code. The source code is organized to be portable and modular, making it easy to use, extend, and improve. All validations are available in the source code and can be easily performed by the front-end user.

II. THE PHYSICS OF THE IMPACT PARAMETER BASED INTRANUCLEAR CASCADE MODEL
We now discuss in detail the physical content of our INC and how several ingredients were incorporated in the model. The INC starts by generating a nuclear configuration, which describes the spatial distribution of protons and neutrons inside the nucleus. A large number of these configurations are computed from quantum Monte Carlo (QMC) methods [50] that use as input the highlyrealistic Argonne v 18 (AV18) [51] plus Illinois 7 (IL7) [52] Hamiltonian -more details can be found in Sec. II A. All particles in a nuclear configuration are initially labeled as background particles. After randomly picking a nuclear configuration, a random nucleon is selected, struck, and labeled as a propagating nucleon. Propagating nucleons are assumed to be point-like, on-shell particles moving at some velocity given by their four-momentum. The spatial coordinates of background nucleons are kept fixed until they interact with a propagating particle. In that case, their momentum is sampled from either a local or global Fermi gas distribution.
The system is evolved in time steps, which is a parameter of the cascade model, and should be chosen small enough to simulate the cascade accurately but large enough to run in a reasonable amount of time. At each time step, the following procedure is performed for each propagating nucleon. First, it is checked if there are any background particles within the volume perpendicular to the line segment defined by the initial to the final position of a propagating particle; this is schematically shown in Fig. 1. The impact parameter is calculated for all such nucleons, and this list is sorted from closest to furthest. Then an accept-reject step is performed on this list, until either an interaction takes place or the end of the list is reached. This test determines if an interaction occurred according to a probability distribution. If an interaction occurred, the phase space is generated using a fully differential nucleon-nucleon cross section. The Pauli exclusion principle is approximately taken into account by comparing the magnitude of the momentum of the final state particles with the Fermi momentum k F (the user can choose between either the global or the local values). The interaction takes place only if both momenta are above k F . Otherwise, the interaction does not happen, and the propagating particle keeps its original four momentum.
If the interaction took place, the outgoing particles are both treated as propagating particles, and a formation zone is set for each of them [53,54] (see also Ref. [47]). Note that, as a distinctive feature of our INC model, the interaction is finite ranged and the outgoing nucleons do not have to be at the same position in space. Finally, a propagating particle that reaches a sufficiently large radius may exit the nucleus if its kinetic energy overcomes the effective nuclear binding. If there is insufficient energy to escape, the particle is relabeled as a background particle. Otherwise, this effective binding energy is subtracted from the particle's energy, and then the particle is labeled as a final particle and stops propagating. The INC stops when there are no more propagating particles in the nucleus.
The structure of the algorithm is summarized in Fig. 2.
In what follows, we provide the details and expressions that are used in our impact parameter based INC.

A. Nuclear configuration
Nuclei are complicated many-body systems of fermions, whose structure and dynamics emerge from individual interactions among the constituent protons and neutrons. To a remarkably large extent, the latter can be modeled by the non-relativistic Hamiltonian where p i denotes the momentum of the i-th nucleon with mass m N , while v ij and V ijk are the nucleon-nucleon (NN ) and three-nucleon (3N ) potentials, respectively. We employ the highly-realistic Argonne v 18 (AV18) NN interaction [51] that includes spin, isospin, tensor, spinorbit and quadratic momentum-dependent terms as well as isospin-symmetry-breaking corrections and reproduces the Nijmegen nucleon-nucleon database with a χ 2 /datum 1. For the 3N interaction, we use the Illinois-7 (IL7) potential [52], which consists of the dominant standard Fujita-Miyazawa two-pion exchange and smaller multipion-exchange components resulting from the excitation of intermediate ∆ resonances. The IL7 potential also contains phenomenological isospin-dependent central terms.
The parameters characterizing this three-body potential have been determined by fitting the low-lying spectra of nuclei in the mass range A=3-10. The overall AV18+IL7 Hamiltonian then leads to predictions of ≈ 100 groundand excited-state energies up to A=12, including the 12 C ground-and Hoyle-state energies, in good agreement with the corresponding empirical values [50]. Over the last decades significant progress has been made in the development of ab initio nuclear methods, which solve the many-body Schrödinger equation associated with the Hamiltonian of Eq. (1) with controlled approximations [55][56][57][58]. For light nuclei, quantum Monte Carlo (QMC) and, in particular, Green's function Monte Carlo (GFMC) methods have been exploited to carry out calculations of nuclear properties, based on realis-tic NN and 3N potentials, and consistent one-and twobody meson-exchange currents [50]. GFMC begins with the construction of a trial wave function Ψ T that is a symmetrized product of two-and three-body correlation operators acting on an antisymmetric A-body singleparticle wave function that has the proper quantum numbers for the state of interest. The variational parameters in Ψ T are found by minimizing the energy expectation value where E 0 is the true ground-state energy of the system. The calculation of E T requires the numerical solution of a multidimensional integral that is carried out employing standard Metropolis Monte Carlo sampling in configuration space. GFMC then projects out the lowest eigenstate Ψ 0 of the given quantum numbers starting from Ψ T by performing a propagation in imaginary time τ The propagation In addition to binding energies the GFMC provides detailed information on the distribution of nucleons in a nucleus in both coordinate and momentum space, which are interesting in multiple experimental settings. For example, the mixed-estimate of the single-nucleon density is calculated as is the neutron or proton projector operator; and, ρ N integrates to the number of protons or neutrons. The two-body density distribution, yielding the probability of finding two nucleons with separation r, is defined as The positions of the constituents protons and neutrons utilized in the nuclear cascade algorithm are sampled from 36000 GFMC configurations. We employ the socalled constrained-path approximation [59] to make sure that their Monte Carlo weights remain positive, thereby facilitating their usage in the cascade algorithm. As a consequence, the single-proton distribution displayed by the blue solid circles of Fig. 3 is slightly different from the results reported in Ref. [60], which have been obtained performing fully unconstrained imaginary-time propagations. Since we neglect the charge-symmetry breaking terms in the Hamiltonian, and since 12 C is isospin symmetric, the single-neutron distribution is identical to that of the proton. For benchmark purposes, we also sample 36000 meanfield (MF) configurations from the single-proton distribution. The corresponding single-proton densities coincide by construction with the GFMC one, as shown in Fig. 3. However, the differences between GFMC and MF configurations become apparent when comparing the corresponding two-body density distributions represented in Fig. 4. The short-range repulsive core of the N N interaction prevents two nucleons from being close to each other. As a consequence, the pp and np GFMC density distributions are small at short separation distances. Furthermore, the difference between the GFMC pp and np density distributions around r = 1 fm can be attributed to the strong tensor correlations induced by the one-pionexchange part of the N N interaction, which is further enhanced by the two-pion-exchange part of the 3N potential. Note that the short-range behavior of ρ N N , which is largely nucleus independent, does depend strongly on the N N interaction model [61]. On the other hand, the MF ones do not exhibit this rich behavior as the correlations among nucleons are entirely disregarded.

B. Nucleon momentum distribution
As mentioned above, when a nucleon is struck, its momentum is obtained assuming either a local or global Fermi gas distribution. In the case of the local Fermi gas, the magnitude of the three-momentum is randomly sam- is the Fermi Momentum defined in terms of the single nucleon density k N F (r) = (ρ N (r)3π 3 ) 1/3 and N = p, n. In the case of the global Fermi gas, the momentum is determined in the same way, but k N F is position independent. The local Fermi gas model is known to provide a more realistic nucleon momentum distribution for finite nuclei than the global Fermi gas. For this reason, although both models are implemented in our code, we only present results based on the local Fermi gas predictions. In the future, we plan to include more accurate nucleon momentum distribution, based on state-of-the-art many-body calculations that properly account for nuclear correlations.

C. Nucleon-nucleon interaction algorithm
To check if an interaction between nucleons occurs, an accept-reject test is performed on the closest nucleon according to a probability distribution P (b) (see e.g. Ref. [62] for similar considerations) where b is the impact parameter. We impose two conditions on this probability, where the cross section σ depends on the incoming particle content and the center-of-mass energy, which is sampled from the nuclear configuration. The second condition ensures that the mean free path of a nucleon traveling in a medium of uniform density is λ mfp = 1/σρ, whereρ is the number density. Two implementations of P (b) have been studied here. The first we dub the cylinder interaction probability, where Θ(x) = 1 if x ≥ 0, else Θ(x) = 0. This probability mimics a more classical, billiard ball like system, where each billiard ball has a radius ≈ σ/π. The second implementation is the Gaussian interaction probability which is inspired by the work of Ref. [62]. Both P cyl and P Gau satisfy the conditions in Eq. (6). We use the nucleon-nucleon cross sections from the SAID database [63] obtained using GEANT4 [64], or from the NASA parametrization [65].

D. Phase space, Pauli blocking and after-interaction
If an interaction occurred, the phase space of the outgoing particles is generated using fully differential nucleon-nucleon cross sections. Note that, at the moment, we only include protons and neutrons in our INC model. Pauli blocking enforces Fermi-Dirac statistics for the nucleons and amounts to testing whether their finalstate momenta are above the Fermi momentum. Two different models of the Pauli exclusion principle have been approximately implemented. The global and local Pauli blocking routines essentially forbid a scattering if the momentum of any of the final state particles is below the average Fermi momentum (for the global Fermi gas model) or the local Fermi momentum (for the local Fermi gas model), respectively. We emphasize again that, although we have implemented the global Fermi gas model, we do not report any results using it.
If the interaction took place, the outgoing particles are both treated as propagating particles, and a formation zone is set for each of them [53,54] (see also Ref. [47]). The formation zone is a length or period of time in which a particle does not interact with any nucleons. This models the coherence and interference of interactions in quantum mechanics. The formation zone is given by where m N = 938 MeV is the nucleon mass, E is the energy of the outgoing nucleon, and p and p are the four momenta of the incoming and outgoing nucleon respectively.

E. Exiting the nucleus
When the radial position of a propagating particle is larger than a certain distance, much larger than the nuclear radius, a test is performed to check if the particle has enough energy to escape the nuclear potential. If its kinetic energy is larger than the the nuclear potential barrier, then the particle escapes the nucleus and is labeled as a final state particle. The momentum of the particle is modified to include the effective nuclear binding. Final state particles do not propagate and are essentially the input that should be given to a detector simulation.

III. RESULTS AND MODEL VALIDATION
In this Section, we present the validation tests we performed and the results we obtained within our INC model. For comparison purposes, we also implement a version of the nucleon-nucleon interaction algorithm that we dub mean free path (MFP). This approach is routinely used in event generators. Within this algorithm, the system is not evolved in time steps but rather in constant position steps that we indicate as λ max . At each step of the loop the mean free path of the propagating particle readsλ where N = n, p refers to the isospin of the propagating particle and r its distance from the centre of the nucleus. Note that in the limit of constant density, ρ p(n) (r) =ρ, we would recoverλ = λ mfp . (see also Sec. III A below). The probability that the struck particle traveled a distance λ without interacting can be written as P (λ) = e −λ/λ . We can then sample λ = −λ·ln(rand[0, 1]) and say that an interaction took place if λ < λ max . Note that λ max has to be chosen small enough in order to satisfy the assumption of a constant density. The simplest validation consists of calculating the mean free path of a nucleon travelling through a medium of uniformly distributed nucleons, with fixed interaction cross sections. Our goal is to obtain where σ is the nucleon-nucleon cross section, and ρ is the number density of nucleons. We simulate a background of static, randomly and uniformly distributed protons and neutrons at rest. For an event, we insert into the medium a test particle with a fixed energy and propagate it through the medium. When the first interaction happens, we record the length the test particle has traveled. Figure 5 shows the distribution of lengths traveled for a 500 MeV test proton in a 0.16 nucleons/fm 3 background density of nucleons, assuming a fixed interaction cross section of 50 mb, and using the Gaussian interaction probability. The expected distributions, according to Eq. (11), are also shown for the fitted (black line) and the expected (orange line) values of the mean free path. The nuclear density and interaction cross section were chosen arbitrarily. Adjusting these values does not change the agreement between the expected mean free path and the calculated mean free path. While not shown here, using the cylinder interaction probability does not change the results either. As we can see, our code correctly reproduces the expected behavior for the mean free path, allowing us to proceed to more complex tests of our INC.

B. Proton-carbon Scattering Data
Reproducing the proton-nucleus cross section measurements is an important test of the accuracy of the INC model. Proton-nucleus scattering probes the nucleonnucleon cross section which is typically divided into two pieces, the reaction and the elastic cross sections, In the elastic part, no energy is transferred into nuclear excitation and the nucleus remains unbroken, that is n + A → n+A. The reaction cross section includes transition to nuclear excited states, n + A → n + A * , as well as inelastic reactions n + A → X. Several experiments have been carried out to determine the total reaction cross section, see for example Refs. [66][67][68][69][70][71]. The latter is typically obtained by measuring the total cross section from the change in intensity of a calibrated proton beam traversing a carbon target and then subtracting the calculated elastic cross section.
We compute σ R neglecting Coulomb interactions, as they are expected to contribute mostly to σ el . We obtain the proton-carbon scattering cross section by the following simulation (with a different setup from the proposed algorithm of Fig. 2). We define a beam of protons with energy E, uniformly distributed over an area A (orthogonal to the proton momenta). Note that A πR 2 , where R is the radius of the carbon nucleus. The carbon nucleus is situated in the center of the beam. We propagate each proton in time and check for scattering at each step. The Monte Carlo reaction cross section is then defined as the area of the beam times the fraction of scattered events, namely, This is not exactly the experimentally measured reaction cross section. Angular and/or momentum acceptances for the attenuated beam are finite, and we do not include these effects in our calculation. Nevertheless, we do not expect such effects to change our results significantly, and thus σ MC should be a good approximation of the reaction cross section. Moreover, imposing Pauli blocking on both outgoing nucleons will effectively suppress the contribution of elastic transitions. The two panels of Fig. 6 display the proton-carbon scattering cross sections as a function of the proton kinetic energy. In the upper panel our Monte Carlo simulations are compared with experimental data in the entire energy region in which data are available [71]  Note that the interaction model MFP does not use any configuration, but rather the local density. The solid lines have been obtained using the nucleonnucleon cross sections from the SAID database [63], obtained using GEANT4 [64], in which only the elastic contribution is retained. The dashed lines used the NASA parameterization [65], which includes inelasticities. The inelastic contribution leads to an enhancement of the cross section which is necessary to reproduce the data correctly for intermediate to large proton kinetic energies. The NASA parameterization allows us to approximate the effects of including pions in the cascade. Results obtained with these two nucleon-nucleon cross sections are indicated in Fig. 6 by "El" for elastic cross section and "Tot" for the total cross section. As expected, the p-carbon cross section obtained using only the elastic nucleon-nucleon cross section is consistently lower than the one obtained using the total nucleon-nucleon cross section. The effect is large for T p > 50 MeV, where pion production becomes relevant. In a future work, we will explicitly include pion degrees of freedom in our INC as these are known to play a crucial role in this region. Additionally, the present results neglect modifications to the nucleon-nucleon cross section arising from the nuclear medium itself. Such modifications have been considered in the literature [72]. We also leave these considerations for a future work.
The dependence of the results on the functional form of the interaction probability and on the nuclear model adopted to generate nuclear configurations has also been investigated. The blue and red lines are obtained from the cylinder probability of Eq. (7) and sampling the initial nucleon configurations from the quantum Monte Carlo (QMC) and mean field (MF) distributions, respectively. To obtain the orange and green curves, QMC and MF configurations were used together with the Gaussian probability of Eq. (8). We observe that in both cases, the QMC and MF curves are almost superimposed, indicating that this observable does not depend strongly on correlation effects among the nucleons. A more detailed discussion of their impact on particle propagation distances will be discussed in Sec. III D.
The purple lines refer to mean free path (MFP) calculations which present conceptual differences with respect to the cylinder and Gaussian cases and yield large discrepancies in the results. The predictions for the p-carbon cross section obtained from the MFP approach underestimate the experimental data and are significantly lower with respect to the other curves obtained with more refined techniques.
In the lower panel of Fig. 6 we focus on the comparison between the different MC simulations and the data in the low energy region, i.e., proton kinetic energies below 200 MeV. The results obtained using the Gaussian and cylinder algorithm are in fair agreement with the experimental data: the curves display the correct behaviour although the position of the peak is not exactly reproduced. We note that Pauli blocking plays a fundamen-tal role for these kinematics, significantly quenching the results. More sophisticated techniques aimed at consistently orthogonalizing the nuclear wave function in the final state are discussed in Ref. [73]. Their implementation in our INC will be the subject of a future work.
The functional forms adopted to determine the probability distribution (Gaussian vs. cylinder) yield predictions which differ by ∼ 15% in the low energy region (compare, for instance, the green and red lines in Fig. 6). For intermediate proton kinetic energies T p the curves converge to the same result and reproduce the data. At low T p , the predicted cross sections for the MFP method are ∼ 50% smaller than those of the Gaussian and cylinder methods. As these differences are also present in nuclear transparency, we will discuss their origins in the next section.

C. Nuclear transparency
The nuclear transparency yields the average probability that a struck nucleon (either a proton or a neutron) leaves the nucleus without interacting with the spectator particles. Measurements of the nuclear transparency to high energy protons in quasielastic e + A → e + p + A scattering have been carried out by a number of experiments [4,6,7,[74][75][76]. These measurements are performed with fixed kinematics (i.e., for a fixed incoming beam energy and scattering angle) and measure the outgoing electron and proton. The proton produced in the primary vertex can be absorbed or deflected while exiting the nuclear medium because of final state interactions with the remnant system, leading to a reduction of the measured e+A → e+p+A cross section. This reduction defines the nuclear transparency, which is given by the ratio of the observed events to events predicted in the Plane Wave Impulse Approximation (PWIA) To obtain N PWIA the initial proton is treated as a bound state with energy and momentum distribution described by a spectral function; while the final proton is a free particle state propagating as a plane wave. In our Monte Carlo simulation, the nuclear transparency has been obtained via the following procedure. We randomly sample a nucleon, its position being generated according to either the QMC or MF distributions. We give the nucleon a kick by assigning it a given kinetic energy T p and three-momentum p, and propagate it through the nuclear medium. The Monte Carlo transparency is then defined as where Pauli blocking has been implemented in the determination of N hits . Note that for a given initial and final   7: Carbon transparency as a function of the proton kinetic energy. The different curves indicate different approaches used as described in Fig. 6. The experimental data are taken from Refs. [4,6,7,[74][75][76] energy and scattering angle of the electron, one can unambiguously define the momentum q transferred to the target nucleus. The direction and the momentum of the nucleon in the final state has to be determined applying energy-and momentum-conservation relations and accounting for the Fermi motion of the struck nucleon in the initial state. It follows that defining the kinematics of the hadronic final state after the hard scattering depends on the nuclear model of choice. However, in the analysis of different experiments, the data are given as a function of the average nucleon momentum (and kinetic energy) given by p = q (T p = |q| 2 + m 2 N − m N ). In Fig. 7 we compare the nuclear transparency data from Refs. [4,74] to our predictions. The different lines are the same as for Fig. 6. We find an overall satisfactory agreement between the Gaussian and cylinder curves with the experimental data once inelastic effects are taken into consideration; this corresponds to the results using the NASA parametrization for the nucleonnucleon cross sections. For moderate to large values of the proton kinetic energy, pions play an important role in quenching the transparency. Moreover, the Gaussian and cylinder curves exhibit correct behavior consistent with the data also for small T p where the simplified MFP model described above fails. As in Fig. 6, we observe very small differences between the QMC and MF calculations. For low and intermediate kinetic energies, the transparency obtained from the MFP approach is much smaller than the corresponding results for the cylinder and Gaussian curves.
Finally, we discuss the origin of the discrepancies between the MFP and the cylinder algorithm with MF configurations for the p-carbon cross section and carbon transparency. Both approaches rely on the single-nucleon density distribution to sample the initial nucleon posi- The distance from the proton to the center of the nucleus is r 1 , and the propagation step is d . The radius of the cylinder is given by σ/π where σ is the interaction cross section between the proton and a background particle; d is also the height of the cylinder. Right panel: same as for the left one, but for a nucleon kicked inside the nucleus.
This follows what is done in the nuclear transparency event simulations.
tions (nuclear correlations are neglected) but use different definitions of the interaction probability. The left panel of Fig. 8 schematically shows one contribution to the p-carbon cross section in which the proton is at a distance r 1 larger than the nuclear radius. In the cylinder algorithm, the interaction probability is equal to one if a particle is present in the volume defined by: Both σ pp and σ np have a maximum for low proton momentum values. Hence, for low momenta, the probability of interaction could be non-vanishing even when the projectile proton is far from the center of the nucleus. On the other hand, within the MFP approach, if the probe is outside the nucleus then the approximation of a constant density ρ(r 1 ) = 0 within the volume V = d · σ yields a vanishing interaction probability. This different behaviour leads to a lower p-carbon cross section using the MFP approach, as observed in Fig. 6. When computing the nuclear transparency we kick a nucleon which is located inside the nucleus as displayed in the right panel of Fig. 8. In this case, assuming a constant density is more likely to overestimate the interaction probability, especially for low momenta where the cross section is larger. This observation is consistent with Fig. 7 where the MFP curves predict a larger number of interactions, and therefore a lower nuclear transparency, for small T p .

D. Correlation effects
The role played by nuclear correlations in final state interactions of the recoiling nucleon has been investigated FIG. 9: The four panels corresponds to histograms of the distance traveled by a struck particle before the first interaction takes place for different values of the interaction cross section. The results in blue and red correspond to MF and QMC initial nucleon configurations, respectively. For each of the panels we also report the fixed cross-section used, the total number of events generated, and the number of hits for each configuration.
in Refs. [72,[77][78][79][80]. As discussed in Ref. [81] the hit nucleon is surrounded by a short-distance correlation hole produced by both the Pauli principle and the repulsive nature of realistic nuclear interactions. Because of this correlation hole, the stuck nucleon is expected to freely propagate for ∼ 1 fm before interacting with any of the background particles. To test the validity of these observations in our INC model, in Fig. 9 we report the histograms of the distance traveled by a struck nucleon before its first interaction occurs-we stop the simulation afterwards-with each panel corresponding to a different value of the interaction cross section. In order to gauge the effect of nuclear correlations, the initial positions of the nucleons are sampled from either MF (blue) or QMC (red) configurations. A random nucleon inside the nucleus is recoiled and assigned a momentum of 200 MeV. Pauli Blocking has been neglected here to isolate the dependence of the results on the spatial distribution of the nucleons. We employ the cylinder algorithm and use a fixed cross section-which determines the cylinder base area-varying between 0.5 and 100 mb.
For σ = 0.5 and 10 mb, the volume spanned by the propagating particle is very small. The first and second panels of Fig. 9 clearly show the MF distribution peaking toward smaller distances than the QMC distribution. This difference primarily originates from the short-range repulsion of the AV18 potential that reduces the probability of finding two nucleons close to each other and allows the struck particle to propagate longer before interacting. This effect is more pronounced for cross sections below about 10 mb = 1 fm 2 since correlations affect nucleon configuration for inter-particle distances within 1 ∼ 2 fm, as can be seen in Fig. 4. On the other hand, larger cross sections yield larger cylinders. In this case, the propagating particle becomes less sensitive to the local distribution of nucleons and more sensitive to the integrated density in a larger volume, reducing the effect of correlations. For these larger cross sections, the MF and QMC event distributions follow the same trend, as can be seen in the lower panels of Fig. 9, corresponding to σ = 50 and 100 mb.
In each panel we also report the number of hits and the total number of events simulated. Obviously, the number of interactions increases for the larger cross sections. It is more interesting to notice that the total number of hits does not present a strong dependence on configurations with (QMC) and without (MF) correlations. This explains why sampling from these two different set of configurations give comparable results for the p-carbon cross section and carbon transparency, displayed in Figs. 6 and 7. These observables are indeed more sensitive to the occurrence of any interaction rather than to its location inside the nucleus. This is especially true for the p-carbon cross section, since the attenuation of a proton beam is sensitive to the total nuclear density rather than its differential profile (compare, e.g., the solid blue and red lines in Fig. 6).
An alternative way to include correlation effects in INC models, based on two-body density distributions, has been recently proposed by the NuWro collaboration and leads to a larger enhancement of the nuclear transparency than we find [49]. This difference can partly be ascribed to correlation terms between two spectator nucleons [79] that are automatically incorporated in our QMC configurations, but cannot be modeled by the twobody correlation functions alone. In addition, as discussed in Sec. III D, when the nucleon-nucleon cross section is larger than about 20 mb, the base of the propagating cylinder in our model covers the correlation hole and reduces the importance of correlations. This effect is absent when only local properties of two-body distribution functions are considered.
Although nuclear correlations in final state interactions do not play a prominent role in the observables that we have considered, correlations may be important in other experimentally relevant quantities or in high-precision studies. These observables include proton multiplicity and the distribution of outgoing protons' direction and energy. We plan to continue exploring these questions in future work.

IV. CONCLUSIONS
We have presented a novel INC model that takes as input realistic QMC nuclear configurations to sample the distributions of protons and neutrons inside the nucleus. Either a cylinder or a Gaussian distribution is used to define the interaction probability of the propagating nucleon as a function of the impact parameter and the nucleon-nucleon cross sections. We have considered the elastic and the total cross sections, corresponding to the SAID database [63] and to the NASA parametrization [65], respectively.
To validate the INC we have compared our simulations with experimental data for p-carbon scattering and carbon transparency. For both observables, we find minimal difference between the cylinder and Gaussian distributions, and both correctly reproduce the trend of the data. We gauge the role of nuclear correlations in the propa-gation of the recoiling nucleons by utilizing nuclear configurations computed within quantum Monte Carlo and comparing with mean field distributions. The distance traveled by a struck particle before the first interaction occurs, for sufficiently small value of the nucleon-nucleon cross section, is significantly longer when nuclear correlations are accounted for. This is due to the short-range repulsion characterizing the AV18 potential that gives rise to a correlation hole surrounding the propagating nucleon, thereby reducing the probability of finding two nucleons close to each other. However, the total number of hits is only mildly affected by the presence of nuclear correlations. As a consequence, they mildly affect our calculations of p-carbon scattering and carbon transparency. The results presented for the latter quantity by the NuWro collaboration indicate a larger importance of nuclear correlations, whose implementation is based on the two-body density distributions [49]. Nucleon-nucleon cross sections are an important ingredient in INC models. Since the scattering between the two nucleons takes place in the nuclear medium, the free nucleon-nucleon cross sections have to be corrected. In our simulation, we do so by approximately implementing the Pauli exclusion principle: in the aftermath of each collision, the magnitude of final nucleons' momenta have to be larger than the Fermi momentum. This constraint reduces the effective cross section and considerably improves the agreement between our simulations and experimental data. As a next step, along the line of Ref. [72], we plan to replace the bare nucleon mass with an effective one that depends on the momentum and the nuclear density, leading to a more complete and consistent treatment of in-medium effects.
At energies larger than the pion-production threshold, inelastic contributions to the nucleon-nucleon cross section become relevant. Using the NASA parametrization, which includes these inelasticities, noticeably improves the agreement with p-carbon scattering and carbon transparency experimental data in the higher-energy region. Our INC model has yet to include pion propagation in the nuclear medium. Existing event generators have to rely on a number of semi-phenomenological approaches constrained by experimental data [33,[82][83][84][85], the validity of which strongly depends upon the energy of the propagating pion. A consistent implementation of pion propagation and absorption is an extremely challenging problem and we leave this aspect for future works.
The transparency measurements used here depend upon theoretical calculations based on the plane wave impulse approximation. To compare directly to the experimental data, we plan on implementing the primary interaction vertex using the spectral function formalism [86][87][88]. Additionally to further validate our model, we will carry out extensive comparison against available semiexclusive electron scattering data [89,90].  [63] and obtained using GEANT4 [64] and the NASA model is given from [65].