The Angantyr model for Heavy-Ion Collisions in PYTHIA8

We present a new model for building up complete exclusive hadronic final states in high energy nucleus collisions. It is a direct extrapolation of high energy pp collisions (as described by PYTHIA), and thus bridges a large part of the existing gap between heavy ion and high energy physics phenomenology. The model is inspired by the old Fritiof model and the notion of wounded nucleons. Two essential features are the treatment of multi-parton interactions and diffractive excitation in each NN sub-collision. Diffractive excitation is related to fluctuations in the nucleon partonic sub-structure, and fluctuations in both projectile and target are here included for the first time. The model is able to give a good description of general final-state properties such as multiplicity and transverse momentum distributions, both in pA and AA collisions. The model can therefore serve as a baseline for understanding the non-collective background to observables sensitive to collective behaviour. As PYTHIA does not include a mechanism to reproduce the collective effects seen in pp collisions, such effects are also not reproduced by the present version of Angantyr. Effects of high string density, shown to be able to reproduce e.g. higher strangeness ratios and the ridge in pp, will be added in future studies


Introduction
At hadron collider experiments at RHIC and LHC, protons as well as large nuclei, are collided, and the results are interpreted to obtain better knowledge about the dynamics of the fundamental interactions at high energies. The strong nuclear force plays a central role, but the studies of proton-proton (pp) collisions and heavy ion collisions respectively, are often carried out in quite different ways.
In the case of pp collisions, so-called "general purpose Monte Carlo event generators", such as SHERPA [1], Herwig 7 [2] and PYTHIA8 [3], have been established as cornerstones in aiding our understanding. These event generators have over the last three decades succeeded in simultaneously simulating the dynamics of strong and electroweak processes from very high momentum transfer scales where perturbation theory is applicable, down to scales around Λ QCD , where one must rely on models inspired by analogies to electrodynamics or results from lattice QCD. This has resulted in a remarkably precise description of the majority of observations in proton-proton collisions, which both further experimental and theoretical developments often rely heavily upon.
In high energy heavy ion collisions, the landscape is quite different. Here efforts are more often directed towards signals for the formation of the Quark-Gluon Plasma (QGP), and studies of its properties. The existence of such a phase is demonstrated in lattice calculations and it is presumed to have existed in the hot, early Universe. In this area event generators also exist, but are usually more "special purpose" than "general purpose", each attempting to describe a specific array of observations ascribed to the formation of a QGP. Event generators generating full exclusive events also exist, and the ones most frequently used in analyses investigating particle production mechanisms are, arguably, EPOS-LHC [4], AMPT [5] and HIJING [6]. At least for the bulk event properties, these three generators have for many years defined the "golden standard" for Monte Carlo comparisons to experimental data. In section 7 we outline some of the main similarities and differences between these models and our own.
Several features, which in heavy ion physics are interpreted as a QGP effect, are also observed in pp collisions at the LHC, which may indicate that the dynamics at play in these two types of collision systems are in fact very similar. Two typical examples are enhanced strangeness [7] and the formation of a "ridge" [8]. This immediately raises a challenge for the general purpose pp event generators and their underlying models. If a QGP is indeed formed even in pp collisions, then the effects of such a formation should be included. On the other hand, if the flow-like effects in pp collisions have a different, non-thermal, origin, then it might be possible to capture the general features of nuclear collisions by adding a nuclear structure "on top" of existing pp models. In the present paper we will primarily address the second of these possibilities, presenting a model, henceforth called "Angantyr", which is an extrapolation of pp dynamics to collisions with nuclei, assuming no formation of a thermalised QGP state.
A generalisation of pp collisions to an event generator for pA and AA collisions will have two separate components. First it is necessary to estimate the number of interacting nucleons and binary N N collisions. This is generally performed using the Glauber formalism [9,10]. This formalism is based on the eikonal approximation in impact parameter space, where the projectile nucleon(s) are assumed to travel along straight lines and undergo multiple sub-collisions with nucleons in the target. The importance of including diffractive excitation was early pointed out by Gribov [11], but has often been neglected also in recent applications (see e.g. the review by Miller et al. [10]). As an example, in many analyses the N N interaction has been approximated by a "black disk model", where diffractive excitation of individual nucleons is completely neglected.
Secondly one has to estimate the contribution to the final state from each interacting nucleon. The Angantyr model is here inspired by the old Fritiof model for pA and AA collisions [12,13], with essential features from the idea of "wounded" (or "participating") nucleons [14]. 1 In the early Fritiof model [12] it was assumed that an interacting nucleon suffers a longitudinal momentum exchange with a distribution ∼ dQ/Q, leading to an excited mass ∼ dM 2 /M 2 . When hadronising like a colour string this gives on average a triangular distribution in rapidity. This behaviour was also obtained by Bia las and Czyż in analyses of dAu collisions at RHIC [15].
The Angantyr model is a generalisation of the model for pA collisions presented in ref. [16]. This model was able to reproduce general features in pA collisions, like multiplicity as a function of (measured) centrality, rapidity distributions, and to a certain degree also p ⊥ distributions. It was also here discussed how these observed features rely on the often used model dependent quantities, such as distributions of wounded nucleons and binary N N collisions, and how they are influenced by cross section fluctuations. In this paper we will describe its generalisation, Angantyr, to AA collisions and detail its implementation into the PYTHIA8 event generator.
One essential feature which distinguishes the model for pA collisions in ref. [16], and the Angantyr model in this paper, from the early Fritiof model, is the inclusion of multiple (semi-)hard partonic sub-collisions. At higher energies these play a very essential role. Another difference is the explicit treatment of diffraction. We note here that, if the mass distribution for diffractive excitation can be approximated by dP ∝ dM 2 /M 2 , then the contribution from a diffractively excited nucleon is very similar to the contribution from an average wounded nucleon in the Fritiof model or from the analysis in ref. [15]. The wounded nucleons in Fritiof can therefore effectively represent both non-diffractively and diffractively wounded nucleons.
High mass diffractive excitation was described within the Regge formalism by the triple-Pomeron diagram by Mueller [17]. At higher energies unitarity constraints imply that multi-Pomeron diagrams become important. This leads to complicated schemes for summing up the different contributions [18][19][20]. Diffractive excitation has also been described in the Good-Walker formalism [21], as a result of fluctuations in the internal structure of the projectile. This formalism was first used in high energy collisions by Miettinen and Pumplin in 1978 [22], and later in refs. [23,24]. In ref. [25] it was demonstrated that both the multi-Pomeron and the Good-Walker formalism can in fact be interpreted as different sides of an underlying dynamics based on QCD and BFKL evolution.
In most applications the real part of the N N amplitude is neglected. Diffraction (elastic scattering and diffractive excitation) is then the shadow of absorption into inelastic (nondiffractive) channels. Absorption is here specified by colour exchange between projectile and target, while diffraction corresponds to colour neutral (Pomeron) exchange. In the Good-Walker formalism diffractive excitation is then part of the diffractive beam, when the projectile mass eigenstate is a (coherent) linear combination of states with different absorption probability. An interacting nucleon can then be classified as either "absorptive" ("inelastic non-diffractive", colour connected), "diffractively excited" or "elastically scattered".
Fluctuations in the projectile proton in pA collisions was studied by Heiselberg et al. [26], for estimates of the number of individual N N sub-collisions. This formalism was further developed in several papers (see refs. [27][28][29][30] and further references in there). It is often referred to as the "Glauber-Gribov" colour fluctuation model (GGCF or just GG), and is used in several experimental analyses, e.g. in refs. [31,32].
As discussed in ref. [16], taking averages over target nucleon states is enough for calculations of cross sections and the number of wounded nucleons in pA collisions, provided diffractively excited nucleons are also counted as wounded nucleons. For a generalisation to AA collisions it is, however, necessary to also take into account individual fluctuations in both projectile and target nucleons.
In ref. [16] we introduced the concept of primary and secondary absorptive interactions, when a projectile nucleon is interacting absorptively with more than one target nucleon. The corresponding N N parton-level event could be generated using the full multi-parton interaction (MPI) machinery in PYTHIA8, for both absorptive and diffractive interactions. To generate fully exclusive final states in AA collisions, we have to calculate all subcollisions between a nucleon µ in the projectile and nucleon ν in the target, study the number of multiple sub-interactions for all nucleons µ and ν, and here separate diffractive from non-diffractive (absorptive) interactions. This process is fully described in section 3.
We now return to the question of QGP formation. In the current version of Angantyr the generated partonic states are hadronised using the string fragmentation model in PYTHIA8, without including any final-state collective effects. In this way the model can be used as a starting point for implementing and analysing new models for collectivity. As an example we showed [33] that an enhanced strangeness production can be expected in (high multiplicity) pp collisions, due to overlapping colour strings forming "ropes", in agreement with experimental observations [7]. Furthermore we demonstrated in refs. [34,35] that the enhanced density also ought to give an outward pressure, which may explain the observed flow-like effects in pp scattering.
In the present version of the model we limit ourselves to general features like distributions of particle density in rapidity and p ⊥ , postponing a discussion of flow-like effect to a coming publication. We would like to emphasise, however, that the model can still be used as an important tool for understanding non-flow effects on experimental observables designed to measure flow and other collective behaviours.
In figure 1 we show how the structure described above is put together and tuned in the concrete simulation. Since all parts of the simulation; GG colour fluctuations to Glauber Figure 1: Flowchart showing the programmatic structure of Angantyr. In order to make predictions for heavy ion collisions, several parts of a normal PYTHIA8 simulation needs to be modified, and tuned accordingly. In the flowchart we illustrate how each separate part is tuned to either e + e − , ep or pp data, while no tuning is done to heavy ion data.
generate the number of sub-collisions, the PYTHIA8 MPI model, the parton shower and the hadronisation model rely on a number of parameters, these parameters need to be tuned, and a large part of this paper describes how this procedure is carried out. We want to emphasise from the beginning that all parts are tuned to data from collisions of smaller systems, e + e − , ep and pp, and no tuning is done to heavy ion data. The results can thus be regarded as real predictions depending only on the chosen extrapolation procedure, and not a specific choice of parameters. The layout of the paper follows the workflow of the generation procedure as shown in figure 1, and implemented in PYTHIA8. In section 2 we discuss how to calculate the number of wounded nucleons and the number of individual N N sub-collisions. Here we include fluctuations both in the distribution of nucleons in the nuclei and in the individual nucleon states, both for nucleons in the projectile nucleus and in the target nucleus. We note that a projectile (target) nucleon is fixed in the same diffractive eigenstate through the the passage through the target (projectile) nucleus. If it is then not absorbed, it may end as diffractively excited, when projected to the system of mass eigenstates. Then in section 3 we discuss how to generate the parton-level sub-events for the different kinds of sub-collisions, and in section 4 we describe the procedure for stacking these sub-events together into complete exclusive hadronic final-states in AA. In section 5 we then make a digression to discuss the details of the generation of secondary absorptive sub-collisions, before we present some sample results in section 6. Finally we discuss differences and similarities between our approach and other heavy ion event generators in section 7, before presenting some conclusions and an outlook in section 8.

