Revealing a deep connection between factorization and saturation: New insight into modeling high-energy proton-proton and nucleus-nucleus scattering in the EPOS4 framework

It is known that multiple partonic scatterings in high-energy proton-proton ($pp$) collisions must happen in parallel. However, a rigorous parallel scattering formalism, taking energy sharing properly into account, fails to reproduce factorization, which on the other hand is the basis of almost all $pp$ event generators. In addition, binary scaling in nuclear scatterings is badly violated. These problems are usually ``solved'' by simply not considering strictly parallel scatterings, which is not a solution. I will report on new ideas (leading to EPOS4), which allow recovering perfectly factorization, and also binary scaling in $AA$ collisions, in a rigorous unbiased parallel scattering formalism. In this new approach, dynamical saturation scales play a crucial role, and this seems to be the missing piece needed to reconcile parallel scattering with factorization. From a practical point of view, one can compute within the EPOS4 framework parton distribution functions (EPOS PDFs) and use them to compute inclusive $pp$ cross sections. So, for the first time, one may compute inclusive jet production (for heavy or light flavors) at very high transverse momentum ($p_{t}$) and at the same time in the same formalism study flow effects at low $p_{t}$ in high-multiplicity $pp$ events, making EPOS4 a full-scale ``general purpose event generator''. I discuss applications, essentially multiplicity dependencies (of particle ratios, mean $p_{t}$, charm production) which are very strongly affected by the saturation issues discussed in this paper.

It is known that multiple partonic scatterings in high-energy proton-proton (pp) collisions must happen in parallel.However, a rigorous parallel scattering formalism, taking energy sharing properly into account, fails to reproduce factorization, which on the other hand is the basis of almost all pp event generators.In addition, binary scaling in nuclear scatterings is badly violated.These problems are usually "solved" by simply not considering strictly parallel scatterings, which is not a solution.I will report on new ideas (leading to EPOS4), which allow recovering perfectly factorization, and also binary scaling in AA collisions, in a rigorous unbiased parallel scattering formalism.In this new approach, dynamical saturation scales play a crucial role, and this seems to be the missing piece needed to reconcile parallel scattering with factorization.From a practical point of view, one can compute within the EPOS4 framework parton distribution functions (EPOS PDFs) and use them to compute inclusive pp cross sections.So, for the first time, one may compute inclusive jet production (for heavy or light flavors) at very high transverse momentum (p t ) and at the same time in the same formalism study flow effects at low p t in high-multiplicity pp events, making EPOS4 a full-scale "general purpose event generator".I discuss applications, essentially multiplicity dependencies (of particle ratios, mean p t , charm production) which are very strongly affected by the saturation issues discussed in this paper.

