A Monte Carlo generator of nucleon configurations in complex nuclei including Nucleon-Nucleon correlations

We developed a Monte Carlo event generator for production of nucleon configurations in complex nuclei consistently including effects of Nucleon-Nucleon (NN) correlations. Our approach is based on the Metropolis search for configurations satisfying essential constraints imposed by short- and long-range NN correlations, guided by the findings of realistic calculations of one- and two-body densities for medium-heavy nuclei. The produced event generator can be used for Monte Carlo (MC) studies of pA and AA collisions. We perform several tests of consistency of the code and comparison with previous models, in the case of high energy proton-nucleus scattering on an event-by-event basis, using nucleus configurations produced by our code and Glauber multiple scattering theory both for the uncorrelated and the correlated configurations; fluctuations of the average number of collisions are shown to be affected considerably by the introduction of NN correlations in the target nucleus. We also use the generator to estimate maximal possible gluon nuclear shadowing in a simple geometric model.


Introduction
The structure of nuclei has been described for a long time by independent particle models in which the nucleus is treated as a collection of fermions freely moving about, this picture being suggested by the fact that nuclear matter is a dilute system. Nevertheless the short range structure of nuclei cannot be accurately described by this simplified picture, and essential aspects of nuclei such as high momentum component of the nuclear wave function are completely missed by the non-interacting model; configurations of nucleons at short separations are known as short range correlations (SRC) [1,2]. NN correlations have been recently unambiguously observed [3] in a series of dedicated experiments, in which high-momentum nucleons were knocked out by high energy probes and their correlated partner were observed with opposite momentum; the probability for a nucleon to belong to a SRC pair was measured to be about 20% in 12 C and 25% in in heavy nuclei; see Ref. [4] for a review. A theoretical description of SRC can be made in a satisfactory way using nonrelativistic Hamiltonians containing two-and three-body realistic potentials. It is currently possible to solve exactly the Schrödinger equation for light nuclei while, for A > 12, approximate methods must be considered, which have been shown to produce a reasonable approximation to the basic ground state quantities such as one-and two-body densities and momentum distributions. These studies demonstrated that NN correlations produce non-negligible effects in a variety of processes including those which are not specially probing SRC. For example, significant ef-fects are found for high energy total nucleon-nucleus cross sections [5], and it has thus been shown that they must be taken into account in the theoretical calculation of such processes. In this Letter we are addressing the challenging problem of including NN correlations in the description of high energy proton-nucleus and nucleus-nucleus collisions. It was shown in Ref. [6] that using correlated twobody densities for the analysis of fluctuations of the mean number of participating nucleons in proton-nucleus collisions significantly modifies the results, especially the tails of the distribution, and it was pointed out the importance of implementing NN correlations in the event generator for nucleon configurations in a way consistent with the single nucleon nuclear density and use it for description of the heavy ion collisions on an event-by-event basis. Simulations of such processes are commonly performed using configurations of nucleons as an input, consisting of nucleons spatial locations and isospins, generated assuming the independent particle model; however, inclusion of NN correlations can provide a more realistic description of initial states of nuclei and appears to be feasible with a good accuracy.
The Letter is organized as follows; in section II, details of NN correlations are discussed and the method adopted for generation of configurations is outlined; we describe the generation of nuclear configurations within a MC approach including central two-body correlations which are consistent with realistic calculations of one-and two-body density distributions. In section III we introduce the probability distribution functions P N for the N -th or-der, N = 1, ..., A, interaction of a projectile nucleon with a nucleus A, and we present results for fluctuations. In section IV we use our event generator to calculate maximal possible shadowing for gluons based on simple geometric considerations.