Nucleon-nucleon sub-collisions in pA and AA
In high energy pp collisions the real part of the amplitude A is small. If this can be neglected, we can define the real quantity If diffractive excitation also can be neglected, the elastic cross section is just the shadow of the absorption, which in impact parameter space is determined by the probability 1 − S 2 .
The elastic, total, and inelastic pp cross sections are then given by dσ NN respectively. The formulations of high energy nucleus collisions in terms of individual nucleonnucleon interactions was carried out by Glauber in a pioneering paper in ref. [9]. In this paper several kinds of fluctuations were neglected. As pointed out by Gribov, and discussed in the introduction, diffractive excitation of individual nucleons is essential, both for cross sections and for final state properties. The Glauber theory is formulated in impact parameter space, where cross sections can be directly interpreted as probabilities. It is then most convenient to include diffractive excitation using the Good-Walker formalism [21], as the result of fluctuations in the nucleon wavefunctions. In this section we shortly discuss the Glauber and Good-Walker formalisms for estimating scattering cross sections and distributions of wounded nucleons and N N sub-collisions. The discussion of effects on the properties of exclusive final states will be presented in section 3.

Glauber formalism
The Glauber formalism is based on the eikonal approximation in transverse coordinate space. Here the projectile nucleon(s) travel along straight lines, and undergo multiple subcollisions with small transverse momenta. Multiple interactions correspond to a convolution of the individual S-matrices in transverse momentum space, which in transverse coordinate space simplifies to a product.
We let b µ and b ν denote the set of positions in impact parameter space for the nucleons in the projectile and target nucleus respectively, and b the separation between the centres of the colliding nuclei. The S-matrix for scattering between nucleus A and nucleus B is the given by (2.2) Here b µν = b µ + b − b ν is the relative separation between the two colliding nucleons N µ and N ν . For pA collisions the product over µ contains only the projectile proton with b µ = 0. As mentioned above, fluctuations were neglected in Glauber's original paper. In the optical limit, with a smooth distribution of nucleons in the nuclei, and where the size of the nuclei are large compared to the range of the N N interaction, the resulting nucleus-nucleus cross sections can be calculated analytically.

Nucleus geometry
The simplest way to include fluctuations in the nucleon positions, b µ , within a nucleus, is to randomly distribute the A nucleons in three-dimensional space according to a Woods-Saxon distribution. More advanced models include correlations in form of a hard repulsive core (e.g. [36,37]), or a more sophisticated description of the two-(or three-) particle correlations between the nucleons within the nucleus [38,39]. Fluctuations in the geometry is taken into account, when new nucleus states are generated for each new event.

Fluctuations in the individual N N interactions, and the Good-Walker formalism
We here shortly describe the Good-Walker formalism for diffractive excitation, assuming for simplicity first that a fluctuating projectile collides with a non-fluctuating target. For a projectile particle with an internal substructure, it is possible that the mass eigenstates differ from the elastic scattering eigenstates. We denote the mass eigenstates Ψ i , with the projectile in the ground state (e.g. a proton) denoted Ψ 0 , while Φ l are the eigenstates to the scattering amplitude T , with T Φ l = t l Φ l . The mass eigenstates are linear combinations of the scattering eigenstates, Ψ i = l a il Φ l . The scattering can now be regarded as a measurement, where the projectile "has to choose" one of the eigenvalues t l , with probability |a 0l | 2 . The elastic amplitude for the ground state projectile is then given by Ψ 0 |T |Ψ 0 = l |a 0l | 2 t l ≡ T , where T is the expectation value for the amplitude T for the projectile. The elastic cross section is then given by We here work in impact parameter space, and the amplitude depends on b. The total diffractive scattering σ diff (including the elastic) is the sum of transitions to all states Φ l : where we have used the fact that Φ l form a complete set of states. Subtracting the elastic cross section we then get the cross section for diffractive excitation, which thus is given by the fluctuations in the scattering amplitude: In a nucleon-nucleon collision both the projectile and the target are fluctuating, leading to single diffractive excitation of the projectile or the target, as well as to double diffraction. The different cross sections are then given by Here · · · p and · · · t are averages over projectile and target states respectively, and subscripts Dt, Dp and DD stand for single diffractive excitation of the target, the projectile, and double diffraction respectively. We note here that while the total cross section depends only on the average of T (b), all other cross sections include also average of T 2 over projectile and/or target states. However, if wounded target nucleons include also diffractively excited nucleons, we see that the corresponding cross section for a wounded target nucleon, σ NN Wt ≡ σ NN abs + σ NN Dt + σ NN DD , can be written (2.7)

Fluctuations in collisions with nuclei
The expression for the amplitude T (b) = (1 − S(b)) in eq. (2.6) can be directly inserted into the amplitude for collisions with nuclei in eq. (2.2) (as before we neglect the real part of the amplitudes). The scattering probability can be regarded as a measurement, after which a projectile nucleon is in one of the eigenstates to the amplitude T , and thus also to the probability for colour connection (the absorption probability) 2T − T 2 . Thus all nucleons are frozen in the same state during the scattering process. (We here neglect the modification when one or a few partons have changed colour in the first encounter.) As a consequence the average of the AA amplitude in eq. (2.2) will include also higher powers of T . However, for pA collisions the multiple sub-collisions imply that the total and wounded nucleon cross sections contain higher moments with respect to projectile fluctuations, but still only the average over the uncorrelated target nucleon states. We also note that these moments should be taken for fixed impact parameters. Thus, to calculate the ratios of e.g. the integrated elastic and total cross sections, it is also necessary to know the b-distribution of the amplitude.
To visualise the effects of fluctuations and diffractive excitation we can study a simple example with a proton colliding with two target nucleons, with and without fluctuations. We assume in both cases that the inelastic N N cross section (including diffractive excitation) is dσ N N inel /d 2 b = 3/4. Case 2: With fluctuations. We neglect the fluctuations in the target, and assume that the projectile state is given by . The states Φ 1 and Φ 2 are here diffractive eigenstates with eigenvalues t 1 = 0 and t 2 = 1. From eq. (2.6) we get for collision with one target nucleon dσ abs /d 2 b = 1/2 and dσ D /d 2 b = 1/4. For two target nucleons we get actually the same result. If the projectile is in state Φ 1 it misses both targets, and if in state Φ 2 , it is absorbed already in the first one. Thus the inelastic cross section is only 3/4 (1/2 for absorption and 1/4 for diffractive excitation) compared to 15/16 in the non-fluctuating case. 2

From cross sections to probabilities
The absorptive cross section in impact parameter space shown in eq. (2.6) is the average where T i,k is the scattering amplitude (and S i,k the S-matrix) for a projectile proton in state i colliding with a target in state k. This expression is always ≤ 1, and it can be directly interpreted as the probability for an absorptive interaction between the projectile and the target. (Such an interpretation is not possible in transverse momentum space, where the cross section has the dimension of momentum to the fourth power.) We note, however, that neither the elastic cross section nor diffractive excitation is the average of an expression depending on only i and j. (The elastic cross section can be written i,j,k,l T i,k T j,l .) When the interaction is driven by absorption, elastic scattering and diffractive excitation is the result of interference between waves, which missed the absorbing target. The cross section for this diffractive scattering is also bounded by 1, and together with absorption it gives a total cross section bounded by 2. A consequence of this feature is that to properly generate events including diffractive excitation for AA collisions in an event generator, it is necessary to, for every projectile nucleon, µ, in state i calculate the average of the amplitude T i,k (b µν ) over all states of each target nucleon, ν, for all impact parameters b µν (and similarly all averages over projectile states i for every target state j). This would give a very slow program, and in section 2.5 we show how to obtain a good approximation.
In pA collisions the picture is, however, much simplified. From eq. (2.7) we note that although the wounded nucleon cross section dσ NN Wt /d 2 b contains one piece from absorption and one piece from diffraction, the sum is always bounded by 1. The question whether a target nucleon ν will be a wounded target (with this definition) in a sub-collision with a projectile in state i can only be answered by yes or no. Therefore the answer yes must have the probability given by the cross section in eq. (2.7). This is used e.g. in applications of the Glauber-Gribov model described in section 2.4.2.