I. SOME INTRODUCTORY REMARKS ABOUT FACTORIZATION, PARALLEL SCATTERING, AND ENERGY SHARING
Two major discoveries made it possible to reliably compute cross sections in high-energy proton-proton (pp) scattering.There is first of all the fact that the coupling constant α s of strong interactions becomes weaker with increasing scale, referred to as "asymptotic freedom" [1,2], which allows the use of perturbation theory to compute parton-parton cross sections.The other crucial issue is called "factorization" [3,4], which amounts to separating short-and long-distance physics at some "factorization scale" µ, which allows one to write the inclusive pp cross section as a convolution of two parton distribution functions (PDFs) and a (calculable) elementary parton-parton cross section.The PDFs contain all the long-distance physics, below the scale µ.Factorization in connection with asymptotic freedom turned out to be an extremely powerful concept, with numerous important applications.Extended to collisions of two nuclei, composed of A and B nucleons, factorization means that the cross section for rare processes is given as AB times the pp cross section.This is usually referred to as "binary scaling".
Factorization is an impressive tool, being very useful when it comes to studying inclusive particle production, but there are very interesting cases not falling into this category, like high-multiplicity events in proton-proton scattering in the TeV energy range, where a very large number of parton-parton scatterings contribute.Such events are particularly interesting, since the CMS Collaboration observed long-range near-side angular correlations for the first time in high-multiplicity protonproton collisions [5], which was before considered to be a strong signal for collectivity in heavy ion collisions.And studying such high-multiplicity events (and multi-plicity dependencies of observables) goes much beyond the frame covered by factorization.Here one needs an appropriate tool, able to deal with multiple scatterings.
The most important fact about multiple parton-parton scatterings is that they must occur in parallel, and not sequentially, as I am going to justify in the following.It is known that parton-parton scatterings are preceded by a series of successive parton emissions according to Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations [6][7][8].In particular, the first emitted partons carry a large momentum corresponding to a large γ factor, so they are "long-lived" particles.Correspondingly, the whole reaction takes a long time, which makes it impossible to have two (or more) successive parton-parton scatterings.Multiple scattering must therefore happen in parallel.In the case of nucleusnucleus scattering, the nucleon-nucleon collisions also happen in parallel, and this is simply due to the fact that at very high energies, the "reaction time" (the time it takes for the two nuclei to pass through each other) is much shorter than the particle formation time.So first all the interactions are realized (instantaneously) and particle production comes later.One has a "double parallel scattering" scenario: the nucleon-nucleon scatterings happen in parallel, and for each nucleon-nucleon scattering, the parton-parton collisions occur in parallel.
In the case of multiple scatterings, energy-momentum conservation is an important issue.Of course, everybody agrees on that and all event generators do conserve it.But this has to be seen in the light of the underlying theory.Multiple scatterings have been incorporated in an S-matrix approach in the Gribov-Regge (GR) theory a long time ago [6,[9][10][11].All the scatterings "are equal" in the sense that there is no sequence, nothing like "a first scattering" and a "second scattering" and so on.However, as discussed in [12,13], there is an inconsistency: the "energy-momentum sharing" is simply not taken into account.In a strictly parallel scenario, the initial energy-momentum has to be shared among n parallel scatterings and a projectile and a target remnant in an unbiased way, in the following referred to as "rigorous parallel scattering scenario", which amounts to using integrands (for pp scattering) like (with p referring to four-momentum).In the case of AA scattering, one has products of δ-functions of the type Eq.(1).I insist on the fact that the theoretical basis (S-matrix theory) on one hand, and the Monte Carlo realization in event generators on the other hand, should deal with energy conservation in the same way, i.e., they should be 100% compatible with each other; this is also what I mean by a "rigorous parallel scattering scenario".This is usually not the case, as in some early work of the author [14], where the underlying theory has no energy sharing, but the Monte Carlo realization does, and this is even in recent event generators the usual method.
Employing the "rigorous parallel scattering scenario", one encounters highly multidimensional integrals that cannot be separated.In [15], nevertheless energy sharing in the sense of Eq. ( 1) and its generalization to AA scattering could be implemented, and the technical difficulties (using Markov chain techniques) could be handled.I am not aware of any other attempt in this direction.
Let me discuss the fundamental differences between the "standard QCD generators" and the "rigorous parallel scattering scenario".All "standard QCD generators" such as PYTHIA [16], HERWIG [17], or SHERPA [18], take as starting point the factorization formula, sketched in Fig. 1.In this plot and all the following ones, I show for simplicity only gluons, in the real calculations all kinds of partons are considered.The two light blue thick lines represent the projectile and the target protons.The proton structure and the so-called space-like parton cascade are taken care of by using parton distribution functions (PDFs) f , which allows writing the jet cross section as a convolution of these PDFs and an elementary QCD cross section for the Born process in the middle.This formula (still based on Fig. 1) serves as a probability distribution, which allows one to generate a sequence of hard processes, which are ordered in "hardness".This is the method to introduce multiple parton scattering.
In the "rigorous parallel scattering scenario", the starting point is a multiple scattering diagram as shown in mathematical formula contains in the case of pp scattering a δ function as in Eq. ( 1) for energy-momentum conservation.Here one also considers parton evolutions from both sides, but for each of the n interactions, so one cannot use the usual proton PDFs.Instead, one considers 2n evolutions, starting always from the first perturbative parton on the projectile side and the target side.But one is nevertheless able to define evolution functions E, which are based on the same DGLAP partial differential equations (see for example [3,4]), but in our case, the initial condition is not a parton distribution f (Q 2 0 , x) in the proton at some initial scale Q 2 0 , but a parton carrying the full momentum fraction x = 1.The Monte Carlo procedure to generate partons will be done in two steps: • Step 1: The multi-scattering formalism allows generating a number of scatterings n and in addition for each of the n scatterings its energy (expressed in terms of light-cone momentum fractions x ± i ), with 100% energy-momentum conservation (the cross section formulas contain a δ function to assure it, so energy-momentum violating configurations are never proposed).
• Step 2: With n and all the x ± i known, one generates for each of the n scatterings the hard process based on a convolution E proj ⊗ Born ⊗ E targ , and then the parton emissions via backward evolution.
The technical problems in the "rigorous parallel scattering scenario" can be handled, but there are conceptual problems.In the classical GR approach, it is known that in the case of inclusive cross sections, the multi-scattering contributions cancel, referred to as Abramovsky-Gribov-Kancheli (AGK) cancellations [11], and considering a Pomeron to be a parton ladder, one may deduce factorization and binary scaling in AA, as discussed in more detail in Sec.IV.
In the parallel scattering scenario with energymomentum sharing, imposed in an unbiased fashion via a delta function as in Eq. (1), one does not get factorization (which requires a single Pomeron contribution), and one violates terribly binary scaling for AA scattering, as I am going to discuss in Sec.VI.
The solution to the problem is related to the treatment of saturation, as I discuss in a very qualitative fashion in Sec.II, and in detail in Sec.VII.

II. SOME INTRODUCTORY REMARKS ABOUT SATURATION AND ITS RELATION WITH ENERGY SHARING
The above sketched "rigorous parallel scattering scenario" is an elegant way to introduce unbiased parallel scattering, but in the end, it does not work, one violates factorization (and binary scaling in AA).So something is still missing.
There is actually another important issue in highenergy scattering: with increasing energy, partons with very small momentum fractions x ≪ 1 become increasingly important, since the PDFs at small x become large.This means that the parton density becomes large, and therefore the linear DGLAP evolution scheme is no longer valid, and nonlinear evolution takes over, considering explicitly gluon-gluon fusion.These phenomena are known as "small x physics" or "saturation" [19][20][21][22][23][24][25][26][27][28][29][30][31][32], the main effect being a screening of low transverse momentum (p t ) particle production (below a "saturation scale").Saturation effects are expected to be even stronger in nucleus-nucleus collisions [20,21], simply because parton ladders emitted from different nucleons may fuse.At high energies, the diagrams for each scattering actually look more like the one shown in Fig. 3 (I do not consider, for simplicity, time-like parton emis- sions, but in the real EPOS4 simulations, they are of course taken care of).At least for scatterings carrying a large value of x + x − , one expects "nonlinear effects", which means that two ladders which evolve first independently and in parallel, finally fuse.And only after that is the (linear) DGLAP evolution realized.
As mentioned above, such non-linear effects lead to strong destructive interference, which may be summarized in terms of a saturation scale [20,21].This is the motivation to treat these "saturation phenomena" not explicitly, but by introducing a saturation scale as the lower limit of the virtualities for the DGLAP evolutions, as sketched in Fig. 4.So the diagrams inside the red ellipses are replaced by two scales Q 2 sat,proj and Q 2 sat,targ , and in pp scattering the two are equal.So the final version of the "rigorous parallel scattering scenario" in EPOS4 is sketched in Fig. 5.One still has DGLAP evolution, for each of the scatterings, but one introduces sat- . Rigorous parallel scattering scenario, for n = 3 parallel scatterings, including non-linear effects via saturation scales.The red symbols should remind one that the parts of the diagram representing nonlinear effects are replaced by simply using saturation scales.uration scales.Most importantly, as discussed in great detail in Sec.VII, these scales are not constants, they depend on the number of scatterings, and they depend as well on x + and x − .A smart choice of these dependencies allows finally to recover factorization and binary scaling.One understands that there is a (so far unknown) very strong relation between factorization, energy conservation (or better energy sharing), parallel scattering, and saturation, see Fig. 6.Let me summarize the reasoning for this statement: • at high energies, multiple scatterings must happen in parallel, and there is nothing like a sequence or an ordering of elementary collisions; • ignoring energy sharing (as in the GR approach), factorization and binary scaling are obtained (see Sec. IV); • implementing energy sharing in the sense of a "rigorous parallel scattering scenario", not only do the technical difficulties increase enormously, but there are conceptual problems: it spoils factorization and binary scaling (see Sec. VI); • the only way out (it seems) is that one introduces saturation scales in a particular way, which recovers factorization and binary scaling (see Sec. VII).
Having solved the "factorization and binary scaling problem", one may consider the "low p t domain"-like studying high multiplicity (collective) phenomena (see Sec. X)-within a framework that has (finally) been proven to be compatible with the factorization approach, allowing one to do "high p t physics" as well (as the generators that are based on factorization).In this sense EPOS4 is meant to be a "general purpose event generator".
All this discussion about saturation and factorization is of fundamental importance since both are considered to be very important issues, but usually, they are discussed independently.But the message of this paper is that they are connected, they affect each other (see Sec. VII), and they are just two aspects in a common approach.
Since saturation is so important in this approach, what about the other models?The above-mentioned "standard QCD generators" (based on factorization) do not explicitly deal with saturation (apart of a constant low virtuality cutoff), but certain features have similar effects.Let me consider AA scattering in the Pythia/Angantyr model [33,34].As in EPOS4, there is first the "basic AA model" for t = 0, and in a second step there are string interactions, happening later.Concerning the basics AA model, the total S-matrix is given as product of sub-S-matrices.But in contrast to EPOS4, there are no energy-momentum arguments, and therefore no energy sharing.But one needs to introduce some kind of "sequence", i.e., one loops over all the N N scatterings, and treats the interaction in two different ways: if a nucleon is already wounded (having had already a scattering before) the current scattering is realized as diffractive scattering, called secondary scattering; otherwise a "normal" scattering according to Pythia happens, called primary scattering.The latter is in general a multiple parton scattering process.Here, one needs to introduce some ordering.The first sub-scattering is "normal", whereas subsequent ones are not connected to projectile/target remnants as the first one, but they are connected to the parton of the previous sub-scatterings.
So EPOS4 and Pythia/Angantyr are fundamentally different, but in the latter there are certain features that have similar effects as the saturation scale in EPOS: In Pythia/Angantyr one needs some ordering of N N collisions in AA scattering, which is needed to distinguish primary and secondary scatterings.This is necessary to avoid an overproduction of charged particles in AA collisions.In EPOS4, the same effect is obtained by treating all N N scatterings equally, but introducing the dynamical saturation scale.
Concerning multiple parton scatterings in N N collisions, in Pythia/Angantyr the first and the subsequent scatterings are treated differently with regard to the color connections.This is needed to get the experimentally observed increase of the mean transverse momentum with multiplicity.In EPOS4, one treats all subscatterings equally, but one has a saturation scale, which increases with multiplicity (as will be discussed later), and this is the main mechanism that leads to the increase of the mean transverse momentum with multiplicity.
It is of course dangerous to generalize based on few examples, so let me take the following statement as conjecture.I believe, based on the work in this paper, compared to other approaches, that one has two possibilities: • either one considers subsequent sub-scatterings (parton-parton or nucleon-nucleon) as strictly equal, with appropriate energy-sharing, which requires a dynamical saturation scale as crucial element; • or one does not consider saturation (other than a simple cutoff in the parton cascade), but one needs to distinguish between primary and secondary scatterings (the first one and subsequent ones), for both parton-parton and nucleon-nucleon collisions, which requires some ordering.
This paper is meant to be an overview, with a minimum of technical details.The latter can be found in separate publications, such as [35][36][37].
After these introductory remarks, I will • in Secs.IX to XI discuss results affected by saturation.

III. EPOS4 S-MATRIX APPROACH TO REALIZE PARALLEL SCATTERINGS
I first consider pp scattering.An appropriate tool to implement parallel scatterings is provided by S-matrix theory (see [6,[9][10][11]15]), where multiple parallel scatterings can be implemented in a simple and transparent fashion.Factorization and binary scaling are not "assumed", they must come out.The S-matrix is by definition the representation S ij = i| Ŝ |j of the scattering operator Ŝ using some basis of asymptotic states, and the corresponding T-matrix is defined via Particularly important is the diagonal element, T ii , representing elastic scattering, where the asymptotic state |i corresponds to two incoming protons.Assuming purely transverse momentum transfer, one may Fourier transform T ii with respect to the transverse momentum exchange k and in addition divide by 2s to obtain some function T(s, b), with the Mandelstam variable s and the impact parameter b, in the following simply named "T-matrix" T.
The EPOS4 S-matrix approach is based on the hypothesis that the T-matrix T can be written as a sum of products of "elementary" T-matrices Pom , the latter ones representing parton-parton scattering by exchanging a "Pomeron" (without specifying its nature for the moment), as with the "vertices" V (±) representing the connection to the projectile and target remnants.The symbol X stands for all integration variables, to be specified in the following.In Fig. 7, I show a graphical representation of Pom , representing parton-parton scattering, and the magenta dots are the vertices V (±) .The elementary T-matrices are characterized by the light-cone momentum fractions x ± i of the incoming partons, in addition to s and the impact parameter b, so one has The precise content of the Pomerons (boxes) and the functional dependencies on these variables will be discussed later, the general discussion in this section does not depend on these details.The vertices depend on the light-cone momentum fractions of the remnants, x + remn (projectile side) or x − remn (target side), i.e., with a simple functional form (power law) of V.The "δ" in Eq. ( 2) stands for to assure energy-momentum conservation, which will be crucial for the discussions in this paper.The integration dX amounts to integrating over all lightcone momentum fractions.Each term (for n > 1) in the sum of Eq. ( 2) represents multiple scatterings happening in parallel, as it should.
The generalization of the multiple parallel scattering picture towards nucleus-nucleus (AA) collisions (including proton-nucleus) is trivial, one simply writes a product of pp expressions, for colliding two nuclei with mass numbers A and B, with at least one n k > 0. Here, one has one vertex V (m) per remnant, and a sum of products of elementary Tmatrices per nucleon-nucleon pair k.The "δ" in Eq. ( 6) stands for where π(k) = i amounts to summing Pomerons connected to projectile i and τ(k) = j to summing Pomerons connected to target j.This formula does not mean at all a sequence of pp collisions: they are perfectly happening in parallel; the crucial ingredient is the appearance of δ functions.The integration dX here means integration over all light-cone momentum fractions and over all transverse position of the nucleons.The generalization Eq. ( 6) is conceptually trivial, but it should be noted that one has (for big nuclei) 10000000 dimensional nonseparable integrals.So far I have discussed only elastic scattering for pp and AA, the connection with inelastic scattering provides the "optical theorem" (in b representation), which is at high energy given as with cut T = 1 i disc T (cut diagram), with disc T being the s-channel discontinuity T(s + iǫ) − T(s − iǫ).So one needs to compute the "cut" of the complete diagram, cut T; for example, for pp, one needs to evaluate expressions like cut Cutting a multi-Pomeron diagram corresponds to the sum of all possible cuts [38], considering, in particular, all possibilities of cutting or not any of the parallel Pomerons, so one has finally sums of products with some fraction of the Pomerons being cut ("cut Pomerons"), the others not ("uncut Pomerons").I define G to be the cut of a single Pomeron,  10) is finally a product of "G" terms (and in addition the vertex terms V), so G is the fundamental building block of the approach.
Let me consider a simple example of a realization of a Pomeron (the real one is much more complicated), namely a simple parton ladder with two gluon and two quark ladder rungs, see Fig. a vertical red dashed line.For a cut diagram, the Feynman rules are modified in the sense that all elements to the left of the cut are treated normally, for all elements to the right one takes the complex conjugate of the normal result, and all propagators crossing the cut line are replaced by a mass shell condition δ(p 2 − m 2 ).The cut diagram corresponds therefore to an inelastic amplitude squared, with all particles on the cut line being final onshell particles.This is true not only for this simple example but always.So the notion of cut diagrams is very useful, in particular for multiple scattering scenarios.Concerning nuclear scatterings, the total cross section is still given by Eq. ( 9), together with Eqs.(6)(7)(8), which gives for a collision of two nuclei with mass numbers A and B: with at least one n k being nonzero.Here, W AB contains all the vertices and the integration over uncut Pomerons, with the symbol dX explicitly given as representing the nuclear geometry, with the nuclear thickness functions T A (b) given as dz ρ A ( √ b 2 + z 2 ), with ρ A being the nuclear density, and with The impact parameter is defined as , where π(k) and τ(k) refer to the projectile and the target nucleons corresponding to pair k.There is no "δ" term, since here the remnant momentum fractions are no independent variables, they are expressed in terms of the momentum fractions x ± kν as For completeness, and since it is needed in Sec.VI, let me note that W AB in Eq. ( 12) can be written as [36] with some known (simple) function G.
Let me close this section with some technical remarks, concerning the impact parameter dependence and the energy dependence of the T-matrices and of the G functions: • As discussed in [35], the (Mandelstam) t dependence of the original T-matrices is given (in all cases) as factors of the form exp(R 2 t), with parameters R 2 .Considering purely transverse momentum exchange, one has t = −k 2 ⊥ , and the twodimensional Fourier transform with respect to the transverse momentum exchange k ⊥ gives a factor exp −b 2 /(4R 2 ) .
• In this paper, all the "G" and "T" expressions refer to "impact parameter representation"; so all b dependencies are simply Gaussian factors.Here, I do not specify the precise structure of the Pomeron; this is done in very much detail in [35], where is shown that the "real" Pomerons are convolutions of several parts (soft pre-evolution, hard part), but the b dependencies are always Gaussians, giving always a final b dependence of the form of a factor exp −b 2 /(4R 2 ) .So the b dependencies are trivial, and easy to handle.In the following, I will not write the b dependencies explicitly.
• All T-matrices and G functions depend on s, with s referring in all cases to the nucleon-nucleon centerof-mass squared energy, because using the explicit arguments x + and x − , the transverse mass of a Pomeron is x + x − s.In the following, I will not write this s dependency explicitly.

