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 8 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 JHEP10(2018)134 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 with a minimum of adjustable parameters. In this way it forms a bridge between heavy ion and high energy hadron phenomenology. Angantyr is a generalisation to AA collisions of the model for pA scattering in ref. [9], which 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. Like PYTHIA8 and the model in ref. [9], Angantyr does not include an assumption of a hot thermalised medium. The model can therefore serve as a baseline for understanding the non-collective background to observables sensitive to collective behaviour.
Before discussing the generalisation to heavy ion collisions, we want to discuss some features of high energy pp scattering, which are important for this generalisation.
First, as will be discussed in more detail below, diffractive excitation is important. At high energies the real part of the pp amplitude is small, and usually neglected in applications to collisions with nuclei. Diffraction (elastic scattering and diffractive excitation) is then the shadow of absorption into inelastic (non-diffractive) 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 (the proton) is a (coherent) linear combination of scattering eigenstates with different absorption probability. These eigenstates have in refs. [10,11] been interpreted as different parton cascades.
Secondly multiple partonic sub-collisions are very important at high energies. Here we use the scheme from ref. [12], as implemented in PYTHIA8, to describe inelastic nondiffractive events. Hard scattering is also seen in diffractive events, and here we use the Ingelman-Schlein formalism [13], which is also included in the PYTHIA8 package.
A generalisation of the formalism for pp collisions to an event generator for pA and AA collisions will have four separate components: (i) It is necessary to determine nucleon positions within the colliding nuclei. Here a number of MCs are already available to generate nucleon distributions, see e.g. refs. [14][15][16][17].
(ii) One has to calculate the number of interacting nucleons and binary N N collisions. This is generally performed using the Glauber formalism [18,19]. 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 [20], but has often been neglected also in recent applications (see e.g. the review JHEP10(2018)134 by Miller et al. [19]). 1 As mentioned above, diffractive excitation is a consequence of fluctuations in the nucleon substructure. An important point is then that a nucleon in the projectile is fixed in the same state during its passage through the target nucleus. (And similarly the state of a target nucleon is fixed through the projectile nucleus.) Fluctuations in the projectile proton in pA collisions was studied by Heiselberg et al. [21], for estimates of the number of individual N N sub-collisions. This formalism was further developed in several papers (see refs. [22][23][24][25] 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. [26,27].
As discussed in ref. [9], 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 take into account individual fluctuations in both projectile and target nucleons. As far as we know, Angantyr is the first model where this condition is satisfied.
(iii) One must 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 [28,29] and the notion of "wounded" nucleons. 2 Bia las, Bleszyński, and Czyż [30] showed that the production of soft particles is determined by the number of wounded (or participant) nucleons, rather than the number of individual N N sub-collisions. (The latter was later seen to be correlated to hard processes, like production of high p ⊥ particles or vector bosons.) In the early Fritiof model [28] 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 later obtained by Bia las and Czyż in an analysis of dAu collisions at RHIC [31].
The Fritiof model did not explicitly include diffractive excitation. We note, however, 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. [31]. The wounded nucleons in Fritiof can therefore effectively represent both non-diffractively and diffractively wounded nucleons.
(iv) At high energies, the hard partonic sub-collisions (scaling with N N sub-collisions rather than wounded nucleons) play a very essential role. It is therefore necessary to account for those specifically in events with multiple N N collisions, e.g. when one projectile nucleon interacts with several target nucleons (or vice versa). In ref. [9] 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 JHEP10(2018)134 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. exclusive final states in AA collisions, we then have to calculate all sub-collisions 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 [32] 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. [33,34] 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 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 JHEP10(2018)134 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 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. In section 7 we discuss model uncertainties, especially related to our treatment of secondary absorptive sub-collisions. Finally we discuss differences and similarities between our approach and other heavy ion event generators in section 8, before presenting some conclusions and an outlook in section 9.
2 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 inelastic cross section is then simply the difference between the two. 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. [18]. 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 [35], as the result of fluctuations in the nucleon wave functions. 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.

JHEP10(2018)134
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. [14,15]), or a more sophisticated description of the two-(or three-) particle correlations between the nucleons within the nucleus [16,17]. Fluctuations in the geometry is taken into account, when new nucleus states are generated for each new event.