Non-fluctuating models
The simplest approximation for the N N amplitude is the "black disk model", where the target acts as a black absorber. This model has been frequently used in experimental analyses (see e.g. the review in ref. [10]). It is then assumed that two colliding nucleons are interacting, if their separation in impact-parameter is smaller than some radius R. The cross sections are here given by σ inel = σ el = σ tot /2 = πR 2 . As there are no fluctuations, the cross section for diffractive excitation is zero. It is then obvious that the model cannot reproduce the experimental results, which satisfy σ NN el ≈ σ NN D ≈ σ NN tot /4. (Here σ NN D denotes the sum of single and double diffractive excitation.) In the literature it is common to set 2πR 2 = σ (exp) tot , which reproduces the experimental total cross section, but neither the elastic nor the inelastic cross section (when the latter includes diffractive excitation). In later studies it has become more common to choose πR 2 = σ (exp) inel , which reproduces the total inelastic cross section, but gives σ tot . For most applications in pA and AA, the elastic cross section is not very important, but we note that it could still be reproduced by introducing a grayness or opacity of the collision, assuming that within a radius R the scattering amplitude is a constant a between 0 and 1. R and a can then always be adjusted to reproduce both the total and the elastic cross sections (and thus also the total inelastic cross section). Diffractive excitation would, however, still be absent.

Models including fluctuations
In a variation of the opacity model above, the projectile is instead fully absorbed with probability a. This obviously includes fluctuations and thus also diffraction. With the value a = 1/2 we get the cross section ratios σ NN el = σ NN D = σ NN tot /4, in reasonable agreement with experiments. As the model describes the combined fluctuations of the projectile and the target, it is here not possible to separate diffractive excitation of the projectile from that of the target or from double diffraction.
In the introduction we mentioned the "Glauber-Gribov" model for pA collisions, developed by Strikman and coworkers [26][27][28][29][30]. It is there assumed that the fluctuations in the projectile can be described by a distribution in the quantity σ The second relation follows from eq. (2.6).) This formalism, has been used in analyses of pPb data from LHC, e.g. by ref. [31], to estimate the number of wounded or interacting nucleons, which in turn has been used to estimate the centrality for the collision. The quantity σ is then normally rescaled so that the integral in eq. (2.8) gives the inelastic rather than the total cross section. We note that, as the fluctuating quantity σ includes the fluctuations over projectile states, but averages over target nucleon states, we see from eq. (2.7) that what is counted as wounded nucleons includes diffractively excited nucleons. As discussed in section 2.3, the cross section in eq. (2.7) also determines the probability distribution for wounded nucleons, but we want to emphasise that the differential cross section T (b) t is needed for all values of the impact parameter b. In ref. [28] this is assumed to be Gaussian ∝ exp(−b 2 /2B(σ)), with a slope parameter B(σ) proportional to σ, in order to satisfy the unitarity constraint T (b) ≤ 1.
In ref. [16] we investigated the fluctuations in the nucleon cross sections using Mueller's dipole approach to BFKL evolution [40,41] as implemented in the DIPSY Monte Carlo program [24,42,43]. The model is formulated in impact parameter space, and includes also a set of sub-leading corrections beyond the leading log BFKL approximation. Non-linear effects are introduced by the "colour swing" mechanism, which suppresses large dipoles, corresponding to k ⊥ below a saturation scale. BFKL evolution is a stochastic process, and the result was here that the fluctuations have a longer tail out to large cross sections compared to the distribution in eq. (2.8). Rather than the Gaussian suppression assumed in [26], we found a distribution more similar to a Log-normal for the b-integrated and target-averaged σ: To also describe the b-dependence of T (b) t , we used a semi-transparent disk approximation with the elastic amplitude The parameters (Ω and σ 0 ) in P tot (σ) and T 0 in eq. (2.10) could here be fitted to σ NN tot , σ NN el and σ NN Wt taken from experimental data, to obtain a Glauber-like calculation for pA. Together with the parton-level stacking also proposed in [16] we then also obtained a fair description of e.g. the observable used by ATLAS in [31] for estimating centrality, as well as the corresponding pseudo-rapidity distributions as a function of that centrality.
We note that the stochastic nature of BFKL evolution has also been studied by Iancu, Mueller and Munier in ref. [44]. When the probability for a dipole splitting is small, the mean field approximation in the Balitsky-Kovchegov equation does not properly describe the probability for rare events with large cross section. In ref. [44] they studied the fluctuations in the saturation scale, Q s , and showed that for asymptotic energies the width of the distribution in ln(Q s ) is growing proportional to ᾱ ln(s), with a tail to large Q s -values in qualitative agreement with eq. (2.9).

Nucleon fluctuations in AA collisions
As mentioned in section 2.2.2, to study pA collisions also higher moments over projectile fluctuations are needed. When we now want to generalise the formalism to AA collisions, both projectile and target nucleons are frozen under the collision (but still uncorrelated). This implies that we must be able to calculate not only T (b) n t p , but any moment T (b) np nt p t . To cope with this situation we need a formalism which can give the amplitude T ik (b) for any combination of projectile state i and target state k.
We noted that the Log-normal distribution in eq. (2.9) is quite similar to a Gammafunction, and for technical reasons and the fact that the sum of two Gamma distributed random variables is also Gamma distributed, we will use that instead to model fluctuations in the radius, r, of a nucleon: We then also use a slightly different elastic amplitude where the opacity of the semi-transparent disk now depends on r p and r t : This introduces one more parameter, σ t , (besides k and r 0 in eq. (2.11)) and by this varying opacity makes it possible to get a reasonable fit to all the cross sections in eq. (2.6), as well as the elastic slope parameter B = −d ln σ NN el /dt| t=0 , for a wide range of energies. The result for √ s NN = 5 TeV is shown in table 1.

Determining the interaction of nucleon sub-collisions
We now want to take all pairs of colliding nucleons in an AA collision, and for each of these select which kinds of interactions are possible. At high energies all nucleons are frozen in their (random) states during the passage through the opposite nucleus. The probability for an absorptive interaction between nucleon µ (in a state i with radius r iµ ) in the projectile and nucleon ν (in state k with radius r kν ) in the target, is then directly given by eq. (2.6) as given by eq. (2.12). To estimate the probability for a diffractive excitation of a given nucleon is more difficult, as diffractive excitation is part of the shadow scattering caused by absorption, to which all encountered nucleons contribute.
We showed in ref. [16] that for a given state of a projectile nucleon the probability that a given target nucleon is absorptively or diffractively wounded in the interaction is given by the average over the possible states of the target (c.f. eq. (2.7)) and that this probability factorises for all nucleons in the target nucleus. However, in AA the symmetry between projectile and target complicates things further, as we need both a specific state and the average over all states for all nucleons.
In Angantyr this is handled by generating two states (one primary, r, and one auxiliary, r ′ ) for each nucleon in the nuclei. The primary one is used to calculate the probability of an absorptive N N interaction, while the secondary is used to statistically sample the average state of each nucleon. The algorithm ensures that on average (over the four possible combinations of states in an N N interaction) we get the correct probability of the projectile and target nucleon being absorptively and diffractively wounded.
The technical details of this algorithm is presented in appendix A, while here we will only show that it works as expected. In table 1 we give an example where we have fitted the parameterisation of the fluctuations according to eqs. (2.11) -(2.13) to the default parameterisation of the semi-inclusive cross sections in PYTHIA8. This default parametrisation [45,46] does not necessarily agree well with cross section measurements from LHC [47], and it is possible for a user to easily supply their own cross sections as input to the fit. The last line denoted "generated" shows the results from generating N N collisions in Angantyr for √ s = 5 TeV. We see that the absorptive cross section comes out close to the input one, and also the wounded cross sections, σ Wp and σ Wt are reasonably well reproduced. However, we see that the individual diffractive excitation cross sections are not reproduced, nor is the elastic ones. However, for the final states in AA collisions, we are mainly interested in getting the absorptive and wounded cross section right, so even if our procedure probably can be improved, we are quite satisfied with this result. 2.2 -- Table 1: Fitting the values of input cross sections for pp collisions at √ s = 5 TeV and using the resulting fluctuations in a generation and different collision types. B is the elastic slope −d log σ el /dt| t=0 . The cross sections used as "input" were taken from the default parameterisation in PYTHIA8. The line "model" shows the results of a fit to the model in eqs. (2.11) -(2.13). The line "generated" finally shows the result of the approximation discussed in this subsection and in the appendix. The fitting procedure assumed a 5-10% uncertainty on the input values, and the statistical uncertainty on the presented output values are around and below 0.5%.

From wounded nucleons to exclusive final states
In the wounded nucleon model, as formulated by Bia las and Czyz [15], each wounded nucleon contributes to the final state multiplicity distribution, according to a single nucleus emission function F (η), giving a total multiplicity of: Here w p|t denotes the number of wounded nucleons from left and right respectively, calculated for a given centrality class, defined by impact parameter. In the wounded nucleon model, F (η) must be extracted from data, and depends on centrality class [48], but a crucial feature of the model is that eq. (3.1) reduces to the pp multiplicity distribution for w p = w t = 1. The Angantyr prescription for generating exclusive final states has conceptual similarities with the wounded nucleon model. But instead of extracting an emission function from data, MPI events from PYTHIA8 are used. We will in this section briefly review the PYTHIA8 MPI model, and motivate the addition of additional MPIs from multiple wounded nucleons to the model.