IV. A SIMPLE CASE: FACTORIZATION AND BINARY SCALING IN A SCENARIO WITHOUT ENERGY CONSERVATION
Before further developing the full EPOS4 S-matrix approach, in order to understand the real importance of energy conservation (or energy sharing), I will discuss in this section the S-matrix approach without energy sharing.
I consider the general situation, where the precise structure of Pomeron is not specified.All the diagrams which contribute to cut T (and therefore to the inelastic cross section) in pp represent an infinite series, composed of all possible cut and uncut Pomerons (boxes) as shown in Fig. 10 up to order n = 3.However, here energy sharing will be dropped, which is realized by removing the vertices and the δ term in Eq. ( 2), so one has This simplifies things enormously, since with dX = dX 1 ...dX n , and defining TPom = dX i T (i) Pom (actually not depending on i), one gets which is precisely the expression used in the Gribov-Regge approach.The sub-T-matrix TPom depends only on s and b, as does the full S-matrix T (in both cases I do not write this dependence, for simplicity).Following Eqs. ( 9) and ( 10), including the subsequent discussion, and using Eq. ( 18), one gets for the inelastic cross section with at least one cut.One usually assumes the sub-Tmatrix to be purely imaginary, i.e.TPom = i a 2 , with some real number a, and a factor 1/2 for convenience.Then one gets for the cut Pomeron cut TPom = 2 Im TPom = a.Concerning the uncut Pomerons, one sums up the contributions where the Pomeron is to the left or to the right of the cut, which gives {i TPom } + {i TPom } * = −a.So cut and uncut Pomerons have opposite signs, and one gets where m refers to the number of cut Pomerons.Let me consider inclusive cross sections, like jet cross sections, where m-cut-Pomeron events contribute m times more than single Pomeron events, so one gets where the term in curly brackets represents the sum over all cuts.For a given number n of Pomerons, an elementary calculation allows to compute the sum over all possible cuts, and one finds an amazing result: known as AGK cancellations [11]: • For a given number n > 1 of Pomerons, the sum of all cuts gives zero, i.e.,one gets a complete cancellation.
• Only n = 1 contributes, which corresponds to the case of a single Pomeron.
Therefore, for inclusive cross sections and only for those, only the single Pomeron events contribute, as indicated in Fig. 11.parton evolutions from both sides and a hard elementary parton-parton scattering in the middle, the corresponding inelastic process is shown in Fig. 12(b).So one can write the inclusive pp cross section as a convolution of two parton distribution functions (PDFs) and a (calculable) elementary parton-parton cross section (Born process), which amounts to factorization.So factorization here is the result of a huge amount of cancellations.
For completeness, also for inclusive cross sections in AA scattering, one observes this phenomenon of cancellations, such that finally only a single Pomeron contributes.Colliding two nuclei with mass numbers A and B, the cross section turns out to be AB times the proton-proton cross section (see for example [15]), which is nothing but "binary scaling".
To summarize this section: • in a simplified picture, dropping energy conservation, one gets (easily) factorization for inclusive cross sections in pp scattering • and binary scaling for inclusive cross sections in AA scattering.
But dropping energy-momentum conservation is not really an acceptable solution, in particular since the Monte Carlo procedures eventually need the implementation of energy conservation, so one risks to introduce inconsistencies, in the sense that the theoretical basis and the Monte Carlo realization are not compatible.
In EPOS4, I insist on the compatibility of the theoretical basis and the Monte Carlo realization (this is an evidence, but widely ignored), so one must include energy conservation in the formulas representing the theoretical basis (S-matrix theory).In Secs.V and VI I will discuss why energy conservation spoils factorization and binary scaling, and in Sec.VII I discuss the solution of the problem, which depends strongly on saturation and which leads to "generalized AGK cancellations".