NN correlations
The aim of the present work is to develop a MC procedure for generating spatial nucleon configurations to be used as an input for simulations of collisions involving heavy nuclei which treats NN correlations in a realistic way. The commonly used approach to this problem is taking a given number of nucleons out of a Woods-Saxon density distribution, which describes the probability density function of a nucleon to be located at a distance r from the center of the system. As a result, the A nucleons are positioned independently from each other, such that all two-particle (and higher) correlations are completely ignored; this is justified by the assumption that inclusive quantities are not very much dependent on these correlations. Alternatively one puts nucleons in consequently one after another imposing the condition that the distance between nucleons should be larger than some minimal one. This procedure, however, results in a wrong single nucleon distribution of the nucleons and must be improved.
We argue that when considering event-by-event observables, NN correlations in the nuclear wave function are relevant. This extends observation of [6] that NN correlations significantly modify the variance of the distribution over the number of collisions. In the present work we analyze several different approximations for correlations. A full implementation of realistic short range correlations [7,8] in the wave function Ψ of a complex nucleus can be made using an independent particle model wave function Φ and a proper correlation operatorF = A i<jf (r ij ) as follows: f being a set of state-dependent correlation functionŝ ij obtained variationally from a nonrelativistic Hamiltonian with realistic two-and threebody interactions;Ô (n) ij is the spin-isospin operator appearing in the realistic two-body potential, usually written In the present implementation of correlations, we did not consider the full realistic set of correlations shown in Eq.(1), which would imply using the full |Ψ| 2 as a probability function for the Metropolis random search. This would be numerically very demanding, especially for heavy nuclei, as it would involve generating a full configuration of A nucleons. Hence as a first step we will implement only effective central correlations; the inclusion of spinand isospin-dependent configurations, which allows to account for most of the properties of one-and two-body densities described in [7,9,10] and references therein, will be discussed later on.
A major problem with using a system of uncorrelated nucleons for the event generators is the occurrence of hot spots: regions in space where two or more nucleons are allowed to overlap and produce unrealistically high local densities. A simple excluded volume model, in which nucleons are prevented from being closer than a fixed distance of about the nucleon size d = 1 f m does not account for the realistic picture of nuclei, since this approach artificially moves density from the center of the nucleus towards the periphery, distorting the nuclear profile and altering the surface properties. The first constraint in our generator is thus to preserve the one body density of the nucleus, defined as follows: We suggest a method for generating the spatial configurations of nuclei with the following essential features: i) the nucleons are distributed according to a given singleparticle distribution; ii) each configuration embodies the NN pair correlations. The resulting two-body densities are by far more realistic than the commonly used independent particle model ones, or even those of the excluded volume approximation. The outlined results can be achieved by first generating a large number of random nucleon positions in a cube with given density (0.17 nucleons/f m 3 ). Next, a Metropolis random search is performed, with a probability function given by the square of the wave function; at this stage, we impose the constraint of the Woods-Saxon probability distribution (Ref. [11]), choosing those nucleons which are relevant to the desired one-body density. The procedure comes to an end when exactly A nucleons satisfy the constraints at the same time.
It was shown in Refs. [7,9] that SRC have sizable effects on the single particle density (Eq.(2)) as well as on those quantities which depend depend on the two-body density of the system, defined as such as the two-body momentum distributions. In many theoretical studies, Eq. (3) is approximated as follows where C(r 1 , r 2 ) vanishes if NN and statistical correlations are disregarded. We can compare the realistic results [7,12] for the radial two-body density, defined as the integral of Eq.(4) over the pair center of mass, R: and the results from our MC generator for the same quantity. If we assume that C(r 1 , r 2 ) depends only on the relative coordinate r = |r 1 − r 2 | (which holds true exactly for infinite nuclear matter and which is a good approximation in finite nuclei, far enough from the edge of the nucleus), the pair distribution function C(r) can be obtained from the correlated ρ U (r) two-body radial densities: The quantities defined in Eq.(5) and Eq. (6) were calculated using the MC code using several approximations for the correlation function f (r) (see Eq. (1)). They are shown in Fig.1; the results of Ref. [7] are also shown. Our procedure was to start with a central correlation function f (r) taken as a step function, i.e. vanishing wave function for r < 1 f m; we then used a Gaussian f (r) = 1 − exp(−β r 2 ) in order to match better the realistic result. It is known that the realistic central correlations has a small overshooting over unity, and it is in general different from our Gaussian correlation; however, in realistic calculations of finite nuclei [7,12] the peak in the radial two-body density is mainly due to the existence of state-dependent, mainly tensor correlations, and three-body diagrams. For this reason we prefer not to introduce an unrealistic central correlation to reproduce such a peak and use a parameter β = 0.9 f m −2 in the Gaussian correlation which gives a healing distance similar to the one exhibited by the realistic approach. This amounts to calculate the wave function with a central correlation operator; the fully realistic implementation of the operator of Eq.(1), involving spin and isospin degrees of freedom, extremely computationally intensive, is thus beyond the aim of the present contribution and will be performed elsewhere. In the present work we randomly assigned isospin degrees of freedom for A nucleons, so that each generated configuration consists of A quads of numbers, containing three spatial coordi- nates and one, randomly assigned, isospin variable; the correlated nuclear configurations are downloadable from the PSU physics department web server at the address[13].