Multiparton interactions in pp collisions
In the PYTHIA8 MPI model [49], all partonic sub-collisions are to a first approximation treated as separate QCD 2 → 2 scatterings 3 . Since the cross section diverges at low p ⊥ , it is regularised using a parameter p ⊥0 which depends on the collision energy, giving: This cross section is then folded with parton densities to get a relative probability for each additional sub-scattering. The densities are rescaled according to an overlap Figure 2: Schematic pictures of multi-parton interactions in a pp collision. The y-axis should be interpreted as rapidity. All initial-and final-state radiation has been removed to avoid cluttering. Each gluon should be interpreted as having two colour lines associated with it, which in the subsequent string hadronisation will contribute to the soft multiplicity. In (a) the colour lines for both sub scatterings stretches all the way out to the proton remnants, while in (b) and (c) the secondary scattering is colour-connected to the primary one.
function using some assumption about the matter distribution in the colliding protons and an assumed impact parameter.
In figure 2a there is an illustration of an event with two sub-scatterings (in red and black) which we have assumed are both of the type gg → gg. Note that in the PYTHIA MPI model all incoming and outgoing partons would be dressed up with initial-and final-state radiation, but these have been left out of the figure to avoid cluttering. With completely uncorrelated sub scattering, one would assume the colours of the incoming gluons would also be uncorrelated, and since each gluon carries both colour and anti-colour one would naively think that in the subsequent hadronisation phase, there would be four strings stretched between the proton remnants and giving rise to particle production over the whole available rapidity range. Again to avoid cluttering of the figures, we ask the reader to simply imagine two colour lines (strings) stretched along each gluon and that the vertical axis can be loosely interpreted as rapidity.
Already in the original paper [49] it was realised that it was basically impossible to reproduce data if each sub-scattering was allowed to add particles in the whole available rapidity range. Especially sensitive to this was the multiplicity dependence of the average particle transverse momenta, and to rectify this the MPI model in PYTHIA was modified so that additional sub-scatterings almost always was colour connected to outgoing partons in previous sub-scatterings. This is illustrated in figure 2b and c, where the colour correlation (a) (b) Figure 3: A schematic picture (c.f. figure 2) of multiple scattering between one projectile and two target nucleons (e.g. in a pd collisions). In (a) the second interaction is directly colour connected to the first one, while in (b) the second nucleon is only diffractively excited by a Pomeron exchange. Both cases give rise to final string configurations that will contribute in the same way to the final state hadron distribution. between the two sub-scatterings gives rise to a colour flow as if they were (perturbatively) connected. In this way the multiple scatterings can give rise to increased average transverse momentum from the partons coming from extra sub-scattering, without increasing the multiplicity of soft particles due to the strings stretched all the way out to the proton remnants.

Multi-parton interactions in a pA collision
We now turn to the case of a pA collisions and imagine the target proton interacting absorptively with two nucleons in the nuclei. To be true to the PYTHIA MPI model we should simply redefine the overlap function using the matter distribution of the two target nucleons. In principle this can surely be done, however, technically we found it almost forbiddingly difficult.
Instead we note that the handling of colour correlations in the pp model would typically result in string topologies corresponding to the sketch in figure 3a. The primary scattering looks like normal scattering between the projectile and one of the target nucleons, while the secondary scattering is now between the projectile and the other target nucleon. Since both target nucleons have been found to be absorptively wounded, the secondary scattering must be colour connected to the second target nucleon, while in the direction of the projectile it looks like a normal secondary scattering.
We also note that we would get the same colour topology, and hence the same distribution of particles, if the second sub-scattering was a separate single (high-mass) diffractive excitation event, which in PYTHIA8 is handled as a Pomeron-proton collision. This is illustrated in figure 3b, where the Pomeron is shown as a green zigzag line. This is the picture where secondary absorptive wounded nucleons contributes to the final state particles as if they are produced in a single diffractive excitation event is what we use in the Angantyr model.
The procedure will therefore be to decide which of the the two absorptive interactions is to be considered the primary one, and treat this as a completely normal non-diffractive multiple scattering event in PYTHIA. The secondary scattering will be generated as a single diffractive excitation event in PYTHIA. Also here there may be additional multiple parton scatterings, but they will be treated as multiple scatterings in the Pomeron-proton system, which is standard in the high-mass diffraction machinery in PYTHIA.
Referring back to eq. (3.1), this means that we are modelling the single nucleus emission function F (η) using high-mass diffractive excitation events. We do not expect them to necessarily look like ordinary diffractive event, but we nevertheless use the diffractive machinery in PYTHIA8. In section 5 we will describe how we modify this machinery in order to try to fulfil the requirement that F (η) + F (−η) (i.e. w p = w t = 1 in eq. (3.1)) would reproduce the distribution in a normal non-diffractive pp event in PYTHIA8.
The two different sub-events are then merged together so that the elastically scattered proton in the diffractive event is discarded, and the momentum of the Pomeron is instead taken from remnants of the projectile proton.
The assumption in [16] was that the momentum fraction of the Pomeron in such diffractive events can be taken to be distributed approximately as dx IP /x IP , which means that the mass of the diffractive system is given by dM 2 This is approximately what one has found for normal high-mass diffractive events and it is the same assumption as in the old Fritiof model. We do not have a solid explanation why this should be the case. In [16] we gave some handwaving arguments based on AGK cutting rules and the similarity between triple-Pomeron diagrams in diffractive N N scatterings and (doubly) non-diffractive proton-deuterium scattering, but in the end the best argument for this choice is that it seems to work very well.

Multi-parton interactions in an AA collision
Going one step further in complexity we now consider AA collisions. In figure 4 we illustrate the situation when two nucleons in one nucleus collides with two nucleons in the other, and all four possible N N interactions are absorptive. We find that there is three ways of doing this which are consistent with our pA model. Either, as in figure 4a, we can model it as two primary absorptive interactions, or as one primary and two secondary interactions, where the second of these can either be coupled to the primary interaction (b) or to the first secondary one (c).
All three cases will give us four absorptively wounded nucleons, and in Fritiof and the original wounded nucleon model there would be no distinction between the cases. In the Angantyr model we do, however, want to differentiate between these, and in the following section we will describe a procedure to classify all N N interactions in a AA collisions.

Generating and combining parton-level N N events
In general, each nucleon in the projectile nucleus may interact with several nucleons in the target nucleus an vice versa. When building up the final state by stacking parton level nucleon-nucleon events we need to concentrate on the most important ones first. This is in line with the general philosophy in PYTHIA8, that harder processes always are considered before softer ones. After having gone through all pairs of projectile-target nucleons and determined their interactions as outlined in the previous section, we therefore order all these interactions in increasing nucleon-nucleon impact parameter, b µν . We will then go through this list several times, treating one kind of interaction at the time, starting with the absorptive interactions, as they will give the back bones around which we will build up the full event. As soon as an N N interaction has been selected a corresponding sub-event will be generated with the standard PYTHIA8 minimum bias model, and the corresponding nucleons are marked as already interacted. If an N N interaction is found in the list where one of the nucleons has already interacted, this will be labelled secondary and the generated sub-event will be added to the sub-event to which the already interacted nucleon belongs, as described in [16] and detailed below.

Selecting primary absorptive collisions
The first pass over the potential N N interactions, we will only look at absorptive inter-actions. This will give us a set of N ′ abs primary absorptive collisions and a set of N ′′ abs secondary ones (where one of the nucleons already has already been absorptively wounded in another interaction).
For each of the primary ones we now generate an inelastic non-diffractive minimum bias event in PYTHIA8, each of which will give a separate sub-events. However, since the procedure takes sub-collisions with small b µν first, the primary absorptive events should typically be a bit harder and have higher multiplicity than the secondary ones. In [16] this was handled by telling PYTHIA8 to generate N ′ abs + N ′′ abs events, but only keeping the N ′ abs ones with smallest impact parameter (as reported by PYTHIA8). For the method described here we have instead implemented directly in PYTHIA8 a way to specify by hand which impact parameter you want a given minimum bias event to have, which makes thing a bit more efficient, and also gives a noticeable improvement on the description of some observables, as discussed below in section 6.
Just as in standard PYTHIA8 it is easy to specify signal processes rather than only consider minimum bias events. This may be used to simulate triggers on hard jets more efficiently, or to e.g. produce Z-tagged jets in central AA collisions [50] or top events in pAcollisions [51] or AA. The way this is done is simply to substitute the hardest absorptive primary event with a corresponding signal event, and reweighting the event with a factor to get the correct cross section. For signal processes with a large cross section the possibility to have additional signal processes in the same event is also taken into account, however for technical reasons at most N ′ abs signal sub-events can be included in each event 4 . In the current implementation we assume that minimum bias processes are basically iso-spin invariant, and all such sub-events are generated as pp events in PYTHIA8, flipping by hand the iso-spin of a remnant quark or di-quark afterwards in case the corresponding nucleon was actually a neutron, to conserve total charge. Signal processes are, however, not necessarily isospin invariant. To account for this, we generate pp, pn, np, and nn collisions separately for all signal processes. To decide what type of collision should be generated, all nucleons in the colliding nuclei are marked as either protons or neutrons, under the assumption that neutrons and protons are distributed evenly in the nucleus.
One should note that measurements of proton and neutron distributions in e.g. lead at low energies [52] have indicated that the neutron distribution reaches further out than the proton distribution, giving rise to a "neutron skin" effect. It has been pointed out [53] that this could give rise to effects at the 10% level in selected observables in peripheral PbPb collisions. It has also been pointed out [54] that one could in principle use this effect to design different centrality measures, especially in the case of asymmetrical collision systems. Currently we know of only one very recent Glauber calculation including such effects [55], and in the present version we have left them out entirely 5 .

Adding secondary absorptive interactions.
Once the back-bone sub-events have been generated we go through the list again, this time only looking at the secondary absorptive interactions, in which one of the participating nucleons has already been included in a generated primary absorptive sub-event. As described in [16] we will generate these secondary absorptive sub-collisions as if they were single diffractive excitation events. We here use the the standard PYTHIA8 diffraction machinery, but with important modifications detailed in section 5 below.
The final state generated for a given secondary absorptive interaction is then added to a primary absorptive sub-event. The elastically scattered proton is removed and the energy and momentum it had given to the excited nucleon is instead taken from the remnants of the nucleon in the primary sub-events.
It may very well happen that there is not enough energy left in the remnants in the primary sub-event to allow for the addition of a diffractively excited state. In that case it is possible to try again and maybe generate a diffractive event with lower M X . There is a parameter in the program the limits the number of tries allowed, and if the maximum is passed, the corresponding secondary absorptive interaction is simply discarded (although the corresponding nucleon still has the chance to become wounded in another secondary interaction).
The way secondary nucleon interactions are selected according to the N N cross sections, does not take into account possible effects of energy-momentum conservation, therefore it makes sense to try to take such effects into account in this a-postiori way. The parameter we introduced should not be taken as the final word in the matter, but at least it allows us to investigate the effects of energy-momentum conservation.