V. HOW ENERGY SHARING DEFORMS POMERON ENERGY DISTRIBUTIONS
In this section, I come back to the full EPOS4 S-matrix formalism, including energy sharing, and try to understand why and how energy sharing affects Pomeron energy distributions, which will be crucial with respect to factorization and binary scaling.
Let me consider a particular Pomeron in AA collisions (including pp as a special case), connected to projectile nucleon i and target nucleon j; see Fig. be other Pomerons, connected to one (or both) of these nucleons.The corresponding Pomeron-nucleon connections are marked as red and blue dots.It is obvious that the additional Pomerons connected to the same nucleons i and j compete with each other: they have to share the initial energy-momentum of the two nucleons.The more Pomerons are connected, the less energy is available for one particular Pomeron.
To quantify this statement, I define the "connection number" with N P being the number of Pomerons connected to i, and with N T being the number of Pomerons connected to j (the variable N conn corresponds to half of the number of red and blue points in Fig. 13).In the following, I will discuss the effect of energy sharing related to the connection number.
As discussed in Sec.III, the fundamental "building block" in EPOS4 is the cut single Pomeron expression G = cut T Pom .As shown in Eq. ( 12), the inelastic cross section for the collision of two nuclei with mass numbers A and B is given in terms of expressions which represent particular configurations {n k } characterized by n k cut Pomerons per nucleon-nucleon pair k (with k between 1 and AB).The indices k, ν refer to the νth Pomeron associated with the pair k.Due to energymomentum conservation, 1 − x + remn i is equal to the sum of all x + kν with k connected to projectile i and 1 − x − remn j is equal to the sum of all x − kν with k connected to target j.The expression Eq. ( 25) represents a multi-dimensional probability distributions for the light-cone momentum fractions x ± kν for a given configuration {n k }.Let me consider, for a given configuration {n k }, a particular Pomeron, which means a particular pair index k and a particular value ν, with the associated momentum fractions x ± kν .Let i and j be the projectile and target the Pomeron is connected to (see Fig. 13).
In the simplest case, one has n k = 1 (only one Pomeron associated to pair k), and one has no other Pomeron connected to i and j, so one has N conn = 1, the case of an isolated Pomeron.Using the known form of W AB given in Eq. ( 16), one can see that the integration of Eq. ( 25) over all variables other than x ± kν gives up to a constant the expression This is the probability distribution of x ± kν for the case N conn = 1, so I name it f (1) (x + , x − ), using simply x ± instead of x ± kν .Using the energy-moment conservation relations, one gets for the N conn = 1 probability distribution with W ′ AB = c W AB with a normalization constant c.Since there may be more than one case with N conn = 1, one averages over them, and one averages over all configurations {n k }.But this does not change anything since they all have the same form, as given in Eq. ( 27).
This two-dimensional distribution Eq. ( 27) allows one to compute the distribution being the probability distribution with respect to the "Pomeron energy fraction", with M Pom being the transverse mass of the Pomeron, with y PE = 0.5 ln x + x − , and with J being the corresponding Jacobian determinant.Having very narrow y PE distributions, one may use and one gets Let me now consider a more complicated situation, corresponding to N conn > 1.The general formula for the probability distribution f (N conn ) (x + , x − ) is given as where δ b a is the Kronecker delta, and where P(K) is the expression Eq. ( 25) for particular multi-Pomeron configurations K = {n k }, {x ± kν } , with given energymomentum sharing.The symbol ∑ {n k } =0 means summing over all possible choices of n 1 , n 2 , ..., n AB , excluding the case where all n k are zero.The symbol dX is explicitly given as d 2 b dT AB dX AB , with dT AB and dX AB defined in Eqs. ( 13) and ( 14).I only consider Pomerons (k ′ , ν ′ ) with connection number N conn (k ′ , ν ′ ) equal to N conn .In principle one does the same procedure as for the case N conn = 1, namely for a given configuration {n k }, one chooses a particular pair index k ′ and a particular value ν ′ , with the associated momentum fractions x ± k ′ ν ′ , which are replaced by x ± after integrating over these variables, because of the delta functions.Then one integrates Eq. ( 31) over all variables other than x ± k ′ ν ′ , which is always possibe.But in this case, these integration variables and the "chosen variables" x ± k ′ ν ′ (now x ± ) can no longer be separated, and one needs to do the integration numerically (via Monte Carlo, in practice).Nevertheless, it is well defined, and one gets the x ± distribution f (N conn ) (x + , x − ).As for N conn = 1, one integrates over y PE , to get f (N conn ) (x PE ).In practice, one defines event classes (according to multiplicity or impact parameter) and computes the average N conn values as well as the average x ± distributions per class, to finally get f ( N conn ) (x PE ).Then one takes the obtained distributions as the basis to compute f (N conn ) (x PE ) for arbitrary values of N conn via interpolation.
The distribution f (N conn ) (x PE ) or the two-dimensional distribution f (N conn ) (x + , x − ) are kind of "master" distributions for all kinds of "inclusive distributions", for example the inclusive p t distribution of partons, or of hadrons if one adds a fragmentation function.For computing just inclusive spectra, the knowledge of f is enough, whereas otherwise the full calculations are needed, using Monte Carlo methods based on Markov chains.In addition to f (N conn ) (x PE ) I also check the inclusive distribution of y PE , which is narrow and strongly peaked at zero (note that the x ± refer to Pomeron momenta, not to those of outgoing partons).This is why I concentrate in the following on f (N conn ) (x PE ).
Although f (N conn ) (x PE ) for N conn > 1 cannot be calculated analytically, one has some idea of how it should look like compared to f (1) (x PE ): the integrand of Eq. ( 31) contains in addition to other G terms and most importantly W AB ({x + remn i }, {x − remn j }), contains factors of the form (x ± remn i ) α with α > 0, with arguments where ∑ ′′ sums over all indices different from k ′ , ν ′ , being connected to i or j.Due to the additional term kν .This is also what one expects without any calculation: energy sharing involving more than one Pomerons leads to a reduction of the energy of the Pomerons, compared to the case of an isolated Pomeron.
But let me be quantitative, and discuss the real calculations.In Fig. 14, I plot f (1) (x PE ) and f (N conn ) (x PE ) for the centrality class 0-5% in PbPb collisions at 5.02 TeV with an average value of N conn of roughly 7.7, obtained after a full EPOS4 simulation.One observes (as expected) for N conn > 1 a "deformation" of the the x PE distributions compared to f (1) (x PE ), due to energy-momentum conservation.Therefore I define the ratio called the "deformation function".In Fig. 15, I show the deformation function for the centrality class 0-5% in PbPb collisions at 5.02 TeV.The functional form is as one Let me quickly summarize the main results of this section: • Imposing energy sharing (as it should be) has a very important impact on the distribution of Pomeron energies.
• A useful variable to quantify the effect of energy sharing is the connection number N conn , counting the number of "other Pomerons" connected to the same remnants as a "given Pomeron".N conn = 1 represents an isolated Pomeron.
• I define a variable x PE = x + x − , representing the squared energy per Pomeron, and the corresponding probability distribution.
• The probability distribution depends strongly on N conn , so I use the notation f (N conn ) (x PE ).The results for N conn > 1 show a suppression at large x PE , as a consequence of energy sharing.This must be so, it is unavoidable, a fundamental feature.
• I therefore define a "deformation function" R deform as the ratio of f (N conn ) (x PE ) over f (1) (x PE ), which drops below unity for large x PE .