Fluctuations in the individual NN 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 :

JHEP10(2018)134
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.

JHEP10(2018)134
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 1: no fluctuations. The N N amplitude and S matrix are T NN = 1/2 and S NN = 1 − T NN = 1/2. The inelastic cross section when hitting two target nucleons is then from eqs. (2.2), (2.7) given by dσ inel /d 2 b = 15/16, (σ D = 0). 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. 3

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

JHEP10(2018)134
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. 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. [19]). 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 σ 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 [21][22][23][24][25]. It is there assumed that the fluctuations in the projectile can be described by a distribution in the quantity σ ≡ d 2 b 2T (b) t of the form

JHEP10(2018)134
(The second relation follows from eq. (2.6).) This formalism, has been used in analyses of pPb data from LHC, e.g. by ref. [26], 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. [23] 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. [9] we investigated the fluctuations in the nucleon cross sections using Mueller's dipole approach to BFKL evolution [36,37] as implemented in the DIPSY Monte Carlo program [11,38,39]. The model is formulated in impact parameter space, and includes also a set of sub-leading corrections beyond the leading-log BFKL approximation. Nonlinear 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 [21], 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 [9] we then also obtained a fair description of e.g. the observable used by ATLAS in [26] 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. [40]. 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. [40] 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 two more parameters, σ t and α, (besides k and r 0 in eq. (2.11)) and 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 P abs = 2T ik − T 2 ik , with T ik = T (b, r iµ , r kν ) 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. [9] 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 (cf. 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. ing 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%.

JHEP10(2018)134
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 [41,42] does not necessarily agree well with cross section measurements from LHC [43], 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.

From wounded nucleons to exclusive final states
In the wounded nucleon model, as formulated by Bia las and Czyz [31], 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 JHEP10(2018)134 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.
model, F (η) must be extracted from data, and depends on centrality class [44], 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 [12], all partonic sub-collisions are to a first approximation treated as separate QCD 2 → 2 scatterings. 4 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 function using some assumption about the matter distribution in the colliding protons and an assumed impact parameter.

JHEP10(2018)134
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 [12] 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 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 collision and imagine the projectile 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. A secondary absorptive wounded nucleon thus contributes to the final state as if the final state particles were produced in a single diffractive excitation. This similarity is what we, in the following, JHEP10(2018)134  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. will exploit to build up a final state from primary absorptive interactions and secondary absorptive interactions, the latter being modelled as single diffractive excitation.
The procedure will therefore be to decide which of 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 [9] 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 JHEP10(2018)134 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 [9] 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.

JHEP10(2018)134 4 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 [9] and detailed below.

Selecting primary absorptive collisions
The first pass over the potential N N interactions, we will only look at absorptive interactions. 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 [9] 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 [45] or top events in pAcollisions [46] 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

JHEP10(2018)134
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. 5 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 [47] 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 [48] that this could give rise to effects at the 10% level in selected observables in peripheral PbPb collisions. It has also been pointed out [49] 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 [50], and in the present version we have left them out entirely. 6

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 [9] we will generate these secondary absorptive sub-collisions as if they were single diffractive excitation events. We here use 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 posteriori way. The 5 For most use cases this should be adequate, as σ NN signal σ NN abs for most processes of interest 6 An interested user can, however, plug in their own Glauber MC including neutron skin effects.