Adding other N N interactions
Having taken care of all absorptive interactions we continue with the other interaction types in much the same way: 6 1. Generate all primary double-diffraction sub-events.
2. Take all secondary double diffraction interactions and generate these as single diffraction and add them to the corresponding primary sub-event.
3. Take any primary single diffractive interactions that are left and add these as corresponding primary sub-events.
4. Add any secondary single diffraction interactions added to previously generated primary sub-events.
In this way we have generated a set of parton-level sub-events, which we now can join together in a single AA event together with the non-interacting projectile and target Charged particle η at 7 TeV, track p ⊥ > 500 MeV, for N ch ≥ 1 η 1/Nev dN ch /dη Figure 5: Illustration of the shape of the multiplicity function in eq. (3.1) using the η distribution as measured by ATLAS in [56]. The black and red lines are the shapes of standard non-diffractive and singe diffractive events from PYTHIA8 respectively. The green dashed and blue dot-dashed lines are single diffractive events generated by PYTHIA8 using the modifications presented in [16] and the modifications presented in this article respectively. For each single diffractive line there is also a pale line corresponding to adding the mirror image to emulate a non-diffractive distributionà la Fritiof.
nucleons, bunched together in two remnant nuclei, 7 and the parton-level AA event is finally handed back to PYTHIA8 for hadronisation and decay of unstable hadrons.

Modifications of single diffractive to secondary absorptive
In section 3 and in ref. [16] we argued that secondary absorptive interactions will contribute to particle production in the same way as a single diffractive (SD) excitation event (c.f. figure 3). Assuming that such SD events produce a simple flat string with mass distributed as dM 2 X /M 2 X , this would naively give a triangular shape of the F (η) wounded nucleon emission function in eq. (3.1).
We will use the SD excitation machinery in PYTHIA8, where at high energies the diffractive systems are much more complicated than a single string. As described in more detail below, it models the diffractive excitation as a non-diffractive (ND) interaction between the target nucleon and a Pomeron emitted from the projectile (in the spirit of Ingelman and Schlein [57]), and this is then treated with the full MPI machinery as if the Pomeron was a hadronic object with parton densities. In figure 5 we show the average multiplicity as a function of pseudo-rapidity for ND events, and compare it to SD events from PYTHIA8 using the default settings. Clearly we get a somewhat triangular shape for the SD events (SD(def) in the figure), and adding the multiplicity from target and projectile excitation, we get a shape similar to the ND shape, fully in accordance with eq. (3.1) in the case of w p = w t = 1.
In ref. [16] we noticed that using the default PYTHIA8 SD machinery for secondary absorptive collisions resulted in too low activity in pA and tried different modifications to increase the multiplicity. One of these modifications included increasing the gluon density in the Pomeron, which is also shown in figure 5 (SD(glu)).
Here we will try to be more systematic in our approach to modify the default PYTHIA8 SD machinery. Looking at figure 3b, it is clear that the rapidity region close to the one close to the direction of the two target nucleons will be our main focus. Here we note that we could equally well have chosen the second nucleon to be in the primary interaction and the first nucleon to be in the secondary, and would then want to have the same distribution of particles. This means that we want the single diffractive event to look as much as possible as a non-diffractive event close to the direction of the two nucleons. We have therefore investigated several different modifications of the SD model and for different diffractive masses we have studied particle distributions in different pseudo-rapidity intervals and compared these with the corresponding particle distributions in the same intervals for ND events.
In the end we settled for a new modification (labelled SD(new) in figure 5), which is the default way of generating secondary absorptive interactions as of version 8.235 of PYTHIA8 8 . To motivate this, we first need to take a closer look at the SD machinery in PYTHIA.

High-mass diffractive excitation and secondary absorptive
There are more than one way of generating diffractive events in PYTHIA8, but here we will only concern ourselves with the soft diffraction used for minimum bias events. Also here there are two treatments depending on the mass, M X . For low masses, 10 GeV, the excited system is modelled as a simple longitudinally stretched string. In an AA collision, such small excitations will typically be mixed up with the nucleus remnants in the very forward and backward regions and we will here mainly concentrate on high-mass diffraction, which contributes also in the central rapidity region as seen in figure 5.
For high-mass diffraction, PYTHIA treats a proton-Pomeron collision as a normal nondiffractive (ND) hadron-hadron collision and uses the whole MPI machinery with initialand final-state parton showers. This means that there will be multiple 2 → 2 semi-hard partonic scatterings given by Here x IP denotes the fraction of the target proton momentum taken by the Pomeron; β is the fraction of the Pomeron momentum taken by the parton j; and x 1 is the fraction of the projectile proton momentum taken by parton i. Furthermore we have the parton densities in the proton, f p i , and the corresponding densities in the Pomeron, f IP j . Finally we have the flux factor F (x IP ) controlling the diffractive mass given by M 2 X = x IP s. In the following we will assume a flat distribution in log (M 2 X ), in which case F (x IP ) is just a constant.
The partonic cross section dσ ij (p 2 ⊥ ) diverges for small p 2 ⊥ , and although it is regularised as in eq. (3.2) the integrated partonic cross section may still exceed the total non-diffractive IP p cross section for a given M X . In the PYTHIA MPI model this is then interpreted as the possibility of having several sub-scatterings in each collision, with the average number of sub-scatterings given by Here the default value of the of the non-diffractive IP p cross section, σ IP p ND (M X ), is just set to a constant 10 mb. The p 2 ⊥ integral is over the full available phase space, all the way down to zero, but with theσ ij regulated as in eq. (3.2). The parameter p ⊥0 here varies as a small power of M 2 X , in the same way as the p ⊥0 in normal pp scatterings varies with s. In figure 6 we show the resulting pseudo-rapidity distribution of charged particles for different values of M X for diffractive events from PYTHIA8 with √ s = 5 T eV . 9 Here we see the expected behaviour with a large rapidity gap for smaller M X , typical for diffraction. When we want to use the diffractive excitation in PYTHIA to model the secondary absorptive interactions, we want to make the event in the target proton direction to look as much as a normal non-diffractive pp event as possible, and in particular we want the whole event to look approximately the same in the limit M 2 X → s. From the figure we see that this is not quite the case for the default diffraction parameters in PYTHIA8. We also see that the modifications we presented in [16] seems to be a bit too forceful.
Looking at eqs. (5.1) and (5.2) it is easy to see that we can increase the multiplicity by either increasing the general activity by modifying the Pomeron parton densities (as is done in SD(glu) in figures 5 and 6), or we can try to increase the number of sub-scatterings by e.g. adjusting the free parameter σ IP p ND (M X ). We will here look at both these options by studying eq. (5.2) more closely. Studying the average number of sub-scatterings for a fixed rapidity, y = log(x 1 /βx IP )/2, we get If we now compare this to the same for standard non-diffractive pp events, we see immediately that if we modify the Pomeron parton density and make it x IPdependent, βf IP j (β, p 2 ⊥ ) → x IP βf p j (x IP β, p 2 ⊥ ), and at the same time make the total nondiffractive pIP cross section as well as the soft regulator, p ⊥0 , independent of M X , i.e., σ pIP ND (M 2 X ) → σ pp ND (s) and p ⊥0 (M 2 X ) → p ⊥0 (s), we will get very similar expressions. They will not be exactly the same, since the kinematical limits p ⊥ will differ, especially for small M X . Also, for technical reasons, PYTHIA8 will adjust the selected p ⊥0 for each M X value to ensure that the average number of scatterings is always larger than one, effectively making low M X events softer.
The resulting modification is shown in figure 6 as the lines labelled SD(new), and we see that the multiplicity in the proton direction is not much improved at small M X , but at large M X it traces the non-diffractive quite well.
In the next section we will look in more detail on the particle distributions in the rapidity regions where we want the secondary absorptive sub-events to resemble normal non-diffractive events in PYTHIA.