VI. HOW DEFORMED POMERON ENERGY DISTRIBUTIONS SPOIL FACTORIZATION AND BINARY SCALING IN CASE OF A "NAIVE" POMERON DEFINITION
In this section, the aim is to understand why and how energy sharing ruins factorization and binary scaling.I showed in the last section that energy sharing leads unavoidably to a "deformation" of the Pomeron energy distribution f (N conn ) (x PE ) compared to the reference f (1) (x PE ), with N conn being the connection number, counting the number of other Pomerons connected to the same remnants as a given Pomeron, which leads to the definition of a "deformation function" R deform as the ratio of f (N conn ) (x PE ) over f (1) (x PE ).In the following, I will show how this deformation spoils factorization.
Actually "the problem" related to factorization and binary scaling depends very much on the precise definition of G in terms of QCD.Sofar a Pomeron is a black box with all the QCD details hidden inside, but now one needs to be more specific.To do so, one introduces in [35] (with many details and all necessary formulas) a "QCD expression" called G QCD representing a QCD calculation of parton-parton scattering.Every "G function" which one uses, including G QCD , is meant to be the cut of the Fourier transform of the T-matrix, divided by 2s.The term G QCD is a sum of several contributions, the most important one being the "sea-sea" contribution G sea−sea QCD ; see Fig. 16.For a precise definition see [35].The vertices F i the projectile and target nucleons.In addition, one has three blocks: the two soft blocks and in between a parton ladder, the latter being a DGLAP parton evolution from both sides, with a pQCD Born process in the middle.I define parton evolution functions E QCD obeying the usual DGLAP equations, but in this case the evolution starts from a parton, not from a nucleon, since a Pomeron corresponds to parton-parton scattering.I compute and tabulate E QCD , and then the convolution with an elementary QCD cross section "Born" and a soft preevolution E soft .In addition to "sea-sea", one has more contributions, named "val-val", "sea-val", "valsea", "soft", and "psoft", as discussed in great detail in [35].Like G, the QCD expression G QCD depends as well on x + and on x − , and in addition there is the crucial parameter Q 2 0 , which is the low virtuality cutoff in the DGLAP evolution, so I use the notation G QCD (Q 2 0 , x + , x − ).Whereas the cutoff is usually a constant of the order of 1 GeV 2 , I consider it as a variable that may take any value, and I compute and tabulate G QCD (Q 2 , x + , x − ) for large ranges of discretized values of all arguments, such that G QCD can be computed via interpolation for any choice of arguments.After these preparations, the functional form of G QCD (Q 2 , x + , x − ) is known.Actually, G and G QCD also depend on s and b, not written explicitly as discussed earlier, so one should always consider "for given s and b".
What is the relation between G (the Pomeron, the main building block of the multiple scattering theory) and G QCD (which contains all the QCD part)?A first attempt might be (and this is what was actually used in [15]) to consider the two to be equal, i.e., with a constant Q 2 0 .Then one gets for the Pomeron energy distribution for an isolated Pomeron corresponding to N conn = 1 [see Eq. (30)] with the Pomeron energy variable x PE = x + x − , and using which means that the x PE distributions will get more and more "deformed", in particular suppressed at large x PE .This is a general feature and is unavoidable, a direct consequence of energy sharing.What does this mean concerning transverse momentum (p t ) distributions of the outgoing particles from the Born process?Here one needs to consider the internal structure of G QCD , first of all G sea−sea QCD (actually similar arguments hold for the other contributions).The important element is the parton ladder, see Fig. [35], given as a convolution (for the formulas, see [35]).The p t of the outgoing partons is related to the factorization scale µ 2 F (one uses µ 2 F = p 2 t ), which is the virtuality of the partons entering the Born process.Large values of p t require large µ 2 F and large squared energy ŝ of the Born process, and this requires a large Pomeron squared energy, and therefore a large value of x PE .The essential points are • Large values of p t of the outgoing partons are strongly correlated with large values of x PE .
• A suppression of large x PE values in f (N conn >1) (x PE ) compared to f (1) (x PE ) will therefore lead to a suppression of large p t values in the case of N conn > 1 compared to N conn = 1.
In Fig. 17, I sketch this situation of a suppression of parton yields at high p t with increasing N conn .Let me first discuss the consequences for pp scattering.In order to get factorization, as discussed in Sec.IV, one would need something like AGK cancellations, such that the full multiple scattering scenario is identical to the single Pomeron (N conn = 1) case, which means eventually for the minimum bias (MB) inclusive particle production result.The latter may be written as a superposition of the different contribution for given values of N conn [the latter being (for pp) identical to the number of cut Pomerons] with the corresponding weights w (N conn ) : where the contributions dn (Nconn) dp t show with increasing N conn more and more suppression at large p t (as indicated in Fig. 17).The average Pomeron number at say 7 TeV is around 2, so one has definitely an important contribution from terms with N conn > 1.This means that also the MB result will be reduced at high p t compared to N conn = 1, and this means one cannot fulfill Eq. ( 39), so factorization is not achieved.
The discussion for the scattering of two nuclei with mass numbers A and B is similar to the pp case.To assure binary scaling, one expects for the minimum bias (MB) inclusive particle yields.For the contributions for different values of N conn , one has a picture similar to that shown in Fig. 17, and also here one concludes that the MB inclusive yield will be reduced at high p t compared to N conn = 1, and therefore one cannot fulfill Eq. ( 41), and therefore binary scaling is violated.Let me summarize this section: • I consider here the case where the cut Pomeron G is identical to G QCD , the latter representing a pQCD result for parton-parton scattering.
• Considering the internal structure of G QCD , one concludes that there is a strong correlation between the Pomeron energy variable x PE and p t of the outgoing partons (large p t corresponds to large x PE ).
• Therefore the suppression of large x PE with increasing N conn amounts to a suppression of large p t , and one can conclude a suppression of yields at large p t for minimum bias results compared to N conn = 1.
• This means one cannot obey Eqs. ( 39) and ( 41), which are necessary conditions for factorization (pp) and binary scaling (AA).

VII. HOW SATURATION ALLOWS ONE TO RECOVER FACTORIZATION AND BINARY SCALING (GENERALIZED AGK CANCELLATIONS)
In the following, I discuss the "key issue" of the EPOS4 approach, namely the appropriate definition of G(x + , x − ), the cut Pomeron, represented so far as a "cut box" as shown in Fig. 18, and used earlier (see for ex- ample Figs. 8 and 10), to develop the multiple scattering scheme.The latter is actually completely independent from the precise definition of G, which is very useful, so one can investigate different options concerning the internal structure of G.
I showed in the previous section that the "naive" assumption G = G QCD (42) (which was also adopted in [15] and [39]) completely spoils factorization and binary scaling.And from the discussion in the previous section, it is known that this is a fundamental, unavoidable problem, and not just a wrong parameter choice.So the assumption Eq. ( 42) seems to be simply wrong.
There is another serious problem with Eq. ( 42): as discussed somewhat in the previous section (and in detail in [35]), the essential part of G QCD is a cut parton ladder, based on DGLAP parton evolutions.But as already discussed in the introduction, this is certainly not the full story: with increasing energy, partons with very small momentum fractions x ≪ 1 become increasingly important, since the parton density becomes large, and therefore the linear DGLAP evolution scheme is no longer valid and nonlinear evolution takes over, considering explicitly gluon-gluon fusion.These phenomena are known as "small x physics" or "saturation" [19][20][21][22][23][24][25][26][27][28][29][30][31][32].
At least for scatterings carrying a large value of x + x − , one expects "nonlinear effects", which means that two ladders which evolve first independently and in parallel, finally fuse.And only after that is the (linear) DGLAP evolution realized.Such nonlinear effects lead to strong destructive interference at low transverse momentum (p t ), which may be summarized in terms of a saturation scale [20,21].This suggests treating these "saturation phenomena" not explicitly, but by introducing a saturation scale as the lower limit of the virtualities for the DGLAP evolutions, as sketched in Fig. 19.

nonlinear effects nonlinear effects
Figure 19.Nonlinear effects (inside the red ellipses) also referred to as saturation effects are "summarized" in the form of saturation scales, which replace these non-linear parts.

So one has two problems:
• a wrong identity G = G QCD , • a missing treatment of saturation.
But fortunately, the two problems are connected, and there is an amazingly simple solution that solves both problems.Instead of the "naive" assumption G = G QCD , one postulates: with Here, R (N conn ) deform (x PE ) is the deformation function discussed in Sec.V, and n is a constant, not depending on x PE .The independence of G on N conn is absolutely crucial (as I will show later) and to achieve this, one first parametrizes G based on the (very reasonable) assumption that G has a "Regge-pole structure" as G ∝ α x PE β , where the s and b dependences of α and β are parametrized with the parameters being fixed by comparing simulation results to elementary experimental data, and then use Eq. ( 43) to determine Q 2 sat .In this way, Q 2 sat depends on N conn and on x ± : which means that this Q 2 sat , being the low virtuality cutoff for the DGLAP evolutions in G QCD , is not a constant, but its value depends on the environment in terms of N conn and on the energy of the Pomeron.I will refer to this as "dynamical saturation scale".
But why does Eq.( 43) work?One gets for the Pomeron energy distribution for an isolated Pomeron, corresponding to N conn = 1 [see Eq. ( 30), using where Eq. ( 43) with N conn = 1 was used to replace G.In the case of N conn > 1, one has [see Eq. ( 33)] ∝ R Here I will use again Eq. ( 43) to replace G, but with deform expressions cancel, and one gets The crucial point here is the fact that thanks to Eq. ( 43) and since G does not depend on N conn , the R ) This equation is very interesting, it means that the N conn dependence of x PE distributions is guided by the saturation scale, and nothing else.This is the only difference between the numerator and the denominator.The Eq. ( 51) also means that the partonic structure is given by G QCD , and therefore also the p t distribution of the outgoing partons is encoded in the single Pomeron expression G QCD , for any N conn .Only the saturation scales Q 2 sat depend on N conn , and these saturation scales suppress small p t particle production, but do not affect high p t results, as sketched in Fig. 20.
What does this mean concerning factorization?The minimum bias (MB) inclusive parton yield may be written as a superposition of the different contribution for given values of N conn (for pp identical to the number of cut Pomerons) with weights w (N conn ) : At large p t , all contributions are equal, as indicated in Fig. 20, so one can replace dn  so only a single Pomeron contributes.This will allow one to define parton distribution functions f PDF and to compute cross sections as convolutions f PDF ⊗ Born ⊗ f PDF .
The discussion for the scattering of two nuclei with mass numbers A and B is similar to the pp case, the Pomeron connections do not affect high p t , as shown in Fig. 20.The only difference compared to pp scattering is the fact that one has to sum over all possible nucleonnucleon pairs, which gives for the minimum bias (MB) inclusive particle yields, which amounts to binary scaling.
Equations ( 53) and ( 54) state the following: • The inclusive pp cross section is equal to the one of a single Pomeron contribution, and the inclusive cross section of the scattering of two nuclei of mass numbers A and B is equal to AB times the single Pomeron contribution, leading to factorization and binary scaling.
• I refer to this as the "generalized AGK theorem", valid at high p t , in a scenario with energy sharing.
One recalls that the classical AGK cancellations [11] are based on a scenario without energy sharing, as discussed in Sec.IV.
Let me summarize this section: • One tries to find the relation between G (the multiple scattering building block) and G QCD which represents a QCD result concerning single parton parton scattering.
• Two problems are identified: (1) the naive expectation G = G QCD , having been used so far, does not work, and (2) an appropriate treatment of saturation is missing.
• Both problems are solved by postulating G(x ± ) ∝ G QCD (Q 2 sat , x ± )/R deform , which means that a saturation scale, depending on the Pomeron connection number N conn and on x ± , replaces the virtuality cutoff Q 2 0 usually used in DGLAP evolutions.In this way one incorporates saturation.
• A direct consequence of the above postulate is the fact that the Pomeron energy distribution f (N conn ) (x PE ) is for any value of N conn given in terms of a single Pomeron expression G QCD , with an only implicit N conn dependence via Q 2 sat .• As a consequence, N conn affects low p t (suppression) but not high p t , and one recovers factorization and binary scaling (generalized AGK theorem).
As a final remark: within a rigorous parallel scattering scenario (which seems mandatory), and respecting energy conservation (which seems mandatory as well), the only way to not get in contradiction with factorization and binary scaling seems to be the consideration of saturation via G = k × G QCD (Q 2 sat ) with k being inversely proportional to the deformation function.In this sense, parallel scattering, energy conservation, saturation, and factorization are deeply connected.