JHEP10(2018)134
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 diffractive interactions
Having taken care of all absorptive interactions we continue with diffractive interactions in much the same way. For each type we again go through the impact-parameter-ordered list of N N interactions twice. In the first round, we only consider primary interactions, i.e. where neither of the nucleons have previously been included in a sub-event, and generate a sub-event which could be a single or a double diffractive excitation. These are treated as (soft) diffractive events in PYTHIA8, as discussed in section 5.
In the second round we also consider secondary interactions, where one of the nucleons has already been treated, and an appropriate contribution from the other nucleon (which we will here call a half event) is generated and added to the corresponding previous sub-event.
As an example consider an already wounded nucleon in the projectile nucleus, which interacts with a previously unwounded nucleon in the target. The wounded nucleon is already connected to another target nucleon, and cannot be further excited. There are three possibilities for the diffractive interaction: 1. The new interaction is a single diffractive excitation of the target nucleon. The interaction is then treated as a normal single diffractive excitation of the target nucleon.
2. The new interaction is a single diffractive excitation of the projectile (already wounded) nucleon. In this case the target nucleon is elastically scattered.
3. The new interaction is a double diffractive excitation. In this case the already wounded projectile nucleon is not modified, and the interaction is again treated as a single diffractive excitation of the target nucleon.
In a final iteration 7 also purely elastic interactions are considered, and here again the half events are single elastically scattered nucleons. In each case energy and momentum conservation is handled in the same way as for secondary absorptive interaction.
Modulo the effects of secondary interactions being discarded due to energy-momentum conservation, this procedure will correctly handle the probability that a given nucleon is wounded in some way. Note however that, as discussed in section 2.5, although some nucleons in the program are classified as elastically scattered, elastic scattering is not included properly. As elastic scattering is a coherent effect of shadowing due to absorption, the Good-Walker formalism can be used to calculate the cross section for elastic scattering of the incoming nuclei, but not for individual nucleons in a nucleus. 8 Diffractive excitation of individual nucleons can, however, be calculated via the trick described in section 2.5.
In the end we have generated a set of parton-level sub-events, which we now can join together in a single parton-level AA event. This event is then handed back to PYTHIA8 for Charged particle η at 7 TeV, track p ⊥ > 500 MeV, for N ch ≥ 1 η 1/N ev 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 [51]. The black and red lines are the shapes of standard non-diffractive and single diffractive events from PYTHIA8 respectively. The green dashed and blue dash-dotted lines are single diffractive events generated by PYTHIA8 using the modifications presented in [9] 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 distributioǹ a la Fritiof.
hadronisation and decay of unstable hadrons. Finally the non-interacting projectile and target nucleons are bunched together in two remnant nuclei. 9 5 Modifications of single diffractive to secondary absorptive In section 3 and in ref. [9] we argued that secondary absorptive interactions will contribute to particle production in the same way as a single diffractive (SD) excitation event (cf. 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 [13]), 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,

JHEP10(2018)134
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. [9] 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. 10 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. 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 pIP 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 pIP cross section, σ pIP 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 TeV. 11 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 [9] 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 11 The kinematics is given by the LHC pPb run, giving a slightly tilted distribution in η.