Comparing primary and secondary absorptive sub-events
From figure 6 we see that SD final state particles only populate the rapidity region corresponding to the colour exchange between the Pomeron and the proton (c.f. figure 3b). We will here investigate further to what extent the SD events generated by PYTHIA (with or without modifications) look the same as the ND events in this region. To do this we will study the distribution of particles in different pseudo-rapidity slices for different values of the diffractive mass, M X . In these slices we have looked at standard minimum bias observables based on charged particles, such as average multiplicity (shown in figure 6), the distribution in multiplicity (N ch ), the transverse momentum distribution (p ⊥ ), the distribution in summed ( p ⊥ ) and average ( p ⊥ ) transverse momentum for particles within one unit of η, and average transverse momentum as a function of multiplicity ( p ⊥ (N ch ) ). Naturally, we do not expect these observables to look the same for a diffractively excited system and a full non-diffractive event. Close to the rapidity gap, we are in the fragmentation region of the Pomeron remnant, and here the transverse momentum of final state particles are severely restricted by the kinematics. Also close to the proton fragmentation region, the transverse momenta are limited by kinematics, but here we expect the SD and ND events to look very similar, and indeed we find that they do.
Here we will concentrate on the rapidity regions around the plateau of each M X , and in  figure 6) for mass bins with M X ≈ 70, 500 and 900 GeV respectively. As for the overall multiplicity we find that the default SD machinery, (SD(def)), is quite far from the ND observables in the same rapidity slice. The SD(glu) modification is much closer, but overshoots quite significantly at large M X in the multiplicity distribution (figure 7) and p ⊥ (figure 8). The SD(new) curve typically lies between the two, except for the average transverse momentum the new modification gives a distribution fairly close to the non-diffractive result. In particular the dependence of the average transverse momentum on the multiplicity is known to be very sensitive to the handling of the multi-parton interactions [49], and here we see that SD(new) is quite close to the ND curves here, as may be expected from comparing eqs.  The fact that the p ⊥ distribution in SD(glu) in figures 8 is much harder than in standard ND events would be a problem for the description of the centrality observables used in pA and AA, which are often based on the total transverse activity in the forward/backward region (see section 6.2 below).
In section 4.1 we explained how the impact parameter obtained for each N N subcollision is used as input to PYTHIA8. Here small impact parameters will lead to more multiple scatterings for primary absorptive sub-events. The same impact parameter dependence is also used for secondary absorptive sub-events. It is therefore interesting to compare the SD events with ND events for a specific impact parameter. In figure 9 we show typical examples of such comparisons for impact parameters slightly smaller and larger than average. Comparing with the corresponding distributions in figures 7 and 8, we Sum p ⊥ , M X ≈ 500 GeV and -3< η <-2 Figure 9: Multiplicity, p ⊥ and p ⊥ of charged particles for different modifications of SD events with M X ≈ 500 GeV compared to ND events in the pseudo-rapidity interval −3 < η < −2. All events were generated at fixed impact parameter, b/ b = 0.9 (top panel) and 1.3 (bottom panel. The lines are as in figure 7. see that the difference between the SD and ND curves tend to diminish with increasing impact parameter, which is good, since by construction the secondary absorptive interactions are at larger impact parameter than the primary ones. In figure 9 we did not show curves for SD(glu), and in the following we will disregard this option completely. The modifications there are too severe and somewhat ad-hoc, resulting in far too large effects especially on particle production at high M X (figure 6). We will also disregard the SD(def) option, as it produces too few particles in the nucleus' fragmentation region in pA collisions [16]. SD(new) does not give a perfect reproduction of the ND distributions, and we do not expect any SD model to do that, due to phase space constraints.
The conclusion from the analyses in this section is that SD(new) provides an overall fair description, as well as being more theoretically appealing than the other variations. The SD(new) is therefore, since version 8.235 of PYTHIA, the default model for secondary absorptive sub-events in Angantyr.

Sample results
All results presented here are generated with PYTHIA8 version 8.235 using default settings 10 . 10 Since Angantyr is the default heavy-ion model in PYTHIA, it suffices to specify suitable nuclei as beam This means in particular that: • the nucleon distributions in the nuclei are generated according to the formulae in [37] using the hard-core option, where parameters are tuned to low-energy eA data; • the impact parameter is sampled using a Gaussian distribution with a width large enough to have fairly uniform weights; • the fluctuations in the nucleons were modelled according eqs. (2.11) -(2.13), fitted the default parameterisation of semi-inclusive cross sections in PYTHIA8; • the different N N interactions were classified using the procedure described in sections 2.5 and 4; • the sub-events were generated with the default PYTHIA8 minimum-bias machinery, except for the secondary absorptive ones, where the modifications in section 5 was used.
As with most things in PYTHIA8, there are many options beyond the default behaviour in Angantyr, and there are also so-called user hooks where the user can implement alternative models for e.g. the nucleon distribution, impact-parameter sampling and modelling of fluctuations. There are also a number of parameters in Angantyr that influences the generation of collisions involving nuclei, but most of these can be fitted to pp data. In fact, there are only two parameters that clearly influences the results presented here, which cannot be tuned to pp data. One is the distribution of diffractive masses used in the generation of secondary absorptive sub-events. Here we have assumed a distribution ∝ dM 2 X /M 2(1+ǫ) X where we have simply chosen ǫ = 0 as in the original wounded nucleon model as implemented in Fritiof. The other was mentioned in section 4.2 and is related to energy-momentum conservation when adding secondary sub-events. The default is to simply veto a secondary N N interaction if there is not enough energy left in the corresponding remnant nucleon in the primary sub-event. An alternative is to instead generate a new secondary sub-event (regenerating M X ) to see if that one can be included. 11

pp results
We begin by using the Angantyr generation for the simplest of nuclei, i.e. for pp collisions. Since we actually use the PYTHIA8 minimum bias machinery, we need to make sure that typical minimum-bias observables are reproduced as well when using Angantyr. We expect some differences since all semi-inclusive cross sections are not exactly reproduced in the generation, as explained in section 2.5. Furthermore the distribution in impact parameter is not the same, and since this directly affects the amount of MPI it is important make sure that make the translation between the two works, at least on average.
In PYTHIA8, the impact parameter is by default chosen according to an exponentially falling overlap function, while in Angantyr it is determined by the fluctuations and opacity particles to reproduce the results presented here. 11 The number of attempts allowed for this is governed by the parameter Angantyr:SDTries.  Figure 10: The default PYTHIA8 description of some typical minimum-bias observables in pp, compared to the description using the Angantyr machinery. The latter is given for a range of values of b scale . For comparison we show data from ATLAS [56] as implemented in Rivet [58]. functions in eqs. (2.11) -(2.13), and it is not straight forward to translate directly between the two. In principle one could try to implement the Angantyr distribution as an option in the PYTHIA8 MPI machinery, which then would require a full retuning to pp data. Here we have decided to instead implement a simple scaling factor, b scale , so that for absorptive (non-diffractive) events, which is set to a value ensuring that Angantyr gives approximately the same results as PYTHIA8 for typical pp minimum-bias observables. In figure 10 we see that our tuned value of b scale = 0.85 fairly well reproduces the PYTHIA8 results and gives approximately the same level of agreement with data. 2) for pPb collisions at √ s NN =5 TeV. Data from ATLAS [31] is compared to results from Angantyr. The blue and green lines correspond to allowing Angantyr to retry adding secondary sub-events that fail due to energy-momentum conservation, allowing 2 and 4 attempts respectively. The table shows the resulting bin edges when dividing up in percentiles for the experimental and generated data respectively.

pA results
Comparing to pA data means that we need to consider the concept of centrality, which is used in almost all published experimental heavy ion results. Centrality is based on a final-state observable that is assumed to be correlated with the overall impact parameter of a collision. Typically, this observable involves the activity (multiplicity, transverse energy) close to the direction of the nuclei, and other observables are then conventionally presented in bins of percentiles of this centrality observable.
We will here use the centrality observable defined by ATLAS in [31], which is based on the summed transverse energy in the pseudo-rapidity interval [−4.9, −3.2]. As seen in figure 11, Angantyr is able to reproduce fairly well the measured distribution. However, it should be noted that the experimental distribution has not been corrected for detector effects, so it is difficult to draw firm conclusions about the performance of the model. In the figure we also show what happens if we allow Angantyr to retry adding secondary sub-events that fail due to energy-momentum conservation, as discussed in section 4.2. We see that it does have an impact on the most central collisions.
When we want to use this centrality measure we now have the option to divide it into percentile bins using the measured distribution or the generated distribution, and since they do not exactly agree we will get somewhat different bins, as is shown in table in figure 11.
In figure 12(a) we show the average charged particle multiplicity as a function of pseudo-rapidity measured in the centrality bins defined in figure 11. It is important to remember that even if this is presented as the centrality dependence of the pseudo-rapidity distribution, what is in fact measured is the correlation between the transverse energy flow in the direction of the nuclei and the central multiplicity. In the figure we therfore show two sets of lines generated with Angantyr with the two different binnings presented in figure 11. Clearly the difference is here not significant, which is an indication that Angantyr fairly well reproduces the centrality measure. And the fact that neither curve is far from the experimental data 12 gives a strong indication that the Angantyr is a reasonable way of extrapolating pp final states to pA.
Comparing to the results we presented in [16], the description of data has been much improved. The main reason for this is the more careful treatment of secondary absorptive sub-events, but the new handling of the impact-parameter dependence in the primary absorptive events has also somewhat improved the description of data. Within our model it is possible to look at the actual centrality of an event in terms of the generated impact parameter, and in figure 12(b) we show a comparison between the pseudo-rapidity distribution when binned in percentiles of the generated impact parameter and when binned in the generated E P b ⊥ distribution. Clearly, in the Angantyr model, the binning in E P b ⊥ is not very strongly correlated with the actual centrality in impact parameter. This is especially the case for the most central collisions. The reason for this is the fluctuations modelled in Angantyr, both in the number of wounded nucleons and in the correlation between the number of wounded nucleons and the activity in the direction 12 The η-distributions in figure 12(a) has been corrected for detector effects. of the nucleus.
To study the fluctuations further we show in figure 13 the average number of wounded partons as a function of E P b ⊥ -centrality, both for Angantyr and for three Glauber-model fits performed by ATLAS in [31]: one using standard calculation without fluctuations, and two using the fluctuating cross sections in eq. (2.8) with different Ω-parameters. Clearly we see that Angantyr has larger fluctuations than these standard calculations. In figure 13 we also show the number of wounded nucleons in percentile bins of generated impact parameter. As expected the dependence is very weak for the most central bins (0 − 30%), confirming here that the ATLAS centrality measure mainly picks up the fluctuations in the number of wounded nucleons in this region, and does not correlate very well with the actual impact parameter. The number of participant nucleons is a thus highly model dependent quantity, especially considering pA collisions.
Another way of studying possible nuclear effects in pA is to study particle production as a function of p ⊥ . In figure 14 we show a comparison to CMS data. The model is clearly not perfect, but nevertheless gives a fair description of the shape over the ten orders of magnitudes shown. Comparing to the results in [16] we again see an increased agreement due to the more careful treatment of secondary absorptive sub-events.