VIII. REMARKS CONCERNING DEFORMATION FUNCTIONS AND SATURATION SCALES FOR GIVEN EVENT CLASSES
Let me come back to the deformation function R (N conn ) deform (x PE ), which plays a fundamental role in the new approach.As explained earlier, this function can be computed based on Monte Carlo simulations.But to do so, one first needs to define the Pomerons.This is done using a parametrization of G in "Regge pole form" αs β , and based on this, one computes the deformation functions.Then one uses Eq. ( 43), to do full simulations, and compare with data.If needed, the initial parametrization of G, and as a consequence also the deformation functions, are changed, and one repeats the procedure.In practice, I have found a very simple functional form for the deformation functions which accommodates all systems, centrality classes, and energies [36].I determine and tabulate the parameters, and then use parametrized deformation functions.In Fig. 21, I show as an example the deformation function for the centrality class 0-5% in PbPb collisions at 5.02 TeV, with an average value of N conn of roughly 7.7.I plot the parametrized function (black dots) and the "computed distribution" which comes out after a simulation (red curve).The two curves agree.
One recalls that in the above iterative procedure, Eq. ( 43) is used based on a G already known and parametrized in "Regge pole form" αs β .As a historical side remark, including saturation effects via a Regge pole form αs β was introduced first in EPOS1 [40], where the term "parton ladder splitting" was used rather than "saturation", but it refers to the same phenomenon.A dynamical saturation scale, assuming some functional form of Q 2 sat (N conn , x PE ), was introduced first in EPOS3 [41] for proton-lead scattering, where also for the first time the expression "saturation" was used.Also in [41], real simulation results show (Fig. 3) the suppression of parton yields at high p t in case of the "naive" assumption G = G QCD , and it is shown (Fig. 4) that the suppression can be avoided by introducing a saturation scale.Finally in [42] it was proposed to use a parametrized G as in [40] and use it to determine Q 2 sat (N conn , x PE ).However, at the time, the role of the deformation was not yet understood, but this is crucial to ensure the correct asymptotic behavior at large x PE and large p t .
As was done for the case of central PbPb in Fig. 21, one computes the deformation functions always for event classes and associates the obtained function to the mean N conn for the corresponding class.There are several ways to define event classes; one possibility is to do it based on the number of cut Pomerons (or simply the Pomeron number) N Pom , which is related to the multiplicity dn/dη(0).I consider simulations for pp at 7 TeV and PbPb at 2.76 TeV (because I will later come back to these two systems to compare simulation results with data).In Fig. 22, I show the multiplicity dn/dη(0) as a function of the Pomeron number N Pom , for pp at 7 TeV (red line) and PbPb at 2.76 TeV (blue line), together with the dotted line representing the function 2.9(N Pom ) 0.9 , which provides a simple conversion formula between these two quantities.This might be useful when I analyze later observables as a function of dn/dη(0) .One gets a continuous curve when going from pp to PbPb.In Fig. 23, I plot the saturation scale Q 2 sat as a function of x PE , for several N Pom event classes.The most striking result is the fact that in pp the Q 2 sat values change very  strongly with N Pom , whereas for PbPb the variation is quite moderate, and towards central collisions Q 2 sat even "saturates" (no variation anymore).This discussion will be important to understand the results in Sec.X.
As a first result of EPOS4 simulations, I am going to show that binary scaling really works in practice.In Fig. 24, I show the inclusive p t distribution of partons for a full simulation (simu) divided by N coll and the "reference curve" (theo), which is the corresponding distribution for a single Pomeron, calculated analytically.I show results for minimum bias PbPb collisions at 5.02 TeV (red curve) as well as results for different centrality classes.One can see that the ratio is close to one for large values of p t , whereas low p t values are suppressed.Also in pp, the full simulation over the "reference curve" (single Pomeron) is close to unity at large p t and comparing pp and AA, one gets R AA ≈ 1.An iterative procedure is employed that relies very much on experimental data: One starts with a parametrization of G in "Regge pole form", already constrained by basic experimental data like energy dependence of cross sections.Then one computes the deformation functions, which allows finally to determine saturation scales, which allows in this way including saturation effects.Based on G QCD , with a (now) known saturation scale, one can generate partons, and then make very detailed comparisons with all kinds of data, and if needed redo the procedure with an improved parametrization of G.So to some extent one has a data-driven method to obtain saturation scales, based on a fully self-consistent pQCD based multiple scattering framework, which is complementary to efforts of computing saturation scales.

IX. EPOS4 FACTORIZATION MODE (SINGLE POMERON) AND EPOS4 PDFS
Since in the case of inclusive spectra at large p t everything can be reduced to the single Pomeron case (generalized AGK cancellations, see Sec.VII), one may use "a shortcut" and compute inclusive particle production simply by using a single Pomeron, without any need to use complicated Monte Carlo procedures.This is referred to as "EPOS4 factorization mode".In this case, one simply needs to evaluate the cut single Pomeron, corresponding to G QCD , which is composed of several contributions; see [35].The most important one is G sea−sea QCD ; see Fig. der with parton evolutions (E QCD ) from both sides and an elementary QCD Born process in the middle.In addition, the QCD parton evolution is preceded by a soft evolution (E soft ).The vertices F sea 1 and F sea2 couple the parton ladder to the projectile and target nucleons.The complete expression is a convolution of several elements which in addition needs to be convoluted with the vertices V as V ⊗ ... ⊗ V.This expression may be regrouped in several ways.One possibility is to convolute first the vertices, the soft evolution, and the QCD evolution on the projectile side, representing the parton distribution function (PDF) of the projectile, and correspondingly on the target side.The two PDFs represent actually the upper and lower part of the graph in Fig. 25, plus a vertex V, but excluding the Born process.So far I only consider the so-called "sea-sea" contribution G sea−sea QCD , with a sea quark (after a soft evolution) being the first parton entering the partonic cascade on both sides.But as shown in [35] there is, in addition, a "valval" contribution, where valence quarks enter the partonic cascade, and correspondingly "val-sea" and "seaval" contributions.Since the parton distribution function is just half of the Pomeron diagram, there are two contributions, the "sea" and a "val" one.For a precise definition of the PDFs, see [35].
One computes (and tabulates) the PDFs f k PDF (x, µ 2 F ), with x being the light-cone momentum fraction of the parton of flavor k entering the Born process, and µ 2 F being the factorization scale.After this preparation, one may express the di-jet cross section (where di-jet simply refers to the two outgoing partons of the Born process) in terms of the PDFs, as with p 1/2 and p 3/4 being the four-momenta of the incoming and outgoing partons, and M kl→mn being the corresponding matrix element.In order to get the complete expression corresponding to Fig. 25, one needs to integrate over the differential cross section Eq. ( 55), whereas to obtain the inclusive jet (=parton) cross section one needs to integrate . In any case, thanks to the four-dimensional δ function, the remaining numerical integration can be done, as discussed in detail in [35].
At least the quark parton distribution functions can be tested and compared with experimental data from deep inelastic electron-proton scattering.The structure function F 2 is given as 2pq), with p being the momentum of the proton, q the momentum of the exchanged photon, and Q 2 = −q 2 .In Fig. 26, I plot F 2 as a function of x for different values of Q 2 .The red curve refers to EPOS PDFs, the green one to CTEQ PDFs [43], and the black points are data from ZEUS [44] and H1 [45][46][47].The two PDFs give very similar results, and both are close to the experimental data.
x 5000 Having checked the EPOS PDFs, I will use these functions to compute the jet (parton) cross section, using Eq. ( 55), integrating out the momentum of the second parton and the azimuthal angle of the first parton, for pp at 13 TeV.I define the parton yield dn/dp t dy as the cross section, divided by the inelastic pp cross section, showing the result in Fig. 27.I show results based on EPOS PDFs (red full line), CTEQ PDFs [43] (green dashed line), the full EPOS simulation (blue circles), and experimental data from ATLAS [48] (black triangles).At large values of p t , all the different distribution agree, whereas at low p t the EPOS Monte Carlo simulation results (using the full multiple scattering scenario) are significantly below the PDF results, as expected due to screening effects.