JHEP10(2018)134
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 σ pIP 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 , 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 (cf. 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 figures 7 and 8 we show some distributions in the slices η ∈ [−5, −4], [−3, −2] and [−1, −0] (the shaded regions 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 gives a slightly better description of p ⊥ in figure 7 and the average p ⊥ observables in figure 8, but no improvement -or even a slightly worse performance -in the two remaining observables. The choice of which option to use can therefore only be based on an assessment of what types of observables are deemed most important to reproduce correctly. 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 [12], 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 ⊥ distributions in SD(glu) in figure 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).
It is, however, clear that we could have put more emphasis on charged multiplicity and p ⊥ in the regions where the SD(glu) option outperforms SD(new), and thereby made another choice of recommended option. In section 7 we will compare the three different choices against each other for pA results.
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 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 [9]. 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.

JHEP10(2018)134 6 Sample results
All results presented here are generated with PYTHIA8 version 8.235 using default settings. 12 This means in particular that: • the nucleon distributions in the nuclei are generated according to the formulae in [15] 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. 13 Below in section 7 we will study the effects of these choices.

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 12 Since Angantyr is the default heavy-ion model in PYTHIA, it suffices to specify suitable nuclei as beam particles to reproduce the results presented here. 13 The number of attempts allowed for this is governed by the parameter Angantyr:SDTries. Charged particle η at 7 TeV, track p ⊥ > 100 MeV, for N ch ≥ 2  is not the same, and since this directly affects the amount of MPI it is important to make sure that 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 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,

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 [26], 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.
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 therefore show two sets of lines generated with Angantyr with the two different binnings presented in figure 11. Clearly the difference between the two is 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 14 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 [9], 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 of the nucleus.
To study the fluctuations further we show in figure 13 the average number of wounded nucleons as a function of E P b ⊥ -centrality, both for Angantyr and for three Glauber-model 14 The η-distributions in figure 12(a) has been corrected for detector effects. fits performed by ATLAS in [26]: 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 ref. [9] 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 must 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. [53].
principle use the ALICE experimental definition of centrality, rather than the one from ATLAS used in the previous section. 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 [53] 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. [53]. 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. 15 The ALICE definition is at its heart very similar, but has been described in more detail in ref. [54]. This definition has been conveniently implemented in Rivet [54], 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 [53] 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 16 for Xenon-Xenon collisions at √ s NN = 5.44 TeV compared to the ALICE data that were published in [55].
In figure 17 we show the charged multiplicity compared to ALICE data [56][57][58] 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 15 This means e.g., that a pair of π + π − which comes from the decay of a K 0 S , will not be included in the charged multiplicity. 16 Although we present this after the data was published we still consider it a prediction, as the program was released before the data was analysed.

JHEP10(2018)134
Pythia8/Angantyr ALICE PbPb S NN = 5.02 TeV η; the multiplicity is systematically 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 [59] 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 figures 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, 17 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 [60][61][62], 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 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 20 we show v 2 {2} as function of centrality 18 measured with and without a ∆η gap of 1.0, by ALICE [64,65] and CMS [63] 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 [34,66]. (without ∆η-gap) and ALICE [64,65] (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.

Model uncertainties
The main idea behind Angantyr is to extrapolate pp dynamics, as described by the model for MPIs/underlying event in the PYTHIA8 MC, to heavy ion collisions, retaining as much as possible from pp. This principle was outlined already in the introduction, especially figure 1, but as the model has now been presented, as well as results from pA and AA collisions, we will here also discuss the model uncertainties related to this extrapolation procedure. Primary interactions correspond directly to inelastic non-diffractive pp collisions. Here PYTHIA8, is known to reproduce most features of both soft and hard pp collisions at LHC fairly well, and the extrapolation to primary interactions in a heavy ion collision is therefore mainly a source of model uncertainty up to PYTHIA8's shortcomings in describing such collisions in pp. We already discussed some of those shortcomings in the previous section, but as they are not uncertainties directly related to the Angantyr model (but rather the underlying PYTHIA8 model) we will not discuss them further here.
The largest uncertainty comes instead from our treatment of secondary absorbed nucleons. The main reason is that secondary absorption has no pp equivalent. In section 5 we outlined the procedure of modifying single diffractive collisions to describe secondary absorbed nucleons, and we will investigate uncertainties related to this treatment in section 7. of nucleons can in principle be determined from pp collisions, but as we will discuss in section 7.2, this is not straight forward.

Uncertainties in treating secondary wounded nucleons
A core feature of the Angantyr model, is that the contribution from a secondary absorbed nucleon is similar to the contribution from an excited nucleon in a single diffraction event. This corresponds to the black pieces in figures 3a and 3b respectively. This assumption has two components: (i) The distributions in the rapidity range covered, ∆y, and the corresponding mass, M ≈ exp(∆y/2) × (1GeV), are similar.
(ii) The distribution of partons from the projectile nucleon, involved in the interaction with the secondary absorbed nucleon in figure 3a, is similar to the partons in the Pomeron in figure 3b.
Naturally none of these assumed similarities can be exact. Extracting the relevant properties in diffractive excitation in pp collisions from data at LHC has also large uncertainties, as we will discuss further in section 7.2. We also note that: (iii) Energy-momentum conservation has generally important effects in high energy reactions, and has to be satisfied when nucleons suffer multiple N N sub-collisions.
Also this point is associated with some model uncertainty, as discussed in section 4.2 and in section 7.1.3 below. In the following we will discuss the uncertainties associated with all three choices in the treatment of secondary wounded nucleons, and their impact on model predictions. We will focus on pA collisions, where there can at most be a single primary interaction, and the treatment of secondary interactions consequently has a relatively larger effect. This is illustrated in figure 21, where we see that secondary absorbed nucleons correspond to about 80% of all wounded nucleons in central pPb collisions, but only about 25% in central PbPb collisions.

Mass distribution
We begin by discussing point (i), the mass distribution of secondary wounded nucleons. The picture in figure 3a has the structure of a triple-Pomeron diagram. This similarity is somewhat symbolic, as each chain in this figure includes the multiple parton scatterings in figure 2, which correspond to Pomeron loops in a Reggeon field theoretical approach (see e.g. refs. [67][68][69]). The triple-Pomeron diagrams shown in figure 22 would have a weight proportional to: for secondary absorption. we see a noticeable effect already below 50 GeV. The effect follows the expectation that a larger ∆ will give larger M A values and thus more activity. However, we also see that above 50 GeV the distributions for larger ∆ seem to run out of steam, which we attribute to the fact that higher M A values mean that the energy available from the projectile proton is used up faster. This means that fewer secondary absorptive interactions are accepted.
In figure 23b we also show the resulting pseudo-rapidity distribution of charged particles for two centrality bins (using the experimentally determined bin edges in E P b ⊥ ). The larger values of M A are also reflected in the η-distributions, where the effect is that the distribution becomes too flat to describe data, especially for central events.

Parton distribution in the projectile
As discussed in section 5, the secondary absorptive interaction in figure 3a may involve several partons coming from the projectile nucleon, in a way similar to how diffractive excitation is described by a Pomeron PDF in the Ingelman-Schlein model. Point (ii) concerns the distribution of these partons. In section 5 we studied three different distributions, SD(new) (which is the default for secondary absorption), SD(def) (which is the PYTHIA8 default for diffractive excitation), and SD(glu) (which is the modified PDF for increased gluon activity introduced in ref. [9]). In figure 24a we show the effect on the the differences found in section 5. The resulting pseudo-rapidity distributions shown in figure 24b do not show so dramatic differences. It is, however, clear that our default choice gives the best description of data. As discussed in section 5 our default choice is the one that makes most sense on theoretical grounds, and it is satisfying to see that it also makes sense in comparison to data.

Energy-momentum conservation
Energy-momentum conservation is frequently seen to have a very large impact in high energy reactions. Here its effect could be seen in figure 23a. It is not clear from first principles if energy-momentum conservation should prohibit a sub-collision, if a single sampling of the M A distribution turns out to require more than what is available, or if it is possible to simply try again. To further study the effects of this ambiguity, we show in figure 25, what happens if we allow Angantyr to retry adding secondary sub-events, which 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 in the E P b ⊥ centrality measure, while the effect on the resulting η-distribution is barely visible. It is interesting to note that the effect of allowing more attempts seem to saturate quickly, and going from 2 to 4 attempts makes a much smaller change than allowing two attempts instead of one (which is the default).

Diffractively excited nucleons
In contrast to the secondary absorbed nucleons, a positive ∆ in eq. (7.1) means lower masses for diffractively excited nucleons. In principle the M D -distribution could be measured in pp collisions at LHC, but it is quite challenging to isolate single diffraction from the experimental distribution in the size of a gap in rapidity (see refs. [70,71]).
In collisions with nuclei, multiple N N interactions imply that, the probability for absorption is enhanced and, as a consequence, the probability for diffractive excitation is  Figure 25. Same as figure 23, but now varying the number of attempts (parameter Angantyr:SDTries) allowed to generate secondary sub-events that can be added without violating energy-momentum conservation before giving up and vetoing the secondary interaction. The default version allows only a single attempt and is shown as the blue lines, while allowing two or four attempts is shown as dashed red and dotted green lines respectively. reduced. From figure 21 we see that in pA collisions about 20% of the wounded nucleons are diffractively excited, dropping to 10% in central pPb collisions. In AA collisions this fraction is further reduced to an average about 10%, and below 4% for central PbPb collisions. This implies that a reasonable variation of the diffractive component will have comparatively small effect. For this reason we have here chosen to keep the default setting in the PYTHIA8 MC, with a distribution ∝ dM 2 D /M 2 D . One could here also imagine including a Reggeon contribution ∝ dM 2 D /(M 2 D ) 1.5 . This contribution is concentrated to low masses, and would not affect the results in most of the rapidity range, including the forward detectors used to measure the centrality. It could, however, give a contribution in the very forward region, and thus it might be of importance e.g. for interactions with cosmic rays.

Uncertainties in AA collisions
Above we have discussed model uncertainties in pA collisions. We have also pointed out that the corresponding uncertainties are significantly smaller in AA collisions, in particular in central AA collisions. In figure 21 we showed that the fraction of wounded nucleons which are secondary absorbed is about 70% in pPb but about 35% in PbPb collisions. For central collisions these ratios are about 80% in pPb and only about 25% in PbPb. We have checked that a corresponding reduction of the uncertainties is obtained in the MC results for AA collisions.

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],

JHEP10(2018)134
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 [72].
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 [29,73], gluon radiation is added using the soft radiation scheme [74] implemented in Ariadne [75]. 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 [76]. 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 [77] or EPOS3 [78]. 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 26 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 26, 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 JHEP10(2018)134 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 [52], 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 in this way it bridges the gap between heavy ion and high energy physics phenomenology. It does not assume a hot thermalised medium, and 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. [9], 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.

JHEP10(2018)134
• The Glauber model is formulated in impact parameter space. Diffractive excitation is then most conveniently described by the Good-Walker formalism, as the result of fluctuations in the nucleon substructure. We here for the fist time account for fluctuations in both the projectile and the target nucleons, in a Glauber calculation.
(As frequently in MC simulations, fluctuations in the position of nucleons in the nuclei are also included.) 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 gives a good description of general final state properties. This includes not only multiplicity and transverse momentum distributions both in pPb and PbPb collisions, but also its dependence on centrality. We note, however, that this dependence is very sensitive to the experimental definition of centrality. Thus we see that for low centrality the correlation between central multiplicity and "centrality" is more a correlation between central and forward activity, rather than between central activity and impact parameter. The model predictions for XeXe collisions are also in good agreement with ALICE data published later.
The model underestimates somewhat central particle production, when p ⊥ is integrated down to p ⊥ = 0. This may be not surprising, as it is an extrapolation of PYTHIA's description of pp dynamics, which is too low for small p ⊥ below 200 MeV. Future work is needed to improve the hadronisation models in this region, including their interface to the perturbative shower.
The description of data is quite sensitive to the handling of, in particular, secondary absorptive sub-events. We have investigated several different choices relating to this treatment, relating to (i) the distributions in the covered rapidity ranges, (ii) distributions of partons in the projectile nucleon, and (iii) energy-momentum conservation. For visualization we performed this investigation in pPb collisions, noting that they will be significantly smaller in PbPb collisions. Although our final choices may not be based on completely solid theoretical grounds, the fact that alternatives investigated give a poorer description of data tells us, that the choices are reasonable. Certainly there are other variations to investigate, but we postpone such studies to a future publication.
In PYTHIA8 all strings decay into hadrons independently. Thus it does not include a mechanism to reproduce the collective effects seen in pp collisions. Such effects are therefore also not reproduced by the present version of Angantyr, and the model should be thought of as a baseline for understanding the non-collective background to observables sensitive to collective behaviour.
Also in high energy pp collisions the number of strings is quite large, in particular in events with high multiplicity. In ref. [33] we showed that overlapping strings forming "ropes" can qualitatively reproduce the increased strangeness in pp [7], as well as in pPb and PbPb [79] collisions. In ref. [34] we further showed that the transverse pressure due to JHEP10(2018)134 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 the observed collective effects in nucleus collisions. Besides the angular correlations, the transverse expansion may affect the p ⊥ distributions, which are less accurately reproduced in pPb and PbPb collisions.
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%. This emphasises the importance of correlation studies, 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.

A Generating absorptively and diffractively wounded nucleons
Here we will go through the technicalities of choosing the interactions between projectile and target nucleons. In [9] 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 JHEP10(2018)134 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 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.

JHEP10(2018)134
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. 19 Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.