AA results
When we now turn to AA collisions, we expect the fluctuations to have less influence on the centrality measure, since at small impact parameters there are so many N N sub-collisions that most fluctuations will average out. It is therefore reasonable to assume that basically any centrality observable based on multiplicity or energy flow in the nuclei directions will be well correlated with the number of wounded nucleons and the actual impact parameter. Since we will now compare simulation to results from the ALICE experiment, we have in principle to use the ALICE experimental definition of centrality, rather than the one from ATLAS used in the previous chapter. In ALICE centrality is defined as percentiles of the amplitude distribution obtained in the two V0 detectors, placed at −3.7 < η < −1.7 and 2.8 < η < 5.1. Since this amplitude is not unfolded to particle level, and cannot be reproduced by Angantyr without realistic detector simulation, we instead construct a reasonable particle level substitute for this measure. We assume that the V0 amplitude is proportional to the total E ⊥ from charged particles with p ⊥ > 100 MeV in that region. In figure 15 we compare the measured V0 amplitude [59] with the substitute observable, scaled to match the bin just before the distribution drops sharply at high amplitudes. The shape of the distribution is described quite well, while the normalisation is a bit off. This is likely due to difficulties extracting the data for very low amplitudes. We will throughout this section use this as a centrality observable, combined with the trigger setup described in ref. [59]. Furthermore, all experiments have some definition of what a primary particle is. In figure 12 we used the ATLAS definition where all particles with cτ > 10 mm are considered as primary 13 . The ALICE definition is at its heart very similar, but has been described in more detail in  Figure 15: Scaled E ⊥ of charged particles at −3.7 < η < −1.7 and 2.8 < η < 5.1 from Angantyr, compared with the ALICE V0 amplitude, data taken from ref. [59]. ref. [60]. This definition has been conveniently implemented in Rivet [60], and we use this definition instead of a cut on cτ . In order to finish the discussion on the centrality measure, we show in figure 16(a) the ALICE results on the centrality dependence of the average charged multiplicity in the central pseudo-rapidity bin for PbPb collisions at √ s NN = 2.76 TeV [59] using the measured centrality, and in figure 16(b) with impact parameter bins. The agreement between these two results are clearly much better in PbPb than for pPb, confirming the initial statement in this section.
In figure 16(a) we also show our predictions 14 for Xenon-Xenon collisions at √ s NN = 5.44 TeV compared to the ALICE data that were published in [61].
In figure 17 we show the charged multiplicity compared to ALICE data [62][63][64] over a much wider η range, for both √ s NN = 2.76 TeV and √ s NN = 5.02 TeV. The trend, also visible in figure 16, is that Angantyr produces somewhat too few particles at central η; systematically the multiplicity is 5-10% too low. We regard this as surprisingly good, considered that no tuning of any kind to AA data has been done.
We now turn to transverse momentum spectra in AA collisions. In figure 18 we show results from ATLAS [65] compared to our model. The published p ⊥ spectra was scaled with the average number of wounded nucleons, calculated using a black disk Glauber model. We have not used the number of wounded nucleons as input to Angantyr, just scaled our result with the same number (as published in the article) to obtain comparable spectra. Hence, the results are not scaled to match, as both are simply scaled with the same number.
Finally we want to add a comment about the low multiplicity in the central region, shown in figs. 16(a) and 17. One of the main features of Angantyr is that tuning of MPI model, shower and hadronisation should only be carried out using e + e − , ep and pp data. However, looking at the comparison to pp in figure 10, we see that even the pp model undershoots the multiplicity at very low p ⊥ (below 500 MeV). Since ALICE measures charged particle multiplicity all the way down to zero transverse momentum 15 , it is not clear if the default PYTHIA8 behaviour should even be applicable here. The transverse momentum of such low-p ⊥ particles does not origin in the (perturbative) parton shower, but rather in the dynamics of string breakings. As seen from the comparison to pp this is not yet fully understood. The validity of this point is underlined by comparing to the ATLAS data shown in figure 18, where multiplicity is measured with low-p ⊥ cut-off of 500 MeV. In figure 19 we show the multiplicity distribution obtained by integrating the distributions measured by ATLAS, and see that the description improves. We want to make clear that (part of) this discrepancy could of course be due to a faulty comparison to data, where triggers, centrality measure etc. is not implemented in exactly the way as it is done by experiments. But if it is not, it points to an interesting point for improvement of the underlying model for soft particle production, also in pp. We will return to this subject in a future paper, but meanwhile we note that it would be interesting if experiments like ALICE, who can measure very near zero p ⊥ , will extend their publications to also include data with a minimum p ⊥ cut-off, which could serve as an important aide in further understanding.

Collectivity and non-flow estimation
One of the primary goals of the heavy ion programs at RHIC and LHC, is to investigate the collective behaviour of final state particles produced in collisions of nuclei accelerated to relativistic energies. The anisotropic flow measures the momentum anisotropy of the final state particles. As such, it is sensitive to both the initial geometry of the nuclear overlap region, as well as the transport properties of the final state before hadronisation.
The anisotropic flow is quantified in flow coefficients v n and corresponding symmetry planes Ψ n , defined by a Fourier series decomposition of the azimuthal distribution of final state particles: In practise, the flow coefficients are calculated using cumulants [66][67][68], which we also employ here. When flow coefficients are calculated using two-particle cumulants, the calculated coefficient also picks up azimuthal correlations not related to collectivity, but from  [70,71] (with ∆η = 1), compared to the non-flow contribution calculated by Angantyr. In the ratio plot it is seen that the non-flow contribution without ∆η-gap is nearly 40%. This is reduced to 20% when applying a gap. e.g. resonance decays and intra-jet correlations. Such "non-flow" effects can be suppressed by requiring a gap in η between particle pairs.
In figure figure 20 we show v 2 {2} as function of centrality 16 measured with and without a ∆η gap of 1.0, by ALICE [70,71] and CMS [69] respectively. Since Angantyr produces a full final state, it allows for the construction of the same observable, even in the absence of collective effects, giving an estimate of the non-flow present. We see that the non-flow contribution in the most central collisions is negligible (as one would expect), but rise to about 40% of the measured result for v 2 {2} without gap for peripheral collisions. This number falls to 20% when a gap is included, indicating that the method of applying a gap can remove some non-flow effects, but not all.
We want to emphasise that at this point, Angantyr does not make any attempt at modelling collective effects, and can therefore be used to estimate the contribution of nonflow. It is our plan to introduce a microscopic model for collectivity, based on string-string interactions to Angantyr, which has shown promising results in pp. The increased energy density from overlapping strings would here give a transverse pressure, leading to strings "shoving" each other before hadronisation [35,72].