X. FULL EPOS4 (CORE+CORONA, HYDRO, MICROCANONICAL DECAY): CHECKING MULTIPLICITY DEPENDENCIES
The "factorization mode" as discussed in the last section is very useful to investigate inclusive cross sections at high p t .But this represents only a very small fraction of all possible applications, and there are very interesting cases outside the applicability of that approach.A prominent example, one of the highlights of the past decade, concerns "collective phenomena in small systems", following many discoveries showing that highmultiplicity pp events show very similar "collective" features as earlier observed in heavy ion collisions [5].
High multiplicity means automatically "multiple parton scattering", and as discussed earlier, this means that one has to employ the full parallel scattering machinery developed earlier, based on S-matrix theory.But this is not the full story.The S-matrix part concerns "primary scatterings", happening instantaneously at t = 0.As a result, in the case of a large number of Pomerons, one has correspondingly a large number of strings, which may overlap and "fuse".In the EPOS framework, a core-corona procedure [49][50][51] is employed, where the strings at a given proper time τ 0 are first cut into "string segments", which are separated into "core" and "corona" segments, depending on the energy loss of each segment when traversing the "matter" composed of all the other segments.Corona segments (per definition) can escape, whereas core segments lose all their energy and constitute what is called the "core", which acts as an initial condition for a hydrodynamic evolution [51].The evolution of the core ends whenever the energy density falls below some critical value ǫ FO , which marks the point where the fluid "decays" into hadrons.It is not a switch from fluid to particles; it is a sudden decay, called "hadronization".
In EPOS4, as discussed in detail in [37], a new procedure was developed concerning the energy-momentum flow through the "freeze-out (FO) hypersurface" defined by ǫ FO , which allows defining an effective invariant mass, which decays according to microcanonical phase space into hadrons, which are then Lorentz boosted according to the flow velocities computed at the FO hypersurface.Also new and very efficient methods for the microcanonical procedure [37] were developed.Also in the full scheme, including primary and secondary interactions, energy-momentum and flavors are conserved.All the technical details about the new hadronization procedures can be found in [37]; the aim of this paper is to present an overview and some important results.
As an alternative, in order to better understand the different components, I also consider • the "core+corona" ("co-co") contribution, i.e. primary interactions + secondary interactions but without hadronic afterburner; • the "core" contribution, i.e. primary interactions + secondary interactions but without hadronic afterburner, only considering core particles; • the "corona" contribution, i.e. primary interactions + secondary interactions but without hadronic afterburner, only considering corona particles.
One needs to exclude in these cases the hadronic afterburner, because the latter affects both core and corona particles, so in the full approach, the core and corona contributions are no longer visible.
In the following, I will present particle ratios, always relative to pion yields, as well as mean p t results, for the different contributions ("core", "corona" etc), in pp and PbPb collisions at LHC energies.In all cases, the results depend strongly on the relative weight of core to corona.It is clear that for low multiplicity pp scattering corona will dominate, whereas, for central PbPb collisions, the core will dominate.To be more quantitative, I compute the "core fraction", defined as the ratio of core to core+corona for pion production (with pions being the most frequent particle species).In Fig. 28, I show results for pp (thin lines) and PbPb (thick lines), and one sees an almost continuous curve, going from zero (for low multiplicity pp) up to unity (for central yields over π yields versus multiplicity.I show results for pp at 7 TeV (thin lines) and PbPb at 2.76 TeV (thick lines), compared to ALICE data [54,55].The different line styles refer to different contributions: the yellow dashed line refers to "core+corona" ("co-co"), i.e. primary interactions + hydro but without hadronic after-burner, the blue dotted line refers to the "corona" and the green dashed-dotted line refers to the "core" part.The red line is the "full" contribution, i.e. core + corona + hadronic afterburner.One sees an almost flat line for the corona contribution, similar for pp and PbPb, which is understandable, since "corona" means particle production from string fragmentation, which does not depend on the system.One observes also a flat curve for the "core" part at high multiplicity, which is again expected since the core hadronization is determined by the freeze-out energy density, which is system independent.However, when the system gets very small, one gets a reduction of heavy particle production due to the microcanonical procedure (with its energy and flavor conservation constraints), whereas a grand canonical treatment would give a flat curve down to small multiplicities.It is remarkable that the "core" curve is far above the "corona" one, which simply reflects the fact that Ω production is much more suppressed in string decay, compared to statistical ("thermal") production.This explains why the "core+corona" contribution increases by one order of magnitude from low to high multiplicity, because simply the relative weight of the core fraction increases from zero to unity.The effect from hadronic rescattering (difference between "full" and "co-co") is relatively small, some suppression due to baryon-antibaryon annihilation can be seen.Whereas the Ω over π ratios are essentially smooth curves, from pp to PbPb, the situation changes completely when looking at the average transverse momentum p t versus multiplicity, as shown in Fig. 29(lower panel), where I also show results for pp (thin curves) and PbPb (thick curves), for the different contributions.Here one sees (for all curves) a significant discontinuity when going from pp to PbPb.The "corona" contributions are not flat (as the ratios), but they increase with multiplicity, in the case of pp being even more pronounced than for PbPb.This is a "saturation effect": the saturation scale increases with multiplicity, which means with increasing multiplicity the events get harder, producing higher p t .The situation is different for PbPb, where an increase of multiplicity is mainly due to an increase of the number of active nucleons, with a more modest increase of the saturation scale with multiplicity.Also, the "core" curves increase strongly with multiplicity, and here as well more pronounced in the case of pp, due to the fact that one gets for high-multiplicity pp high energy densities within a small volume, leading to strong radial flow.
Again, the core+corona contribution is understood based on the continuous increase of the core fraction from low to high multiplicity.
In Figs. 30 and 31, I show the multiplicity dependencies of ratios and mean p t for different hadrons, which are qualitatively similar to the Ω results, just the difference between the corona and the core curves are smaller.The data are from ALICE [54][55][56][57][58].
It is very useful (and necessary) to consider at the same time the multiplicity dependence of particle ratios and of mean p t results, since their behavior is completely different (the former is continuous, the latter jumps).De- spite these even qualitative differences between the two observables, the physics issues behind these results is the same, namely saturation, core-corona effects which mix flow (being very strong) and non-flow, and microcanonical hadronization of the core.Another very important and useful variable is the multiplicity dependence of D meson production, where D stands for the sum of D 0 , D + , and D * + .This is much more than just "another particle", since the D meson contains a charm quark, the latter one being created ex-clusively in the parton ladder and not during fragmentation or in the plasma.In Fig. 32, I [59].It is interesting to see in which way the simulations and the data deviate from the reference curve, which is the dashed black line representing identical multiplicity dependence for D mesons and charged particles.Considering the EPOS results without hydro (green lines), for low p t (1-2GeV/c) the curve is slightly above the reference, but with increasing p t the green curves get steeper, which is due to the fact that with increasing multiplicity the saturation scale increases, and the events get harder, producing more easily both high p t and charmed particles.Considering EPOS with hydro (red curves), the increase compared to the green curves is much stronger, which is due to the fact that "turning on hydro" will reduce the multiplicity (the available energy is partly transformed into flow).The red curves are close to the experimental data, both showing a much stronger increase compared to the reference curve, with the effect getting bigger with increasing p t .So one may conclude this paragraph: to get these final results (the strong increase), two phenomena are crucial, namely, saturation which makes high multiplicity events harder, and the "hydro effect" which reduces multiplicity and "compresses" the multiplicity axis.
Concerning earlier EPOS versions, there are no "real publications" concerning these multiplicity dependencies, only plots based on on "preliminary versions" shown at conferences or given to experimental colleagues.But none of the preliminary versions were able to fit reasonably well at the same time all the data shown in this section.

XI. CHARMED HADRONS
Having already discussed the multiplicity dependence of charm production in the last section, I will show here some basic charm results (a detailed discussion about charm production can be found in [35]).I consider here just primary interactions, no hydro and no hadronic cascade, so the charm quarks originate from cut Pomerons, more precisely from the parton ladder.Cut parton ladders correspond in general to two chains of partons q − g − ... − g − q identified as kinky strings, with q referring to light flavor quarks, and g to gluons.The Born process or branchings in the space-like or the time-like cascade may lead to Q Q production, where Q refers to "heavy flavor" (HF) quarks, i.e. charm or bottom.In this case, one ends up with parton chains of the type q − g − TeV.The red lines refer to EPOS simulations, the green points to FONLL calculations [60], and the black points to ALICE data [61].In Fig. 34, I plot transverse momentum spectra of Λ c and Ξ c baryons (upper panel) and their ratio with respect to D 0 mesons in pp collisions at 5 TeV.The red lines refer to EPOS simulations, and the black points to ALICE data [62,63].The production of charmed baryons is in principle straightforward, they are also coming from q − g − ... − g − Q and Q − g − ... − g − q strings (with Q being a c quark in this case).The only difference compared to charmed meson production is the fact that here a diquark-antidiquark breakup occurs, which results in an essentially flat baryon / meson ratio, whereas the data show an increase towards small p t .A similar "baryon/meson enhancement" in the region around 2 -6 GeV/c has already been observed in the light flavor sector, where one possible explanation is collective flow.Since charm is produced (as everything else) in an event-by-event manner, one produces for each charm quark the corresponding charm antiquark, and depending on the production details they have characteristic correlations, which are also visible in D meson pair correlations.In Fig. 35, I show two-hadron correlations for D 0 D 0 (upper left plot), D 0 D0 (upper right), D + D − (lower left), and D + D − s (lower right) as a function of the rapidity difference ∆y in pp collisions at 7 TeV, with p t values between 3 and 12 GeV/c and rapidities between 2 and 4. Red lines represent EPOS4 simulations and black dots data from LHCb [64].In Fig. 36, I plot the correlations of these pairs as a function of the azimuthal angle difference ∆φ, again compared to LHCb.
It should be noted that D 0 D 0 represents a c − c correlations, whereas the three other combinations D 0 D0 , D + D − , and D + D − s represent c − c correlations.For the latter ones, the situation is quite simple: the c and the c are always produced as pair from the same process, and therefore one expects them to be close in rapidity, with a preference of ∆φ = 0 (in case of a time-like g → c c), or at ∆φ = π in case of a Born process gg → c c.This is precisely what is seen: The rapidity correlations have maxima at ∆y = 0 and then drop quickly, the ∆φ correlations have maxima et ∆φ = 0 and ∆φ = π, observed in both EPOS4 simulations and data.
Surprisingly, the D 0 D 0 correlations (corresponding to a c − c pair) look very similar, which suggests that also  c − c pairs originate from the same process, like a timelike g → gg → c c c c or a Born process gg → gg followed by g → c c, g → c c. Since EPOS4 creates charm always in terms of c − c pairs, it is quite tempting to look into the possibility to produce charmonium.It is easy to implement the idea of the color evaporation model [65][66][67], where charmonium is created with a certain probability in the case of a c − c pairs being in the appropriate mass range.So one considers all c − c pairs from the same Pomeron (fully evaluated, including time-like emissions), and compute the invariant mass M c c. Whenever this mass is less than the sum of two D meson masses and bigger than the J/Ψ mass, the c − c pair is with a certain probability w J/Ψ considered to be a J/Ψ.In Fig. 37, I plot prompt J/Ψ (not coming from beauty decays) from EPOS4 simulations compared to ATLAS data [68].

XII. SUMMARY
I reported on new ideas, implemented in EPOS4, which provide a new understanding of a deep connection between four basic concepts in pp and AA collisions: rigorous parallel scattering, energy conservation, factorization, and saturation.It is mandatory to treat multiple scatterings in parallel, and a "natural" framework is S-matrix theory, with an S-matrix being given as a product of several entities representing individual scatterings referred to as Pomerons, and with energymomentum conservation being implemented in an unbiased fashion via δ functions without imposing any ordering of collisions (this is what is meant by rigorous parallel scattering).The fundamental quantity of the multiple scattering approach is the cut single Pomerons expression G, representing inelastic parton-parton scattering.The fundamental question discussed in this paper is how to relate G to G QCD (Q 2 0 ), where the latter refers to parton-parton scattering in the framework of pQCD, having as basic elements parton evolutions with constant virtuality cutoff Q 2 0 and a hard 2 → 2 elementary QCD scattering.I refer to G = G QCD (Q 2 0 ) as "naive choice".
One recalls that factorization and binary scaling, often mentioned in this paper, amount to reducing the inclusive cross sections for pp and AA scattering to single Pomeron results, although the underlying physical processes involve multiple parallel scatterings.I showed in Sec.IV that neglecting energy conservation leads perfectly to factorization and binary scaling.But in the Monte Carlo procedures eventually one needs the implementation of energy conservation, so one introduces inconsistencies, in the sense that the theoretical basis and the Monte Carlo realization are not compatible.
On the other hand, as shown in Secs.V and VI, considering energy conservation (or energy sharing) and using the "naive choice" G = G QCD (Q 2 0 ), one completely spoils factorization for hard processes, contradicted by data.I have shown that the problem is due to a "deformation" of the inclusive energy distribution of Pomerons connected to many other Pomerons, compared to isolated Pomerons: the probability of carrying a large fraction of the total energy is reduced, which is unavoidable.These deformations can be quantified in terms of deformation functions R deform depending on the number N conn and the squared energy fraction x PE .
In Sec.VII, one takes note of two problems: (1) spoiling factorization when using the naive choice G = G QCD in case of respecting energy conservation, and (2) not considering saturation effects which are known to be important.The solution of these two problems has been shown to be a dynamical saturation scale Q 2 sat , defined via G = k × G QCD (Q 2 sat ) with k being inversely proportional to the deformation function, with a G which must be independent of the connection number N conn .In that case, even having multiple scattering, all inclusive pp and AA cross sections are reduced to a single Pomeron result, but only for hard processes as it should be.This is referred to as "generalized AGK cancellations", which holds at large p t , even in a scenario respecting energy conservation.The dynamical saturation scale works, because even a large number of parallel scatterings will not affect high p t particle production, it will only make the saturation scale big and thus suppress small p t particle production.
Since in the new formalism, the full multiple scattering scenario converges to the single Pomeron result in case of inclusive cross sections (generalized AGK cancellations), one may use the single Pomeron (or factorization) mode, based on EPOS parton distribution functions.So one can now, with the same formalism, treat extremely high p t particle production in factorization mode, and as well collective effects in high multiplicity events using the full simulation.
I discussed several examples, essentially multiplicity dependencies (of particle ratios, mean p t , charm production) which are very strongly affected by the saturation issues discussed in this paper and core-corona effects mixing flow (being very strong) and non-flow contributions.

Figure 3 .
Figure 3. Nonlinear effects: ladders which evolve first independently and in parallel, finally fuse.

Figure 4 .
Figure 4. Nonlinear effects (inside the red ellipses) are "summarized" in the form of saturation scales.

Figure 6 .
Figure 6.Factorization, energy conservation, parallel scattering, and saturation: four concepts that are deeply connected.

Figure 7 .
Figure 7. Double scattering diagram.a double scattering (n = 2), where the blue and green boxes are the elementary T-matrices T (m)

)Figure 8 .
Figure 8. Sum of all possible cuts of a two-Pomeron diagram.