Probability distribution functions for hadronnucleus collisions
In order to illustrate the present implementation of correlations, we have focused on the 16 O, 40 Ca and 208 P b nuclei.
For a given configuration of nucleons, the probability of interaction with the i-th nucleon for an incoming projectile with impact parameter b is given by the corresponding probability of no interaction being 1 − P (b, b i ). The Γ function in Eq.(7) for high energy incident nucleons can be parameterized as where σ tot N N is the total cross section of NN scattering and the t dependence of the cross section dσ/dt ∝ exp(Bt), neglecting small corrections due to the real part of the amplitude. Let us define the probability of interaction with N nucleons as a function of the impact parameter as The inelastic cross section due to collisions with N nucleons is then given by and the total inelastic cross section is σ in N A = A i σ in i , and it can also be calculated as σ el N A and σ qe N A being the elastic and quasi-elastic cross section, respectively; this calculation is being performed with realistic wave functions (Ref. [14]) and a comparison with the present approach will be presented elsewhere.
We can evaluate the probabilities given by Eq. where √ s = 5500 GeV /c as well. We have checked the independent particle model results with the analytic approximation of Ref. [15], where the probability distributions are given by: where ρ(b, z), which neglects the finite radius of NN interaction. The comparison is made for the  (14), from Ref. [6], with C(r) extracted from the MC configurations and shown in Fig.1. 16 O nucleus in Fig.2, showing a more narrow distribution if the approximation of Eq.(12) is used; the same holds true for the other considered nuclei. The quantities N and N (N − 1) , can be written as follows [6]: which can be calculated in the present framework using Alternatively we can calculate N , N (N − 1) using Eqs.(13), (14) as As a result we can write the variance of the mean number of collisions as a function of the impact parameter b, as follows . (18) In order to check the accuracy of our method, D(b) was calculated both using Eqs.(13), (14) and Eq.(18) using analytical one-body densities and C(r) extracted from the MC calculation; we then used Eqs. (16), (17), calculated with the P N (b) functions obtained with the MC. The results are compared in Fig.4 for 16 O, 40 Ca and 208 P b; a small discrepancy is exhibited, between the calculation with C(r) (shown with symbols) and the corresponding Gaussian correlated result, but the overall agreement is satisfactory. The small difference between the results of the two methods, especially for lead, is to be ascribed to the fact that the quantity C(r) has been extracted from the MC configurations using Eq.(6), assuming the function C(r 1 , r 2 ) to be a function of the relative distance only, and not to depend on the particular region of the nucleus considered, while surface effects should be present on the level of few %. We have compared the performance of the Glauber type interaction with the approach of Ref. [16], dxP (x) not to deviate from unity for more than 5%.
where the probability of interaction in Eq.(9) was taken as a step function, which vanishes for nucleons sitting outside the cylinder centered at the given |b| from the center of the nucleus. We took for the cylinder the value which gives an area corresponding to the one given by where Γ is calculated with the Glauber parameters we have used in the present work. The results, as compared with our Glauber approach showed to differ substantially from both the uncorrelated and correlated MC results.