Relation to other models
As Angantyr is a new model, it is instructive to compare it to existing models, and we here discuss the most commonly used ones, also mentioned in the introduction, HIJING [6], AMPT [5], and EPOS-LHC [4]. Here HIJING is most similar to Angantyr. Like Angantyr it is constructed as an extrapolation of pp dynamics, with the explicit motivation that differences between the model and experimental results may indicate effects of collective behaviour. In contrast AMPT and EPOS are both assuming collective expansion of a thermalised medium.
The HIJING generator is built with a similar starting point as Angantyr, thus it is inspired by the Fritiof model, using PYTHIA for generating multiple hard partonic subcollisions and the Lund string model in PYTHIA for the hadronisation. Similarly to Angantyr, HIJING relies on a Glauber calculation to determine the number of inelastic subcollisions, which are of two types: soft nucleon-nucleon collisions treated as in Fritiof, and hard parton-parton collisions treated as in PYTHIA. A new version written in C++ was recently presented [73].
In contrast to Fritiof the interacting nucleons are in HIJING excited to higher masses, covering most of the available rapidity range, but just as in the later Fritiof version [13,74], gluon radiation is added using the soft radiation scheme [75] implemented in Ariadne [76]. The hard partonic scatterings are determined via nucleus PDFs, where the parton density is suppressed by a shadowing factor R a/A , compared to A independent nucleon PDFs. To avoid double counting, emitted gluons in the soft component is allowed only for p ⊥ below a scale p 0 (chosen to be ≈ 2 GeV), while the hard partonic collisions have a lower cut at p ⊥ = p 0 .
Another difference between Angantyr and HIJING is that in HIJING fluctuations are neglected both in the initial states of the individual nucleons and in the position of nucleons within the nuclei. The soft N N amplitude is then chosen to reproduce the inelastic cross section including diffraction. The probability for multiple scattering is determined by the nuclear overlap function in impact parameter space. In Angantyr we find that fluctuations plus the distinction between primary and secondary absorptively wounded nucleons, have a quite significant effect for the final state multiplicity. In HIJING, the same effect may partly be due to the introduction of the shadowing factor R a/A . The shadowing factor is a geometry dependent "k-factor", which accounts for nucleons shadowing for each other during the collision, thus reducing to nucleon-nucleon cross section from the result obtained from pp collisions, to a lower, effective cross section. This suppresses the hard partonic cross section with up to 50% in AuAu collisions at √ s NN = 200 GeV [6]. In the end all partons are in HIJING connected by strings, and hadronised with PYTHIA. As an option it is possible to include a model for jet quenching, and also a jet trigger, enhancing the rate for events with high-p ⊥ jets. As mentioned above, AMPT presumes that a hot dense medium is formed. It uses the parton state obtained in HIJING as initial conditions. The partons then evolve in a partonic cascade up to freeze-out. After freeze-out the partons are connected in strings, which hadronise according to the Lund model in PYTHIA. Finally the obtained hadrons form a secondary cascade until the density is low enough, when they continue as free particles. As an option the hadronisation can also be calculated via quark-antiquark coalescence.
Finally the EPOS model works on different principles than the other two, as no explicit Glauber calculation is performed. Instead partonic sub-collisions are calculated using parton-based Gribov-Regge theory [77]. An elementary scattering is here represented by a cut Pomeron or "parton ladder". This ladder is interpreted as a flux tube, or a string, where the intermediate gluons provide a transverse motion. The strings then break up into segments by quark-antiquark pair production. In the central region with high density, the "core", the segments within a bin in η form a cluster, which expands longitudinally and radially until freeze-out. In regions of low density, called the "corona", the strings fragment instead directly to hadrons. This is mainly the case in the fragmentation regions. In a recent version, called EPOS LHC [4], a new flow parametrisation is introduced, which does not take advantage of the complete hydrodynamical calculation followed by the hadronic cascade as in EPOS2 [78] or EPOS3 [79]. One consequence is here that the time for one PbPb event is reduced from one hour to a few tenths of a second. According to the authors, this also implies that this version should not be used for a precise study of p ⊥ distributions or particle correlations in HI collisions.
In figure 21 we compare the multiplicity spectra at √ s NN = 2.76 TeV from figure 17(b) with Angantyr and the three generators discussed above. (For HIJING jet quenching is disabled in figure 21, but this should not have a major impact on the result.) We note that with all differences mentioned above, all four generators produce quite similar results for the centrality dependence of the charged particle distribution.
Comparing first HIJING to Angantyr, we see that while Angantyr undershoots at mid-η, HIJING overshoots on the full interval, and produces a too wide shape for the distribution. The likely source of this difference is the different way of handling secondary absorptive events, as described in section 3. HIJING treats all absorptive events on a similar footing, but the nuclear shadowing included in HIJING implies that it produces an overall lower amount of hard sub-collisions.
AMPT uses HIJING for initial conditions, but compared to HIJING the overall multiplicity is reduced by the partonic and the hadronic cascades. However, although the central density agrees with data, the distribution is too wide. We also note that AMPT reproduces multiplicity at mid-η better than Angantyr, and refer to our discussion about possible retuning of Angantyr to low-p ⊥ pp data in section 6.3. Finally EPOS-LHC also does a better job than Angantyr at mid-η, but worse away from the central region. We note that AMPT and EPOS-LHC, which both include the hydrodynamic expansion of a hot medium, do not describe data better than Angantyr over the full η-range.
It is, however, clear that if one wants to pin down the physics of a possible plasma phase, more exclusive observables than particle production must be used. This is indeed also the case in contemporary studies at the LHC and RHIC. Considering the precision obtained by the current tools, we see that there is a need for improved tools for comparing theory to data in heavy ion physics. To account for the final 10% discrepancy shown by all four generators, analysis specific effects like choice of centrality measure, trigger selection, primary particle definition etc. all play a major role. In the present paper all comparisons of Angantyr to data are carried out using the Rivet tool [58], which has proved highly successful for this task in pp. This has, however, been done using our own implementation of the experimental analyses. It is crucial for the further development of Monte Carlo event generators for heavy ion physics, that present and future heavy ion data is released using Rivet (or a similar tool), and we are pleased to note that experiments are now starting to commit to this task, also in heavy ion physics.

Conclusion and Outlook
We have introduced a new model called Angantyr for generating exclusive final states in proton-nucleus and nucleus-nucleus collisions. It extrapolates pp dynamics with a minimum of free parameters, and without assuming a hot thermalised medium. The aim is to see how well such an extrapolation can reproduce experimental data, thus exposing effects of collective behaviour. The model is a generalisation of the model for pA collisions in ref. [16], and is based on the following points: • The basic pp interaction is described by the PYTHIA8 event generator, based on multiple partonic sub-collisions and string hadronisation.
• The generalisation to nucleus collisions is inspired by the Fritiof model, and the notion of "wounded" or "participating" nucleons.
• The number of wounded nucleons is calculated from the Glauber model in impact parameter space, including Gribov corrections due to diffractive excitation of individual nucleons.
• Diffractive excitation is treated by the Good-Walker formalism, as the result of fluctuations. We here account for fluctuations in the states of individual nucleons, both in the projectile and in the target, and also fluctuations in the position of nucleons in the nuclei.
• The nucleon-nucleon cross section is parametrised, so as to reproduce the pp total, elastic, and inelastic cross sections, specifying diffractive excitation of the projectile, the target, and double diffraction, as well as the forward elastic slope.
• The model is implemented in an event generator, which generates exclusive final states. It is included in the PYTHIA8 package, where the user simply specifies a nucleus instead of a hadron as projectile and/or target. The possibility to add a signal process (of electroweak or other origin) is also included, enabling the user to study every process one could normally study in a pp collision.
We have shown that Angantyr does well in reproducing centrality dependent pseudorapidity distributions in pPb and PbPb collisions, at several energies, and have provided predictions for XeXe, in good agreement with ALICE data published later. We noted that the charged particle multiplicity at mid-rapidity is off by 5-10%, when comparing to measurements which also include low-p ⊥ particles. This problem is, however, not limited to nuclear collisions, though more visible here. In pp collisions PYTHIA8 has similar problems giving precise descriptions of multiplicities when p ⊥ approaches the soft region. An important future effort will include further development of the non-perturbative models describing the dynamics at this scale, as well as their interface to the perturbative shower.
Transverse momentum distributions in PbPb are not described equally well, especially for central collisions. The similar distributions in pPb are described slightly better, but still not at the level of precision known from pp collisions. We expect this to improve with the future inclusion of flow-like effects due to high density of strings. In ref. [34] we showed that overlapping strings forming "ropes" can reproduce the increased strangeness in pp [7], as well as in pPb and PbPb [80] collisions. In ref. [35] we further showed that the transverse pressure due to the increased energy density provides a transverse expansion and a qualitative description of the "ridge" observed in pp collisions. An important future direction will be to fully include these models in Angantyr, and test to what degree they provide a description of observed collective effects.
To conclude we think that it is notable that a direct extrapolation of pp dynamics can reproduce general features of inclusive particle production in AA collisions to better than 10%. It emphasises the importance of correlations, and in a future version of Angantyr we plan to include the collective effects from string-string interactions in the description of collisions with nuclei. In the future we also want to find observables sensitive to the fluctuations related to diffractive excitation and the internal substructure of nucleons. This is an essential feature which distinguishes Angantyr from other event generators available for nucleus-nucleus collisions.
Acknowledgments CB wants to thank C. H. Christensen for assistance in comparing to ALICE data. This work was funded in part by the Swedish Research Council, contracts number 2016-03291, 2016-05996 and 2017-0034, in part by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme, grant agreement No 668679, and in part by the MCnetITN3 H2020 Marie Curie Initial Training Network, contract 722104.

A. Generating absorptively and diffractively wounded nucleons
Here we will go through the technicalities of choosing the interactions between projectile and target nucleons. In [16] we showed that for a fixed nucleon-nucleon impact parameter, b, and a fixed projectile state, the cross section for the target nucleon to be wounded is given by the average of the fluctuations in the target nucleon. Writing the imaginary part of the scattering amplitude for given projectile and target states, p and t, in terms of the corresponding S-matrix, T pt (b) ≡ 1 − S pt (b), we have This works well for pA collisions, but for AA we also want to look at the probability for the projectile nucleon being wounded, and on top of this we want to be able to separate between absorptively and diffractively wounded nucleons.

A.1 Absorptively wounded nucleons
We expect the absorptively wounded nucleons will give the most important contributions to the final state particle production, and we therefore want to take special care to capture cross section fluctuations in this case and at the same time make sure we correctly reproduce the absorptive nucleon-nucleon cross section, The procedure will therefore be to generate one state for each nucleon in the projectile and target nuclei and for each pair of nucleons calculate and declare the nucleon-nucleon interaction absorptive with this probability. This will clearly give the correct absorptive nucleon-nucleon cross section. If we find the interaction is not absorptive we want to go on and check if either the target or the projectile or both are diffractively wounded, but this will then require us to consider averages over the possible states of the projectile or target or both. In the following we will consider a diffractively wounded target, but the corresponding treatment of the projectile is completely analogous.

A.2 Diffractively wounded nucleons
In general it is not necessarily straight forward to analytically calculate the average S 2 pt (b) t needed to get the the correct cross section for diffractive excitation. Instead we will estimate the fluctuations by generating a secondary, or auxiliary state for each projectile (p ′ ) and target (t ′ ) nucleon. We will still calculate the probability for absorptive interaction using only the primary states, but to get the probability of the target nucleon to be wounded we note that the product S pt (b)S pt ′ (b) will on average yield the correct value for S 2 pt (b) t p , so naively we could use the probability P Wt = 1 − S pt (b)S pt ′ . However, it is clear that we will then have a negative probability for having a diffractively wounded target P Wt − P abs < 0, for S pt < S pt ′ . Therefore we also need to consider the statistically equivalent situation where the absorptive interaction probability is given by while the corresponding wounded probability is still where the probability for a diffractively wounded target is then positive. The procedure we have chosen to handle this, is to shuffle probabilities between the two situations so that we always get non-negative probabilities for diffractively wounded nucleons according tõ P Wt = S tp < S pt ′ : 0 S pt > S pt ′ : P Wt + P ′ Wt − P ′ abs = 1 − 2S pt (b)S pt ′ (b) + S 2 pt ′ (b) (A.6) P ′ Wt = S tp < S pt ′ : P ′ Wt + P Wt − P abs = 1 − 2S pt (b)S pt ′ (b) + S 2 pt (b) S pt > S pt ′ : 0 (A.7) which will give the correct cross section for the target nucleon being wounded. By considering the auxiliary state for the projectile, p ′ , we can then also find the probability for the projectile being diffractively wounded. And if both are wounded we say that the interaction is a double diffractive excitation 17 .