2 Figure 9 .
Figure 9.A simple example of an uncut and the corresponding cut diagram.

Figure 10 .
Figure 10.All the diagrams which contribute to cut T in pp scattering up to order n = 3. Red dashed lines refer to cuts.

Figure 13 .
Figure 13.A Pomeron connected to projectile nucleon i and target nucleon j, together with other Pomerons connected to one (or both) of these nucleons.

Figure 16 .
Figure 16.The contribution G sea−sea QCD , which is the convolutionE soft ⊗ E QCD ⊗ Born ⊗ E QCD ⊗ E soft .

Figure 17 .
Figure 17.Sketch of the suppression of parton yields at high p t with increasing N conn .

Figure 20 .
Figure 20.Sketch of the suppression of low p t partons with increasing N conn .

Figure 22 .
Figure 22.The multiplicity dn/dη(0) as a function of the Pomeron number N Pom for pp (red line) and PbPb (blue line), together with the dotted line representing the function 2.9(N Pom ) 0.9 .

Figure 23 .
Figure 23.The saturation scale Q 2 sat as a function of x PE , for several N Pom event classes.

Figure 24 .
Figure24.The inclusive p t distribution of partons for a full simulation (simu) divided by N coll and the "reference curve" (theo), which is the corresponding distribution for a single Pomeron, calculated analytically.I show results for minimum bias PbPb collisions at 5.02 TeV (red curve) as well as results for different centrality classes.

Figure 26 .
Figure 26.F 2 as a function of x for different values of Q 2 , the latter one indicated (in units of GeV 2 ) in the upper right corners of each subplot.The red curve refers to EPOS PDFs, the green one to CTEQ PDFs, and the black points are data from ZEUS and H1.

Figure 27 .
Figure 27.Parton yield dn/dp t dy for pp at 13 TeV.I show results based on EPOS PDFs (red full line), CTEQ PDFs (green dashed line), the full EPOS simulation (blue circles), and experimental data from ATLAS (black triangles).

Figure 33 .
Figure 33.Transverse momentum spectra of c quarks and charmed mesons in pp at 7 TeV.

Figure 35 .
Figure 35.Two hadron correlations for D 0 D 0 , D 0 D0 , D + D − , and D + D − s as a function of the rapidity difference ∆y in pp collisions at 7 TeV.Red lines represent EPOS4 simulations and black dots data from LHCb.

Figure 36 .
Figure 36.Two hadron correlations for D 0 D 0 , D 0 D0 , D + D − , and D + D − s as a function of the azimuthal angle difference ∆φ in pp collisions at 7 TeV.Red lines represent EPOS4 simulations and black dots data from LHCb.

Figure 37 .
Figure 37. Transverse momentum spectra of prompt J/Ψ in pp at 7 TeV.
Figure 32.Normalized D meson multiplicity as a function of the normalized charged particle multiplicity for different p t ranges in pp scattering at 7 TeV.I show EPOS results with and without hydro, compared to ALICE data.D meson multiplicityd 2 N dydp t / < d 2 N dydp t > as a function of the normalized charged particle multiplicity d 2 N c dydp t / < d 2 N c dydp t > for different p t ranges in pp scattering at 7 TeV, compared to ALICE data