Lower limit on the parton nuclear shadowing
In this section we present application of our MC code to the study of maximal possible nuclear shadowing in nuclei. We use the well known observation that shadowing in the scattering off the deuteron cannot reduce the cross section to the value smaller than cross section of scattering off one nucleon. Basically it is due to one nucleon screening another one but not itself.
Similarly, it is natural to expect that in any dynamics of parton interactions the gluon (quark) density at a given impact parameter in a particular configuration cannot be less than the maximum of individual transverse gluon densities, g N (x, ρ) in the nucleons at given b. Here g N (x, ρ) = g N (x)F g (ρ) is the generalized diagonal gluon density in the nucleon. leading to In the limit of A → ∞ this leads to The onset of the limiting behavior depends on the transverse shape of the gluon GPD. For the same g N (x) a shape more peaked at small ρ will asymptotically lead to larger value of g A (x, b) min though the approach to the asymptotic value will require larger values of A as chances that there is a nucleon with ρ small enough that g N (x, ρ) is close to g N (x, 0) is smaller in this case. In Ref. [17], two parameterizations for the two-gluon form factor were discussed, fitted to the J/ψ photoproduction data ( [18]); they were taken in the forms of an exponential and a dipole. The corresponding transverse spatial distributions of gluon GPDs are thus where K 1 is the modified Bessel function, m 2 g = 0.6 GeV 2 for x ∼ 10 −4 and B g = 3.24/m 2 g . Using the configurations described in the previous sections, we have calculated for given impact parameter b the maximum transverse gluon density normalized to its peak value, i.e. the maximum value of the probability density P (x, b) of the quantity as a function of ρ j = b−b j , with b j the j-th nucleon transverse coordinate; here P (x, b) is normalized according to dxP (b, x) = 1, for given b. The results of the calculation for 208 P b are shown in Fig.5. It can be seen that the two densities of Eq.(23) produce very similar results.
We next apply these results to determine minimal value of the gluon shadowing at given b defined as the ratio of lower limit on g A (x, b) and its value in the impulse approximation, g A (x, b) = g N (x)T A (b); to this end, we take the ratio Results are presented for the two used models of the gluon GPDs of Eq.(23) in Fig.6 and Fig.7. For A → ∞ the limits using the two different fits differ by F g (0)/F (2) g (0) = 2/3.24 = 0.62; however, this limit is reached at extremely large A, as shown in Fig.8. This is due to a rather small radius of the transverse gluon density r tr g ≤ 0.5 f m and low nuclear density, leading to a small probability for more than three nucleons to significantly screen each other up to very large A. To investigate the onset of asymptotic, we evaluated the A dependence of the quantity x = dxP (x) at zero impact parameter for the two considered gluon GPDs for two values of m 2 g corresponding to x ∼ 0.01 and 0.6 GeV 2 corresponding to x ∼ 10 −4 . We modeled nuclei with very large A by generating random nucleons with constant density ρ 0 = 0.17 f m −3 in a cylinder centered at b = 0, with height 2R A , where R A = 1.25A 1/3 ; the results, shown in Fig.8, are independent from the radius of the cylinder, provided it is larger than about 0.8 f m. It can be seen that, the GPDs being very peaked at ρ = 0, the increase of x is very slow and the maximum value is not reached even with A as large as 10 5 . Onset of asymptotic is somewhat faster for smaller m g due to a smoother behavior of GP Ds at ρ ≃ 0.

Conclusions
We developed a MC event generator for nucleon configurations in nuclei which correctly reproduces single nucleon densities and central NN correlations in nuclei. The generator can be used in modeling a wide range of processes. In particular, it would be interesting to explore how it would modify effects of fluctuations of the number of wounded nucleons in the heavy ion collisions.