Perturbative QCD concerning light and heavy flavor in the EPOS4 framework

We recently introduced new concepts, implemented in EPOS4, which allow to consistently accommodate factorization and saturation in high energy proton-proton and nucleus-nucleus collisions, in a rigorous parallel scattering framework. EPOS4 has a modular structure and in this paper, we present in detail how the"single scattering module"(the main EPOS4 building block) is related to perturbative QCD, and how these calculations are performed, with particular care being devoted to heavy flavor contributions. We discuss similarities and differences compared to the usual pQCD approach based on factorization.


Introduction
Factorization [1,2] (in connection with asymptotic freedom [3,4]) is a powerful concept to reliably compute inclusive cross sections for high transverse momentum (p t ) particle production in proton-proton (pp) scattering at very high energies.However, 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 highmultiplicity proton-proton collisions [5], which was before considered to be a strong signal for collectivity in heavy-ion collisions.Studying such high-multiplicity events (and multiplicity dependencies of observables) goes much beyond the frame covered by factorization, and we need an appropriate tool, able to deal with multiple scatterings, which must happen in parallel at high energies.
Although everybody will agree that multiple scatterings should conserve energy-momentum, the way it is implemented is fundamentally different in EPOS compared to other models.Concerning momentum sharing, the EPOS4 method employs multiple scattering laws of the form f 1 (p 1 ) × f 2 (p 2 ).. × δ(p − ∑ p i ).The delta function here is crucial, like the probability law in a microcanonical ensemble.We refer to this method as "rigorous parallel scattering scenario", for unbiased momentum sharing.This is very different compared to a structure like f 1 (p 1 |p) × f 2 (p 2 |p − p 1 ) × f 3 (p 3 |p − p 1 − p 2 ) with conditional probabilities (as it is usually done), which may perfectly conserve momentum but in a sequential manner.
In [6], we show that such a "rigorous parallel scattering scheme" can be constructed based on S-matrix theory (see also [7,8,9,10,11]), which we will briefly review in the following.We show -and this is far from trivial -how to accommodate factorization and saturation in an energy conserving parallel scattering scenario.The starting point is the elastic scattering T-matrix (1) expressed in terms of "elementary" T-matrices 1 T Pom , the latter ones representing parton-parton scattering, with the "vertex" V representing the connection to the projectile and target remnants.The elementary T-matrices depend on the light-cone momentum fractions x ± i of the incoming partons, in addition to the energy squared s, and the impact parameter b.The vertices depend on the lightcone momentum fractions of the remnants x + remn (projectile side) or x − remn (target side).Most important: the δ which stands for to assure energy-momentum conservation in an unbiased fashion (this is what we mean by a rigorous parallel scattering scenario).The symbol dX stands for integrating over all these momentum fractions.The expression can be easily generalized for nucleus-nucleus scattering [6].So far we discussed only elastic scattering, the connection with inelastic scattering provides the optical theorem in impact parameter representation, with cut T ≡ 1 i discT, where disc T is the s-channel discontinuity T(s + iǫ) − T(s − iǫ), and with T being per definition the Fourier transformed T-matrix divided by 2s.So we need to compute the "cut" of the complete diagram, cut T, i.e., for pp we need to evaluate cut {V × iT Pom × ... × iT Pom × V} , (4) which corresponds to the sum of all possible cuts, considering, in particular, all possibilities of cutting or not any of the parallel Pomerons.We have finally a sum of products with some fraction of the Pomerons being cut (cut Pomerons), the others not (uncut Pomerons), referred to as cutting rules.We define G to be the cut of a single Pomeron, G = cut T Pom . ( Cut Pomerons represent inelastic scattering ( a + b → many partons) and uncut Pomerons elastic scatterings (a + b → a + b).
The uncut Pomerons are finally summed over and the corresponding variables integrated out, so eventually the total cross section is expressed as a sum over products of G, very similar to Eq. ( 1) with the vertices V(x + remn )V(x − remn ) replaced by some effective vertex function W(x + remn , x − remn ).A great advantage of this Tmatrix formalism is its modular structure.Multiple scattering is expressed in terms of modules (G) representing a single scattering each, as shown in the case of two inelastic scatterings in Fig. 1.The precise content of these modules will be discussed later -this is where all the perturbative QCD physics is hidden, and discussing that will be the main purpose of this paper.
The important new issue in [6] is the understanding of how energy conservation ruins factorization, and how to solve this problem via an appropriate definition of G.The cut Pomeron G = cut T Pom is the fundamental quantity in the EPOS formalism.For the moment, we consider the Pomeron as a "box", with two external legs representing two incoming particles (nucleon constituents) carrying light-cone momentum fractions x + and x − , so G = G(x + , x − , s, b), with the energy squared s, and the impact parameter b, see Fig. 2. Let us define the "Pomeron energy fraction" Pom /s, with M Pom being the transverse mass of the Pomeron, which is the crucial variable characterizing cut Pomerons: the bigger x PE , the bigger the Pomeron's invariant mass and the number of produced particles.Large invariant masses also favor = G(x + , x − , s, b) high p t jet production.We also define, for a given cut Pomeron connected to projectile nucleon i and target nucleon j, a "connection number" N conn = (N P + N T )/2, with N P being the number of Pomerons connected to i, and with N T being the number of Pomerons connected to j.The case N conn = 1 corresponds to an isolated Pomeron, which may take all the energy of the initial nucleons, whereas in the case of N conn > 1 the energy for a given Pomeron will be shared with others.To quantify the effect of energy sharing, we define f (N conn ) (x PE ) to be the inclusive x PE distribution , i.e., the probability that a single Pomeron carries an energy fraction x PE , for Pomerons with given values of N conn .The distribution for N conn > 1 will be deformed with respect to the N conn = 1 case, due to energy sharing, and we define the corresponding "deformation function" R deform (x PE ) as a ratio of f (N conn ) (x PE ) over f (1) (x PE ).As shown in [6], this function can be calculated and tabulated.As discussed in very much detail later, we also calculate and tabulate some function G QCD (Q 2 , x + , x − , s, b), which contains as a basic element a cut parton ladder based on DGLAP parton evolutions [9,12,13] from the projectile and target side, with an elementary QCD cross section in the middle, Q 2 being the low virtuality cutoff in the DGLAP evolution.The latter is usually taken to be constant and of the order of 1 GeV, whereas we allow any value.With all this preparation, we are now able to connect G (used in the multi-Pomeron diagrams) and G QCD (making the link to QCD), as follows: For each cut Pomeron, for given x ± , s, b, and N conn , we postulate: with n being independent on x ± .Most importantly, G does not depend on N conn , but Q 2 sat does, it "compensates" the N conn dependence of R (N conn ) deform .Considering even a large value of N conn , we can prove that the distribution f (N conn ) (x PE ) is proportional to a product of vertex functions times G QCD (Q 2 sat , x + , x − , s, b), and so the deformation cancels out!So, for any number of N conn , the inclusive particle production is governed by a single Pomeron represented by G QCD , which does not depend on N conn , apart from the implicit N conn dependence of Q 2 sat .Consequently, the p t dependence of outgoing partons will be independent of N conn in the hard domain (high p t ), see [6] for more details.Computing a p t distribution according to G QCD amounts precisely to factorization or binary scaling in AA scattering, but -againonly obtained for hard processes.In other words, computing R AA in a full AA simulation, we get unity at large p t .
Let us come back to our earlier statement saying that "energy conservation ruins factorization".The statement actually depends on how one relates G and G QCD .In earlier versions, we adopted what we thought at the time to be "the natural Pomeron definition", but what we now call "the naive Pomeron definition", namely G = G QCD .As discussed in very much detail in [6], this leads unavoidably to the following problem: An inclusive pp p t distribution will be a superposition of contributions for different connection numbers.With increasing connection numbers, these contributions get more and more deformed (suppressed at large x PE ), which corresponds to a suppression of the yields at high p t .Therefore the minimum bias p t distribution will be suppressed at large p t compared to the single Pomeron distribution.On the other hand, factorization requires inclusive cross sections to be given by a single cut Pomeron, since based on this diagram one obtains formulas as f PDF ⊗ σ ⊗ f PDF , corresponding to factorization.
The fact that multiple Pomeron interactions reduce for inclusive cross sections to a single cut Pomeron (leading to factorization if one assumes "the naive Pomeron definition"), refers to "AGK cancellations" [10], shown to be valid in a scenario without energy conservation.As discussed above, including energy sharing (as it should be), ruins first "AGK cancellations" and, as a consequence, factorization.With a "proper Pomeron definition" as employed in EPOS4, we recover "AGK cancellations" in the sense that inclusive cross sections can be expressed in terms of a single Pomeron expression G QCD , although in reality multiple scatterings contribute.But this statement is only true for p t values bigger than the relevant saturation scales of the different multiple scattering contributions, where "relevant" refers to the relative weight of the contributions.
The new EPOS4 framework is able to recover factorization at large p t (a difficult task in the parallel scattering formalism).This allows us to compute, tabulate, and employ "EPOS PDFs", and based on these, we may compute inclusive cross sections as a simple convolution of two PDFs and an elementary pQCD cross section, and the result will be identical -at large p t -to the simulation results.But to make this work in practice, we need high precision and appropriate methods for both simulation and cross section calculations.In older versions, for example, we used in the simulations frequently a "redo" whenever "kinematic problems" showed up, whereas such situations should be avoided.
As a side remark, the EPOS4 multiple scattering scheme is quite different to what is usually called 'multiple parton interactions' (MPI): the former refers to multiple cut Pomerons with external legs being soft partons, whereas the latter treats the scattering of two hard partons, based on multiple parton distribution functions, generalizing the factorization approach.The review [14] has actually two parts, namely "hard MPI" and "soft MPI", and the two have little in common.
The main purpose of this paper is to provide detailed information about the calculation of G QCD , based on perturbative QCD, with special care concerning heavy flavors.We discuss the implementation (for the first time in the EPOS framework) of the "backward parton evolution method", which allows much better control of the hard processes.We discuss important "kinematic" issues connected to processes involving charm and bottom, taking into account 12 different "reaction classes" for the cross section calculations, since the kinematics is quite different, e.g., for the Born processes gg → gg and gg → b b.A useful side-effect concerning the new strategies, in particular the "backward evolution", is the fact that many formulas are very similar to what is used in models based on "factorization".We may compare our EPOS PDFs with "standard" PDFs, but there are also fundamental differences: in EPOS, we have first to deal with evolutions for each parton ladder, with an initial distribution of the corresponding parton distribution of the type δ(1 − x), and these singularities need to be taken care of.Related to this, dealing with singularities is the major challenge for cross section calculations as well as for parton generation.
We mentioned already the modular structure of the approach, where the multiple scattering is expressed in terms of the cut Pomerons G, and the latter one corresponds to some function G QCD .This function G QCD itself has a modular structure, the modules being vertex functions, a soft evolution function, and most importantly the "cut parton ladder", which is (up to an impact parameter dependent factor) equal to a parton-parton cross section.In section 2, we will discuss how the EPOS4 building block G QCD is related to parton-parton cross sections, and in section 3, we will discuss how the integrated and differential parton-parton cross sections can be computed.The former are needed to compute G QCD , while the latter are necessary for the parton generation.

Relating the EPOS4 building block G QCD to parton-parton cross sections
As discussed in the last section (see also [6]), the multiple scattering contributions to the total cross section are expressed in terms of (products of) cut Pomeron expressions G, and the latter ones are related the "real QCD expressions" G QCD via Eq.(6).In this sense, G QCD is the basis of everything, the fundamental building block of EPOS4.This G QCD depends on the saturation scale Q 2 sat , which is of fundamental importance in the EPOS4 framework [6], but here the focus is on the details for the calculation of G QCD , where Q 2 sat is "only" a constant, represent-ing the low virtuality cutoff (and we use simply symbols like Q 2 1 and Q 2 2 as arguments).We will use different symbols, like G, T, T, and σ, all being related to each other.For clarity, we recall the definitions in Table 1.

T
Diagonal element of the elastic scattering Tmatrix as defined in standard quantum mechanics textbooks, where the asymptotic state is a system of two protons or two nuclei T Fourier transform with respect to the transverse momentum exchange of the elastic scattering Tmatrix T, divided by 2s (formulas are simpler using this representation) G Defined as G = cut T = 2ImT = 1 i discT (where "disc" refers to the variable s), referring to the inelastic process associated with the cut of the elastic diagram corresponding to T σ Integrated inclusive parton-parton scattering cross section, which is useful because T, T, and G may be expressed in terms of σ Concerning the structure of G QCD , we expect on both projectile and target side two possibilities, namely a valence quark or a sea quark/gluon being the first perturbative parton, and in addition soft contributions, correspondingly we have In this section, we will discuss the different G's one by one, and in particular how the G K QCD are related to parton-parton cross sections σ, the latter one discussed in great detail in section 3.

The usual factorization approach for proton-proton scattering
Certain elements of the EPOS4 scheme are similar to the usual factorization approach, which we will briefly review in the following.Based on the factorization hypothesis, the inclusive parton production cross section in proton-proton scattering is given as [2] where the parton distribution function (PDF) f a PDF represents the number of partons of species a entering the hard scattering, and where p 1 and p 2 (p 3 and p 4 ) are the momenta of the incoming (outgoing) partons.The PDFs can be considered to be known from deep inelastic electronproton scattering, the amplitudes M can be computed, and thus one obtains the inclusive cross section as a simple integral.The factorization scale µ 2 F represents the scale which allows to separate "long range (soft)" and "short range (hard)" parts of the interaction, such that the former ones are part of the proton structure, whereas the latter ones show up in the perturbative matrix elements M.This procedure has been shown to successfully describe experimental jet production data.In the EPOS4 framework, a similar "factorization formula" will be used but in a different context.

Cut parton ladder G hard
QCD and the partonparton cross section σ ij hard We first discuss a somewhat simplified situation, namely a simple parton ladder.Let us consider the scattering of two partons with flavors i and j carrying virtualities Q 2 1 and Q 2 2 , both first evolving via parton emissions with ordered virtualities up to µ 2 F before interacting as an elementary 2 → 2 pQCD (hard) scattering (see Fig. 3).The corresponding di-jet production cross section can be written as where the momenta of the outgoing partons (jets) are integrated out.This expression is very similar to the usual factorization formula Eq. ( 8), but with the parton distribution functions ), representing parton evolution starting at virtuality Q 2 K with a distribution δ(x − 1)δ ki , but using the same DGLAP evolution [9,12,13].The corresponding cut diagram, referred to as "cut parton ladder", is shown in Fig. 4. Concerning the corresponding elastic scattering T-matrix, we assume with R 2 hard ≈ 0 [9], which is compatible with the usual relation σ ij hard = 2Im T ij hard (t = 0)/(2s).Furthermore assuming a purely transverse momentum exchange t = −q 2 ⊥ , the Fourier transform and division by 2s gives For the corresponding G = cut T hard = 2Im T hard , we get So the cut parton ladder expression G is simply the product of the dijet production cross section σ ij hard times a Gaussian impact parameter dependence.

QCD to the parton-parton cross section σ ij hard
Here we consider the "val-val" contribution, where on both sides a valence quark is the first perturbative parton.We might imagine in the multiple Pomeron formula (as Eq. ( 1), but using T and not T) for each "val-val" Pomeron an expression like where the indices i and j refer to the flavors of the valence quarks, with some "vertex functions" F a , and with s being the proton-proton squared cms energy.However, the external legs of our Pomerons should always be colorless objects, carrying light-cone momentum fractions x + on the projectile and x − on the target side.So each valence quark has a partner (antiquark or diquark), with the valence quark carrying a light-cone momentum fraction z ± and its partner x ± − z ± , the sum of both being x ± .The partner is emitted (as a timelike parton) immediately, and the valence quark starts the parton evolution.The corresponding cut diagram is indicated in Fig. the valence quark and of its partner, and the Mandelstam t.Let us first look at the T-matrix.Instead of Eq. ( 14), we use for each "val-val" Pomeron an expression (including integration) like

We imagine vertex functions
which can then be written as an integration dx + dx − /(x + x − ) referring to the "white" legs, with an integrand that can be interpreted as the corresponding T-matrix, i.e.
with the "val-val" T-matrix = ∑ i,j dz For F i val 1 , F j val 2 , and T ij hard , the t dependence can be factored out as e R 2 t , with parameters R 2 val 1 , R 2 val 2 , and R 2 hard .We compute G as usual as twice the imaginary part of the Fourier transform of the T-matrix divided by 2sx + x − , and we get (using Eqs.(11,12)  For the "sea-sea" contribution, we expect a "soft block" preceding the first pertubative parton, as indicated in Fig. 6.The vertices F sea 1 and F sea 2 couple the parton ladder to the projectile and target nucleons.In addition, we have three blocks, the two soft blocks and in between the parton ladder discussed earlier.The corresponding elastic scattering T-matrix for the latter is T jk hard and those for the soft blocks are per definition T j soft and T k soft , with j and k being the flavors of the two parton which connect the soft blocks to the parton ladder.In principle we have for each of these partons a four-dimensional loop integral, but based on the assumption that transverse momenta and virtualities are negligible compared to longitudinal momenta, they can be reduced to one-dimensional integrals [11],

Relating
What is the reason for getting the imaginary part of the soft T-matrices?The loop integrals may be written as dk + dk − d 2 k t ..., and the k − variable can be related to the Mandelstam variables s and u (for the soft block).Then the usual branch cuts translate into cuts for k − , and a rotation of the integration path, as shown in Fig. 7,

Re k
with as usual T being the Fourier transform divided by 2ŝ = 2x + x − s of the T-matrix, we note that for all T-matrices, the t dependence can be factored out as The t dependence of F sea 1 and F sea 2 can also be factored out as e R 2 sea K t , with parameters R 2 sea 1 , R 2 sea 2 , which gives a simple overall t dependence of the r.h.s of Eq. ( 20) as e R 2 t , so that we get easily and defining so-called "soft evolution functions" we get finally with R being given explicitly as The precise form of E i soft and F sea K will be discussed later, the former is based on a parametrization (in Regge-pole form) of the soft T-matrix, the latter has a simple power law form.Having discussed the two contributions "val-val" and "sea-sea", referring to contributions with respectively two valence quarks and two sea quarks as "first perturbative partons" entering the parton ladder, we easily obtain the expressions for "val-sea", namely with R being given explicitly as The expression for "sea-val" is with R being given explicitly as

The vertices F
The main formulas of the preceding sections, Eqs.(18,28,30,32), allow to express the different G J QCD in terms of "modules", among them the vertices F sea K and F val K with K ∈ {1, 2}, which are simple functions, to be discussed in the following.In the case of pp scattering, the functions for K = 1 and K = 2 are identical, whereas for πp or Kp scattering they are in general different (more precisely, the form is the same but not the parameters).The vertices F sea K and F i val K are given as with with q i val being a standard valence quark distribution function for a small value Q 2 0 of the virtuality (see [11] Appendix C.2). From Eq. (35), and using V(z) = z α remn as remnant vertex, we obtain as parton distribution having used 1 0 dy y −α IR (1 − y) α remn = Γ(1−α IR )Γ(1+α remn ) Γ(2−α IR +α remn ) .So we get f i (z) = q i val (z), as it should be, which justifies our choice of F i val K .

The soft elements E soft and G soft
Based on the asymptotic Regge-pole expression, the soft T-matrix T i soft (Q 2 , s, t) is (for the moment) assumed to not depend on Q 2 and to be proportional to to i s 1+B soft exp (2R 2 soft + α ′ ln s s 0 )t , with a "Regge-pole intercept" 1 + B soft = α soft (0), where B soft is a parameter close to zero.The "soft evolution functions" E i soft (Q 2 , z) are defined in Eq. ( 27) as the imaginary part of T i soft for t = 0 and s = s 0 /z, so they are up to constants equal to (s 0 /z) 1+B soft .We add a splitting into a quark, as well as a soft-Pomeron-ladder coupling of the form (1 − z) C soft , to get (we drop the Q 2 argument) with with a parameter w split ∈ [0, 1], and where we use f split (ξ) = P q g (ξ), since we are "close to the perturbative domain".
Still based on the soft T-matrix expression ∝ i s 1+B soft exp (2R 2 soft + α ′ ln s s 0 )t and adding the same two vertices as for the "sea-sea" contribution, we get the "soft" G expressions as F sea 1 (x + , 0) F sea 2 (x − , 0) with parameters D soft and B soft , and with R being given explicitly as

The pseudosoft elements E psoft and G psoft
The "soft evolution functions" E i soft are meant to take care of the "non-perturbative part", before entering the perturbative regime represented by E ik QCD (x, Q 2 , µ 2 F ). Originally, we were using Q 2 being equal to some soft scale Q 2 0 (typically 1 − 2 GeV 2 ).In that case, E i soft represents purely soft physics, and we could actually drop Q 2 as an argument, E i soft = E i soft (z).However, now we introduce evolution functions E ik QCD (x, Q 2 sat , µ 2 F ) where the virtuality cutoff is equal to the saturation scale, which is in general bigger than Q 2 0 , which means that the "part not taken care of" by the QCD evolution, is not purely soft anymore.This is the domain of gluon fusions, but there should be emissions as well and a narrowing of the momentum fraction distributions.To take this into account we define "pseudosoft evolution functions" as which is meant to replace the soft evolution functions used so far.This is actually needed to have the same p t distribution of emitted partons at large p t for different values of Q 2 sat .Having replaced the soft evolution functions with the pseudosoft ones, we need a modification of G soft as well.We have little guidance from theory, so we use the same parametrization as for G soft just with different parameters B psoft and D psoft , which may depend on Q 2 sat .This is enough to compute G QCD .But because of Q 2 sat > Q 2 0 , we expect that hard pQCD processes should occur, producing light and heavy flavor hadrons.This will be discussed more in the next section.

Parton-parton cross sections in EPOS4 involving light and heavy flavors
As discussed in section 1 (see also [6]), the multiple scattering contributions to the total cross section are expressed in terms of (products of) cut Pomeron expressions G, each one representing a single scattering.They are related to the "real QCD expressions" G QCD via Eq.( 6), in other words, G QCD is the fundamental building block of the multiple scattering framework of EPOS4.We showed in section 2, that G QCD has several contributions, each one being composed of "modules", with the "key" module being the integrated inclusive partonparton scattering cross section σ ij hard , where i and j refer to the flavors of the two partons.This "module" σ ij hard contains all the pQCD calculations.We refer to Table 1 in section 2 to recall the definitions of the symbols G and σ and their relation to the T-matrix.In this section, we discuss in detail how to compute σ ij hard .Concerning notations: whereas in the previous section the symbol s usually referred to the Mandestam s of nucleon-nucleon scattering, we will use in this section s and t referring to the Born process of the hard scattering, which here plays a crucial role.Although it is not always explicitly written, the pQCD matrix elements M are always considered to be given in term of s and t (that is how the are tabulated).

Integrated and differential partonic cross sections
The main formula concerning the parton-parton scattering cross section σ ij hard is Eq. ( 9).It is useful to split that formula, by first considering "differential cross sections" representing the inelastic scattering of two partons with virtualities Q 2 1 and Q 2 2 at a center-of-mass energy squared s lad , with x 1/2 being the light-cone (LC) momentum fractions (with respect to the LC momenta of the partons before the evolution) of the partons entering the Born process of the elementary parton-parton scattering 1 + 2 → 3 + 4, see Fig.The differential parton-parton cross section The solid lines represent partons (quarks, antiquarks, or gluons).and n.The matrix elements M is considered to be given in term of s and t, with s = x 1 x 2 s lad .Both partons first evolve via parton emissions with ordered virtualities up to µ 2 F before interacting as an elementary 2 → 2 pQCD scattering.The "integrated cross sections" σ ij hard may then be written as in terms of the "differential cross section".Both differential and integrated cross sections are important in the EPOS4 framework: the integrated ones are needed to compute the weights for all possible parallel scattering configurations, and to generate the corresponding configurations, the differential cross sections are needed to compute parton distributions, and generate partons in the Monte Carlo mode, for given configurations (generated in the first step).
Comparing Eq. ( 8) and Eq. ( 44), we see that our differential cross section has the same structure as the inclusive cross section based on factorization, but in our case, we use evolution functions E QCD rather than the usual parton distribution functions f PDF , where E QCD represents an evolution also according to DGLAP [9,12,13], but starting from a parton and not from a proton, with an initial condition . Special care is needed to treat this "singularity".
We will eventually use dynamical saturation scales Q 2 sat 1 and Q 2 sat 2 as endpoint virtualities, so we compute (via numerical integration) and then tabulate 2 ) for a large range of possible values s lad , Q 2  1 , and Q 2 2 for all possible combinations of i and j.During the simulations, we use the predefined tables to compute σ ij hard via polynomial interpolation.As discussed in detail later, all the kinematic variables shown in fig.8 are related to each other: The scale µ 2 F and the Mandelstam variable t are related to p 2 t (the transverse momentum of one of the outgoing partons, say parton 3), the end virtualities Q 2 1 and Q 2 2 represent lower limits to p 2 t , but the precise relations depend on the particular Born process.We deal with this problem by introducing classes of cross sections, with the same Born kinematics per class.The main reason that we need to be very careful concerning parton kinematics is the fact that in our case, the differential cross section formula Eq. ( 44) concerns one single scattering in a multiple scattering configuration as shown in Fig. 9.Here we show an example with Multiple scattering configuration for three scatterings.three scatterings, but in AA collisions, we have configurations with 1000 scatterings.So it is out of question to use the standard method of Monte Carlo programming, namely rejection in the case of forbidden kinematics, we simply have to avoid such cases.

Parton evolution functions E QCD
Our parton evolution functions E ik QCD are similar to the "usual" ones, obeying the same evolution equations, but the initial conditions are different, and this requires different treatments.We present here a new procedure to compute and tabulate these functions, much faster than the techniques we used in EPOS3.The indices i and k refer to parton flavors in the form of integers with 0 for a gluon, 1 − 6 for quark flavors from u to t, and the corresponding negative numbers for the antiquarks.
For a given end parton i, using t = Q 2 , the evolution equation for an evolution from t a to t b may be written as We use P = α s 2π P, with P being the splitting functions without () + prescription, i.e., [2] with quark indices q ∈ {−6, ..., −1, 1, ..., 6} and g = 0, with C F = 4/3, T R = 1/2, C A = 3, and where assures a possible emission of a quark of flavor q if the number of flavors N f (t) is at least as big as |q|.So for large values of t, with N f (t) = 5, the weight for emitting a bottom quark is identical to the one of a u or d quark.The symbols ∆ i refer to the so-called Sudakov form factors, defined as and The fundamental difference with respect to the "usual" PDFs is the fact that we have a singular initial condition, namely so we have to define the evolution function E ab QCD (x, t a , t b ) for at least one emission, such that the evolution function can be written as As proven in Appendix A.2, we have for an arbitrary t c between t a and t b the relation and using Eq. ( 55), we get We tabulate E, which means we compute E ab QCD x, t i , t j with t min ≤ t i ≤ t j ≤ t max (blue area in the t i − t j plane, in Fig. 10) explicitly and for i, j from 1 to M. We first compute E ab QCD x, t i , t j for all i and j = 1 and j = 2 in an iterative fashion.These are the red points in Fig. 10.We then compute the E values for j = 3, 4, 5..., always for i = 1 to i = M.We use Eq. ( 57) as This works, since E ab QCD x, t i , t j−1 and E cb QCD x, t j−1 , t j are already known via interpolation from already tabulated values.In this way, we are able to compute and tabulate the evolution functions E ab QCD , and use them via polynomial interpolation, and thus we know E ab QCD via Eq.(55).

Computing integrated parton-parton cross sections with quark mass dependent kinematics
As discussed in the last section, the integrated partonparton cross section σ ij hard is a crucial element of the EPOS4 framework.It may be written as (see Appendix B.2 Eq. ( 243)) where the indices i, j, k, l, m, n refer to all kinds of flavors (gluons and (anti)quarks up to bottom).The matrix elements M is considered to be given in terms of s and t, with s = x 1 x 2 s lad , with s lad referring to the partonparton scattering, see Fig. 8. EPOS was originally constructed only for light flavors, simply considering two types of partons, massless quarks and gluons.Very often, in the program, the two cases were simply treated explicitly, which required major changes in the code structure.
For instance, the number of active flavors (which appears in α s and the splitting function), which was a constant so far, now depends on Q 2 .And there is a second challenge: In the "massless case" of earlier versions, we could first do the sum over all Born processes and compute "the" Born cross section for a given pair of partons k and l, before doing the integrals, which considerably simplifies the calculations.Now we need explicitly the information about the process kl → mn, since the masses involved in the process affect the kinematics (like the relation between Mandelstam t and p t and integration limits).So the operations ∑ mn and dx 1 dx 2 dt can no longer be exchanged in general, but it can be within classes: We introduce classes K of elementary reactions kl → mn obeying to the same kinematics and then compute all relevant quantities (to be discussed in detail in the next section) depending on of K.In this way we do the kinematics properly and still keep efficient (fast) procedures, which is crucial in the EPOS4 framework.

Born kinematics
In the last section, we were considering processes involving initial partons i and j and two partons k and j entering the Born process (let us note them 1 and 2), producing partons m and n (let us note them 3 and 4), so altogether we have the Born process 1 + 2 → 3 + 4. In the CMS, we write the parton momentum four-vectors as with p t = 0 for the incoming ones.So per definition, E is the modulus of the momentum, E = | p|, in the CMS, and so we have obviously for all four particles and In the following, E (without index) refers to E 1 (= E 2 ) and E ′ to E 3 (= E 4 ).We will use the following definitions: From energy conservation, we find which gives In case of In case of and in case of m 1 = m 2 = 0, we get Corresponding formulas apply for W ′ .Essentially all kinematic relations will be expressed in terms of W and W ′ .

Constraints for p t
Crucial for the following discussions is the relation between the factorization scale and the transverse momentum of the outgoing parton, or it's inverse with M being the (maximum) mass of the partons involved in the Born process, and where the two coefficients κ and λ represent the freedom in defining µ 2 F p 2 t .In EPOS4.0.0, we use κ = 1 and λ = 0.A choice κ = 1 and λ = 1 makes little difference (and only at very small p t ) for charm production, but for bottom one should better use this choice, as we will do in future releases.
The p 2 t values in the integration in the cross section formulas of the last section are restricted, since we have which amounts to In reality, we have to use the Mandelstam t rather than p 2 ⊥ in the integration, so we need to find the limits for this variable.

Constraints for s
The quantity E ′ is precisely the upper limit for the transverse momentum of particle 3, i.e., To get non-zero results, we need which gives and solving the quadratic inequality equation, we get t min .This is first of all a limit for the energy squared of the Born process, but of course as well a limit for the ladder.

Constraints for t
Concerning t, we have which leads to the "∓" referring to respectively p z ≥ 0 and p z ≤ 0. We may invert Eq. ( 83), to obtain which allows us to compute p 2 t for a given t.From Eq. (83), we get Instead of integrating from |t| min to |t| max , one may define the maximum value |t| max+ of |t| with p z ≥ 0, i.e., ) This quantity |t| max+ is the upper limit of |t| actually used, since we may write (87) and a variable transformation t = 2|t| max + − t ′ for the second integral (and then replacing t ′ by t) leads to since the upper limit of the transformed variable in the second integral is 2|t| max + − |t| max + = |t| max + and the lower limit is (using |t| min + |t| max = 2|t| max + ) given as 2|t| max + − |t| max = |t| min .

Cross section classes
From the discussion in the last sections, it is clear that when computing cross sections, we need to identify explicitly contributions from particular classes of elementary scatterings, obeying different kinematics, due to the different masses being involved.We will use the notation representing the reaction 1 + 2 → 3 + 4 with the four masses being m 1 , m 2 , m 3 , and m 4 .The masses for charm and bottom quarks used in EPOS4.0.0 are m c = 1.27GeV/c 2 and m b = 4.18 GeV/c 2 .We will distinguish 12 different classes ("light" refers to massless partons) as follows: 1. ll → ll : scattering of light partons (case 0000) The cross section calculation [see Eq. ( 61)] is then given as where "klmn ∼ K" refers to indices corresponding to the class K, where one then uses the integration limits and kinematic relation according to K.For numerical efficiency, we use the explicit formulas from Appendix C.This is still not yet the final formula, as discussed in the next section.

Cross sections for both-sided, single-sided, and no emissions
Our evolution functions are based on the same equation as the usual PDFs, but we have a singular initial condition, i.e., E ik QCD (x, t a , t a ) = δ(1 − x)δ ik .Therefore, as discussed in section 3.2, we compute and tabulate ẼQCD , where ẼQCD refers to the case of at least one emission, such that This takes explicitly care of the singular initial condition.Correspondingly, we have to deal with four different integrated parton parton cross sections.First, we have which represents the integrated parton-parton cross section with at least one emission on each side, referred to as both-sided emissions.Then we have single-sided emissions, i.e., representing the integrated parton-parton cross section with at least one emission on the upper side and no emission on the opposite side, and representing the integrated parton-parton cross section with at least one emission on the lower side and no emission on the opposite side.Finally, we have representing the integrated parton-parton cross section with no emission on either side.The four formulas are very similar, just in case of no emission, one replaces Ẽjl QCD by δ jl ∆ j and one has therefore one summation less.
In principle there is the possibility, for all these cross sections, to add a so-called K-factor to compensate higher order effects, but presently we use K-factor = 1 (in EPOS4.0.0), since we use already a "variable flavornumber scheme", i.e., the number of active flavors (in α s , in Ẽ etc) depends on the virtuality: In that case, LO calculations with a K-factor being unity already give a fair description of the data, as has been demonstrated independently in [15] (figure 4) and [16]. We As discussed in section 2, the cross section σ ij hard may then be used to compute "the EPOS4 building block" G QCD , the basic quantity in our parallel scattering scheme.

Taming singularities in the cross section integrals
For all integrations in the last section, we need to be aware of the fact that the integrands contain "singularities" like 1/(x − a) λ .Even when the integration domains avoid such singularities, as b a+ǫ f (x)dx with ǫ ≪ 1, the numerical integration (Gaussian integration in our case) may give wrong results, if it is employed in a naive fashion.One actually needs to properly transform the integration variables to have smooth integrands.Let us consider an integral b a g(x)dx, which may be transformed to which may be written as if the primitive G of g is known.
Then we get u u(a) g(x(u ′ )) dx du du ′ ∝ (u + 1)/2, which means the integrand is constant.If we have to compute b a f (x)dx, we are looking for some "simple function" g(x), being close to f (x), in order to use Eq.(98) to define the coordinate transformation u ↔ x.We then compute where we expect the integrand f (x(u)) dx du to vary slowly with u (since for f = g it would be constant).Then based on the u-integration variable, we use 14-point Gaussian integration.It is not necessary to find a g very close to f , it should be "sufficiently close" to give a well-behaved transformed function, in particular taking care of singularities.In the case of several singularities, one needs to split the integration domain and use an appropriate coordinate transformation for each piece.In the case of singularities, there are always cutoffs which make the integrals mathematical well defined, it is a purely numerical issue, without transformations, the integration results could be completely wrong.
To compute σ both ij hard according to Eq. ( 92), we first note dx 1 dx 2 dt = dz dtdx 1 /x 1 with z = x 1 x 2 .Concerning the z integration, we make a coordinate transformation Eq. (98) with because this gives (for f ≈ g) which is what we expect.Actually σ both ij hard ∝ s δ hard defines the empirical technical parameter δ hard (the numerical value is 0.25).This case shows the importance of the transformation to take care of a singular integrand, here of the form δ hard z −δ hard −1 , which allows nevertheless numerical integration with great precision.As mentioned before, it is not important to have an "approximation" g(z) everywhere close to f (z), but it must take care of the singularity.Next, we consider (inside the z integral) the t integration.Here we make a coordinate transformation Eq. (98) with which corresponds to g(t) = 1/t 2 , the expected t singularity from the Born process.Finally (inside the z and the t integral), we consider the x 1 integration.Here we split the integration domain, say a → b, into a → c and c → b, since we expect singularities 1/(1 − x 1 ) and 1/x 1 , from Ẽij QCD (x 1 ) ∼ P j i (x 1 ) close to x 1 = 1 and x 1 = 0. Concerning the integration c → b, we make a coordinate transformation Eq. (98) with which corresponds to g(x 1 ) = 1/(1 − x 1 ), corresponding to the expected x 1 singularity for x 1 close to unity.For the integration a → c, we make a coordinate transformation Eq. (98) with G(x 1 ) = ln(x 1 ), which corresponds to g(x 1 ) = 1/x 1 , being the expected x 1 singularity for x 1 close to zero.
To compute σ upper ij hard according to Eq. ( 93), we have only a x 1 and a t integration, which are treated as the corresponding integration for σ both ij hard .The cross section σ lower ij hard is not computed but derived from σ upper ij hard as i.e., by exchanging i and j as well as Q 2 1 and Q 2 2 .Finally for σ none ij hard according to (95), one needs only a t integration, also done in the same way as in the other cases.

Differential parton-parton cross sections and backward evolution
One of the (new) key elements in EPOS4 is the fact that even though we have in general to deal with multiple parallel ladders (see Fig. 9), which energy-momentum shared between them, it is possible -for each ladderto first generate the Born process (magenta dot in Fig. 9) and determine the corresponding outgoing partons, and then via backward evolution generate the parton emissions.This has many advantages (compared to the forward evolution technique in EPOS3): Not only does it correspond to the "usual" method employed by factorization models, it also allows much better control of the hard processes.We consider parton-parton scattering as shown in Fig. 11.Starting from the differential parton-parton cross sec- tion Eq. ( 44), as shown in Appendix B Eq. ( 244), we may write Compared to Eq. (244), we added a "∑ 12 K=1 " term and the "klmn ∼ K" condition to account for 12 different cross section classes, as discussed in section 3.3.5.The variables x 1 and x 2 are the momentum fractions of partons 1 and 2. This formula is very useful as a basis to generate the hard process in the parton ladder, in the Monte Carlo procedure, when employing the backward evolution.
However, we should not forget that we always need to distinguish four cases, namely both-sided, onesided(lower), one-sided(upper), or no emissions.Correspondingly, we have four equations.For both-sided emissions, we use Eq. ( 106), with Ẽ instead of E. In case of one-sided emission, we replace one of the Ẽjl QCD by δ jl ∆ j , and in case of no emissions, we replace both Ẽjl QCD by δ jl ∆ j .In the case of both-sided emissions, we define, based on Eq. ( 106) with Ẽ instead of E, the function X ij klmn as

Generating the
We want to generate for given i and j the variables z, t and x 1 , according to the probability distribution which is not trivial due to the fact that the expressions X ij klmn contain singularities, as discussed in section 3.3.7,where we discussed the calculation of σ both ij hard as a sum of integrals over X ij klmn .We also showed how to solve this problem.We had integrals of the form X ij klmn dz dt dx 1 , which could be done after three coordinate transformations.Let us name the three integration variables z 1 , z 2 , and z 3 , i.e., The three coordinate transformations z i → u i were defined as for i ∈ {1, 2, 3}, with There are two transformations for z 3 , due to two singularities (at 0 and 1), which required a splitting of the integral.The three functions g i are useful also for the Monte Carlo generation of the variables z i : One may first generate the z i according to g i , by first generating uniform random numbers u i between -1 and 1, and then computing z i = z i (u i ) by inverting Eq. (112) (which is easy).Then, one accepts these proposals with the probability where M has to be chosen such that p accept ≤ 1.In this ratio, the singular parts of prob(z 1 , z 2 , z 3 ) are canceled out by the g i , and p accept is therefore a smooth function of the variables, with an acceptable rate of rejections.Then, for given i, j, z 1 , z 2 and z 3 , we generate K and the flavors k and l according to and for given i, j, z 1 , z 2 , z 3 , K, k, and l, we generate m and n according to The two factors Ẽik QCD (...) and Ẽjl QCD (...) do not depend on m or n and could be dropped but for the numerical procedures it is easier to keep them, since we may call the same functions as in the steps before.
In the case of one-sided emissions (on the upper side), we define We want to generate for given i and j the variables x 1 and t, according to the probability distribution Again, we use the same coordinate transformations as already used for the integrations to compute σ upper ij hard .Here, with the transformations z i → u i were defined via Eq.( 112) for i ∈ {1, 2}, with We first generate the z i according to g i (first generating u i , then inverting Eq. ( 112)) and then accept these proposals with the probability Then for given i, j, z 1 , and z 2 , we generate K and the flavor k according to and finally, for given i, j, z 1 , z 2 , K, and k, we generate m and n according to In case of one-sided emissions on the lower side, we use the same algorithm as for the upper side, just exchanging the upper side and lower side variables and indices.
In case of no emissions, we define We want to generate for given i and j the variable t, according to the probability distribution which can be done by defining z 1 = t, and then defining the transformation z 1 → u 1 via Eq.( 112) for i = 1, with g 1 (z 1 ) = z −2 1 .We first generate z 1 according to g 1 (first generating u 1 , then inverting Eq. ( 112)) and then accept these proposals with the probability Then for given i, j, and z 1 , we generate K according to and finally, for given i, j, z 1 , and K, we generate m and n according to (132)

Backward evolution
Knowing for given end flavours i and j, the Born process variables x 1 , x 2 and t, and the flavors k, l (ingoing) and m, n (outgoing), we generate the parton emission via backward evolution.The variable t (Mandelstam variable) allows computing p 2 t [see Eq. ( 84)], which corresponds in a unique fashion to some factorization scale µ 2 F , which is the starting value of the virtuality for the backward evolution.We therefore define Q 2 0 = µ 2 F .Be k the flavor of the corresponding parton.
In the following, we will use the symbol "t" as virtuality, to treat the parton evolution, it does not refer to the Mandelstam variable.We define t 0 = Q 2 0 .We will in the following consider the evolution from t 1 to t 0 with t 1 < t 0 .We write the evolution equation [see Eq.( 46)], using The integrand of the second term corresponds to a last branching in [t, t + dt] , so the corresponding probability is Using Eq. (169), we get The probability of a last branching after some t < t 0 is t 0 t g(t ′ )dt ′ , which means that t is generated via with r being a uniform random number.The integral can be easily done, and we get which leads to Using 1 − r instead of r, we obtain finally for the generation of t in the interval [t 1, t 0 ].So far we essentially followed the standard procedures, as explained in [2].What makes things more difficult in the EPOS4 framework is the fact that we do evolutions for each of the parallel parton ladders, see Fig. 9, and there we consider the evolution starting from a parton, not from a proton.The initial condition for the evolution is and this singularity δ(1 − x) needs some special attention.In Eq. ( 139), we need the full evolution function, which may be written as which separates the "smooth" and the "singular" part.For the generation of t based on R = r, we consider the two cases x = 1 and x = 1 separately.
For the case x = 1, we have R = 1, which simply reflects the fact that there is no emission between t 1 and t 0 , so we recover the initial condition x = 1, and we are done for this case, the emission process is finished.
For the case x = 1, we have which needs to be solved to get t.Some elementary root finder (like the bisection method) will do the job.Once this is done, for this value of t, the probability distribution for the momentum fraction z of the last branching is obtained from Eq. (134) as which gives ) for z ≥ x.We need to distinguish again two cases, z > x and z = x.The latter corresponds to a parton with a momentum fraction unity before the splitting (which is the initial value), so the emission process will be finished at the next iteration step.The probability for the z = x case is , (146) where the integral can be done using the Gauss-Legendre method after an appropriate coordinate transform, splitting the integration domain into two parts, to treat separately the z → 0 (for small x) and z → 1 regions.This probability increases when (during the iteration) the virtuality approaches t 1 and the momentum fraction approaches unity.With the probability W 2 = 1 − W 1 , we have the case z > x, and here the probability distribution to generate z is Finally, we generate the flavor a according to for the given values of t and z found earlier.
At each iteration step, we have these two cases, z = x and z > x, where in one case the iteration is finished in the next step, with the initial value x = 1, whereas in the other case, the iteration continues.So in any case, we eventually end up with the initial values.

Monte Carlo results versus cross section formulas for a single Pomeron
An elementary distribution (for tests) is the transverse momentum distribution dn/dy dp t for a given rapidity of primary partons, directly emitted in the Born process.Let us consider a single Pomeron, carrying the full energy, i.e., x + = x − = 1, for an energy E = √ s = 1 TeV, and let us consider . We first investigate the "sea-sea" contribution.Combining Eq. (28) and Eq.(256), after integrating out the impact parameter b, we get (with a normalization constant N) √ s lad e y and x 2 = x 1 x p t √ s lad e −y .To compute M, we need s = x 1 x 2 s lad and t = −p t x 1 √ s lad e −y .
In a similar way we get explicit formulas to compute the "val-val", the "val-sea", and the "sea-val" contribution.In Fig. 12, we show transverse momentum distributions  Transverse momentum distributions of primary partons for Pomerons of type "val-val", "val-sea", "sea-val", and "sea-sea".We compare the Monte Carlo results (points) to the corresponding curves obtained with the help of explicit cross section formulas (see text).
(for y = 0) of primary partons for Pomerons of type "valval", "val-sea", "sea-val", and "sea-sea", based on these explicit formulas (lines), compared to the Monte Carlo results (points), using for the latter the methods explained in section 3.4.1.The Monte Carlo results agree with those based on cross section formulas.This may sound trivial, but in older EPOS versions with forward parton evolution, this was not the case.A great advantage of the new method is simply the fact that one can make rigorous tests, since the Monte Carlo uses the same "modules" as the cross section formulas (E soft , E QCD , and M, all of them tabulated and available via interpolation).In both cases one uses evolution functions E QCD representing all parton emissions, whereas in the old method the emissions in the Monte Carlo have been done one by one, and one needs to worry about additional elements like the reconstruction of the momentum four-vectors of the emitted partons.

EPOS PDFs
We may go one step further and provide formulas for inclusive momentum distributions for pp scattering with one single Pomeron (identical to the full multiple scattering results at high p t , where factorization applies), again first considering "sea-sea".Combining Eqs.(4,5,28,44), after integrating out the impact parameter b, we get with s = ξ + ξ − s pp being the squared energy of the Born process, with For "valval", combining Eqs.(4,5,18,44), we get with s = ξ + ξ − s pp being the squared energy of the Born process, with ξ 1 = z + x 1 and ξ 2 = z − x 2 .Corresponding formulas can be found for the "val-sea", and the "seaval" contribution (we do not consider "psoft" here, since it only contribute at low p t ).As a consequence, the complete momentum distribution (the sum of the four contributions) may be written as These are the "EPOS4 PDFs" (parton distribution functions), composed of a "sea" contribution (the first integral) and a "valence" contribution (the second integral).Eq. ( 153) is identical to the "usual" factorization formula, which everybody uses.However, in our case, it represents the single Pomeron case.The full simulation, i.e., the full multiple scattering case, provides the same result, but only at large transverse momenta, as we will show later.
To make it very clear: In the EPOS4 formalism, we cannot use PDFs as an input, we need several "modules" such as vertex functions (V, F sea , F i val ) and evolution functions (E i soft , E ik QCD ) to construct the "building blocks" G QCD , and then G (cut single Pomeron), being the basis of our multiple scattering scheme.But we can use these "modules" to construct EPOS4 PDFs.
Let us first compare the EPOS4 PDFs with other choices and with data.At least the quark parton distribution functions can be tested and compared with experimental data from deep inelastic electron-proton scattering.The leading order structure function F 2 is given as where p is the momentum of the proton and q the momentum of the exchanged photon.In Fig. 13, we plot 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 sub-plot .The red curve refers to EPOS PDFs, the green one to CTEQ PDFs [17], and the black points are data from ZEUS [18] and H1 [19,20,21].The two PDFs give very similar results, and both are close to the experimental data.
Having checked the EPOS PDFs, we will use these functions to compute the jet (parton) cross section for pp at 13 TeV, using Eq. ( 153), integrating out the momentum of the second parton and the azimuthal angle of the first parton, which finally gives (see Eq. ( 256), with E ik QCD replaced by f k PDF , and s lad by s pp ) with with {...}being the form the squared matrix elements are usually tabulated, with α s = g 2 /4π.We define the parton yield dn/dp t dy as the cross section dσ/dy dp 2 t , divided by the inelastic pp cross section, times 2 p t , showing the result in Fig. 14.We show results based on EPOS PDFs (red full line), CTEQ PDFs [17] (green dashed line), the full EPOS simulation (blue circles), and experimental data from ATLAS [22] (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.Let us consider the production of charmed primary partons (c and c) for pp at 13 TeV, still based on EPOS4 PDFs, according to Eq. ( 158).The term "primary" means that here we do not consider charm created in the timelike cascade (in the full simulation we do, of course).We have the elementary Born scattering kl → mn, where k, l are the incoming and m, n the outgoing parton flavors.Let us note light flavor partons as "L" and heavy flavor ones (here charm) as "H".We may produce charm via HL → HL ("flavor excitation") or via LL → H H. We will in the following consider two cases: including flavor excitation (incl FlavEx) or not (w/o FlavEx), for two calculations: based on EPOS4 PDFs and based on CTEQ PDFs [17], see Fig. 15.The EPOS4 results are shown as green (incl FlavEx) and red lines (w/o FlavEx), the CTEQ based results are shown as yellow (incl FlavEx) and blue  lines (w/o FlavEx).The EPOS4 and the CTEQ results are similar, at large p t they are even identical.In both cases, the flavor excitation contribution is largely dominant, over the whole p t range, above some threshold.In Fig. 16, we plot the same four curves (EPOS4 and the CTEQ based, with and without flavor excitation), in a re- duced p t range, together with FONLL and NLO calculations [23], the latter represented by black circles (FONLL) and gray squares (NLO).The NLO results (at least at large p t ) are quite close to the ones based on EPOS4 and CTEQ PDFs.In the EPOS4 framework, we use a k-factor equal unity; however, we use a "variable flavor number scheme", i.e., the number of allowed flavors depends on the virtuality, which means at large p t , which is correlated with large virtualities, we easily produce heavy flavor.Concerning the FONLL results, they drop with increasing p t considerably below the EPOS4 and the CTEQ based results.In order to be close to FONLL, we multiply the flavor excitation contribution by a factor k FE < 1, the result for a numerical value of 0.4 is shown in Fig. 17.
The PDFs are actually a sum of "sea" and "val".There is, however, also a "pseudosoft" contribution, see Eq. (43).We expect the "pseudosoft" Pomeron to have an internal structure, allowing hard processes.We assume some probability distribution P(X) with , with Q 2 being the scale associated with the hard process, and then for given Q 2 , a probability distribution of the form to generate a particular hard process.We expect the production of partons with transverse momenta of the order of Q 2 sat , in the range 1-10 GeV/c.

A result for full EPOS
In the examples discussed in 3.4.4,based on EPOS PDFs, only primary partons were considered (directly produced in the Born process).To get the complete picture (but without considering secondary interactions as hydro evolution and hadronic rescattering), we employ the full EPOS4 approach for primary scatterings, including multiple scattering, and also including the timelike cascade, which includes Born processes like g + g → g + g, where each of the outgoing gluons may split into a c c pair.In Fig. 18 (upper left panel), we show the corresponding transverse momentum distribution of charmed quarks (including antiquarks) in pp scattering at 7 TeV, compared to FONLL results [23].The yellow dashed line represents the "pseudosoft" contribution, which visibly contributes at low transverse momentum, as expected.
In the full (primary scattering) approach, we consider not only the full partonic evolution (initial-state and finalstate cascade) but also hadronization.This will be discussed in detail in section 4, but let us anticipate here some basic features: 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 partons, and g to gluons.The Born process or branchings in the spacelike or the timelike cascade may lead to Q Q production, where Q refers to "heavy flavor" (HF) quarks, i.e., charm or bottom.In this case, we end up with parton chains of the type q − g − ... − g − Q and Q − g − ... − g − q, which will decay (among others) into HF hadrons.In Fig. 18, we show transverse momentum spectra of D + , D − mesons (upper right), D 0 mesons (lower left), and D + * , D − * mesons (lower right) in pp collisions at 7 TeV.The red lines refer to EPOS4 simulations, the green points to FONLL calculations [23], and the black points to ALICE data [24].We also show as blue  lines the EPOS4 results without flavor excitation.Many more results can be found in [6].

From partons to strings: color flow
In the previous section, we were discussing in detail the partonic structure of the EPOS4 building blocks, the Pomerons.This allows us to compute the weights of multiple Pomeron configurations and generate these.In the present section, we will discuss how to transform these partonic structures into strings.The strings are then the basis of hadron production via string decay or of the formation of a "core" (in case of high densities), serving as an initial condition of a hydrodynamical evolution.The formation of strings for a particular partonic configuration is based on the associated color flow.

An example
In order to be as general as possible, let us consider a collision of two nuclei, as shown in Fig. 19, where we consider the partonic configuration of two colliding nuclei A and B, each one composed of two nucleons.Dark blue lines mark active quarks, red dashed lines active antiquarks, and light blue thick lines projectile and target remnants (nucleons minus the active (anti)quarks).We have two scatterings of "sea-sea" type, and one of "val-sea" type.We consider each incident nucleon as a reservoir of three valence quarks plus quark-antiquark pairs.The "objects" which represent the "external legs" of the individual scatterings are colorwise "white": quark-antiquark pairs in most cases as shown in the figure, but one may as well have quark-diquark pairs or even antiquark-antidiquark pairs -in any case, a 3 and a 3 color representation.Let us for simplicity consider the quark-antiquark option, and first of all look at the "sea" cases (on the projectile or the target side).In each case, a quark-antiquark pair is emitted as final-state timelike (TL) parton pair (marked 1,2,3,4,5) and a spacelike (SL) "soft Pomeron" (indicated by a thick cyan line), which is meant to be similar to the QCD evolution, but emitting only soft gluons, which we do not treat explicitly.Then emerging from this soft Pomeron, we see a first perturbative SL gluon (another possibility is the emission of a quark, to be discussed later), which initiates the partonic cascade.In the case of "val", we also have a quark-antiquark pair as external leg, but here first an antiquark is emitted as TL final particle (marked 6), plus an SL quark starting the partonic cascade.
In the case of multiple scattering as in Fig. 19, the projectile and target remnants remain colorwise white, but they may change their flavor content during the multiple collision process.The quark-antiquark pair "taken out" for a collision (the "external legs" for the individual collisions), may be u − s, then the remnant for an incident proton has flavor uds.In addition, the remnants get massive, much more than simply resonance excitation.We may have remnants with masses of 10 GeV/c 2 or more, which contribute significantly to particle production (at large rapidities).
In the following, we discuss the color flow for a given configuration, as for example the one in Fig. 19.Since the remnants are by construction white, we do not need to worry about them, we just consider the rest of the diagram.In addition, colorwise, the "soft Pomeron" part behaves as a gluon.Finally, we use the following convention for the SL partons which are immediately emitted: We use for the quarks an array away from the vertex, and for the antiquarks an array towards the vertex.The diagram equivalent to Fig. 19  terings.Horizontal lines refer to TL partons, which later undergo a timelike cascade, while the vertical lines refer to spacelike intermediate partons.We added integers just to mark the different TL partons.For the leftmost scattering, starting from one "end", say "1", we follow the color flow to "5", and then starting from "6" to "10", so we get two chains: 1-2-3-4-5 and 6-7-8-9-10.The end partons of each chain are always quarks or antiquarks, and the inner partons are gluons.Similar chains are obtained for the second scattering, 11-12-13-14-15 and 16-17-18-19-20, and for the third scattering, 21-22-23-24 and 25-26-27-28-29.These chains of partons will be mapped (in a unique fashion) to kinky strings, where each parton corresponds to a kink, and the parton four-momentum defines the kink properties, as already done in earlier EPOS versions [11].In the following, we simplify our graphs, instead of Fig. 21 we use a diagram as shown in Fig. 22, where we do not plot the gluon lines explicitly, since the double lines with arrows in opposite direction allow us perfectly to identify the gluons.
Considering these examples, it is useful to define "initial partons", being those which are immediately emitted as TL partons before the parton evolution starts.In the case of the type "sea", a quark and an antiquark (or diquark) are emitted, in Fig. 21 the partons 1,10 and 5,6 and 11,20 and 15,16 and 24,25.The partons like 2,3,4 or 7,8,9 are emitted either in the spacelike partonic cascade (like 2,4,7,9) or in the Born process (3,8).In case of "val", a TL antiquark is emitted, namely the antiquark 21.Here, an SL quark enters the partonic cascade (vertical line).Parton 29 is not an initial parton, it is the first parton emitted in the SL cascade.
The above discussion demonstrates for an example, how to obtain "chains of partons" (and then kinky strings), based on an analysis of the color flow.But the picture is not yet complete.In general, each emitted parton initiates a timelike cascade, as shown in Fig. 23, and   (a) (b) this has to be taken into account when constructing the "chains of partons" based on color flow.In the following, we discuss the general rules, applicable for any multiple scattering configuration, in pp or AA, including timelike cascades.

The initial TL partons and the initial SL partons
Even in the case of multiple scatterings, the projectile and target remnants remain color white and fragment independently, and correspondingly the external legs of the individual scatterings are white.They are quark-antiquark pairs (or less likely quark-diquark pairs).These external legs emit immediately TL partons, referred to as "initial partons", as compared to the partons emitted later during the parton evolution.The number and the properties of the initial partons for each individual scattering depend on the Pomeron type.
We have four scattering (or Pomeron) types: "sea-sea", "val-sea", "sea-val", and "val-val", where the first concerns the projectile, the second the target.In the case of "sea" (on the projectile or the target side), we have three cases: "sea(g)" and "sea(q)", and "sea( q)".The notion of "sea" always means that there is first a "soft preevolution".In case of "sea(g)" the first parton entering the spacelike (SL) cascade, is a gluon, in case of "sea(q)" a quark, and in case of "sea( q)" an antiquark.These partons are referred to as "initial SL partons".In Fig. 24, we sea(g) sea(q) sea(q) val 2 Figure 24: Initial TL partons 1,2,3, for the cases "sea(g)" and "sea(q)", "sea( q)", and "val", on the projectile side.The "initial SL partons" which initiate the SL cascade are marked A, B, C, D.
show the four possibilities, for the projectile side.In case of "sea(g)", the "external leg" is a quark-antiquark pair, which emits a TL quark and a TL antiquark (2 and 1 in Fig. 24) and an SL gluon, which starts the SL cascade and is color connected to 1 and 2. In the case of sea(q), in addition to 1 and 2, the TL antiquark 3 is emitted, and an SL quark is initiating the cascade, color connected to 1.The TL antiquark 3 is color connected to the TL quark 2 (so 2 and 3 constitute a "little" string, which will be given a small amount of energy).In case of sea( q), in addition to 1 and 2, a TL quark is emitted (3), color connected to the TL antiquark 1.Here an SL antiquark is initiating the SL cascade, color connected to 2. Finally, in the case of val, a TL antiquark is emitted (1), and we have an SL quark initiating the SL cascade, color connected to antiquark 1.
The discussion for the target side is identical, we just use a somewhat different graphical representation, as shown in Fig. 25.Here the arrows from 1 (emitted antiquark) to 2 (emitted quark) go from right to left.Knowing the color connections of the initial partons among each other and with respect to the SL partons initiating the SL cascade, we need as a next step to extend the construction of the color connections to the SL cascade.sea(g) sea(q) sea(q) val Figure 25: Initial TL partons 1,2,3, for the cases "sea(g)" and "sea(q)", "sea( q)", and "val", on the target side.The "initial SL partons" which initiate the SL cascade are marked A, B, C, D.

The spacelike (SL) parton evolution
Whereas for technical reasons, we employ (for each individual scattering) a backward evolution method to generate the kinematical variables of the partons, we will here consider the SL parton evolution starting from the "initial SL partons" towards the Born process.To have a unique definition of the color flow in the direction of the evolution, we define (instead of "right" and "left") the terms "color side" and "anticolor side", as shown in Fig. 26.In this way, in the direction of the evolution, "color side" corresponds to "right", and "anticolor side" to "left".For the evolution of an SL gluon, the two arrows representing the color flow are always such that the forward arrow (in the direction of the evolution) is on the "color side", and the backward arrow on the "anticolor side", as shown in Fig. 27.The aim is to construct the color flow, i.e., sequences of TL partons j 1 − j 2 − ... − j n (like the sequence 1 − 2 − 3 − 4 − 5 shown in in Fig. 22).To do so we number all the timelike partons j = 1, 2, 3... (starting with the initial ones).Referring to this parton numbering j, we define an integer array, the connection array n η CJ (i, j) for TL partons j, defined such that defines the TL parton in the chain just before and the TL parton just after the parton j.So for the example of the chain 1 − 2 − 3 − 4 − 5 shown in in Fig. 22, considering the subsequence 2 − 3 − 4, for j = 3, we have When we talk about sequences j 1 − j 2 − ... − j n , we have two options (which we use both): following the color flow (→) or the anticolor flow (←).In the first case, we use η = 1, in the second one η = 2.The index η is referred to as "color orientation".At the end of the emission procedures, this array will contain the complete information about which partons are color connected.
For the evolution process, it is useful to define the "current SL parton" (say parton a), based on which the emission process is considered (a → b + c), with an SL parton b and a TL c.In the next step, parton b is considered to be the current SL parton, and so on.During the evolution procedure, we use for the current SL parton some integer array, the connection array n CC (i, k) for the current SL parton, providing information to which TL partons the current SL parton is connected to.Here, refers to the connection on the projectile side and to the connection on the target side.For gluons, the index i = 1 refers to the color side, and i = 2 to the anticolor side.In the case of quarks and antiquarks, only i = 1 will be used, and the value for i = 2 is zero.
Consider for example the first SL parton on the projectile side, as shown in Fig. 28.The first "current parton" is a gluon.It is color connected to 1 (color side) and 2 (anticolor side), so we have n CC (1, k) = 1 and n CC (2, k) = 2, here for k = 1, since we consider the projectile side.
In the next step, let us assume we have a gluon emission (parton 3).Concerning color flow, we have two possibilities, an emission on the color side or the anticolor side, as shown in Fig. 29.In the first case, we have a "subchain" 1 → 3, i.e., for the new TL parton j = 3, we have n (→), we use η = 1.In the second case, we have a "subchain" 2 → 3, i.e., for the new TL parton j = 3, we have n η CJ (1, j) = 2, and since we follow the anticolor flow (→), we use η = 2.In other words, the two choices correspond to the two choices η = 1 and η = 2, and they are chosen randomly.
Having discussed the case of a gluon being the first emitted TL parton (marked 3), let us discuss the case of the first emitted TL parton being a quark (also marked 3).Here, we have no choice, as can be seen in Fig. 30.Here, we have mandatory color orientation η = 1, the quark is emitted on the color side.And in the case of antiquark emission, we have η = 2.The color orientation determines which of the partons, n CC (1, k) or n CC (2, k), is the one the new emitted TL parton is connected to.In the first case (quark emission), we have for the new parton j = 3 the connection n η CJ (1, j) = n CC (η, k) = 1 (using η = 1).In the second case (antiquark emission), we have for the new parton j = 3 the connection n η CJ (1, j) = n CC (η, k) = 2 (using η = 2) always for k = 1 (projectile).
So far we treated the first TL emission of the SL parton cascade, and we explained how to use the two arrays n CJ and n CC .In a similar way we can treat all subsequent emissions.However, an emitted TL parton will, in general, develop a TL parton cascade, and we first need to deal with that.We will discuss this in the next section.

The timelike cascade
Here we consider the situation, at some stage of the SL cascade, where starting from some "current SL parton" we have an emitted TL parton with (some unique) number j, and suppose we know the color orientation η of this TL emission and the color connection array n η CJ (1, j).To simplify the notation (reduce the indices), we define for given η and given j the integer n CR (1) = n η CJ (1, j).There is also the variable n CR (2), referring to the connected parton on the other side, not yet known here, but important in case of TL emissions in the Born process, to be discussed later.
We need to reconstruct the color flow of the TL cascade and connect it to the rest of the diagram, knowing the connected parton n CR (1), which points to the previous TL parton in the color flow chain, and the color orientation η.We assume that we know at this stage the complete TL cascade, with the flavors and four-momenta of all the partons.Let us consider a concrete case as shown in Fig. 31, with  side), and with the previous TL parton in the color flow chain being n CR (1) = 4 (assuming that the indices 1,2,3 have been used already before).We also assume here that parton 4 is a final particle.We actually use a rigorous unique numbering (1,2,3,...) only for the final partons, for the intermediate partons in the plot we use symbols a, b, c . ... The first emitted TL parton a splits into b and c, and these daughter partons split again in partons, and at the end, we have final partons 5, 6, 7, 8, 9.At this stage (after realizing the TL cascade, as discussed earlier), the planar structure is already fixed, and each vertex like a → b + c has a clearly defined "first daughter" D1 (upper) and a "second daughter" D2 (lower), when considering the evolution from right to left.The fact that there is some freedom concerning the exchange of D1 and D2 is taken care of during the evolution.Here, we simply need to identify the final partons, associate a number (the next free integer), and reconstruct the color connections.This may be done as follows: 1. Define the initial parton to be the current parton.
2. Starting from the current parton, loop over all the first daughters until a final parton is reached, now being the new current parton.
3. Update connection information for the current parton itself, its parent, and the second daughter.
4. If the second daughter is not final, consider it to be the current parton, go back to 2, 5.If the second daughter is final, define the parent to be the current parton, go back to 3.
In this way, we move through the whole diagram, from top to bottom in the graph of Fig. 31.For our concrete example, we identify in this way the final partons 5, 6, 7, 8, and 9 and the color connections 4-5-6, 7-8, 9-X, where "X" in the last expression means that the corresponding parton is not yet known.

A complete example including quarks
The above procedures concern one step in the SL cascade, so we need to iterate all steps of the SL cascade on the projectile side and, of course, do the same on the target side.Finally, we need to treat the Born process, where we have two TL cascades, and for each one we employ the above procedures.Let us consider an example of a single scattering, with a diagram with SL and TL cascades, including quarks, as shown in Fig. gram, we remember that in case of a g → g + g splitting, the emitted gluon may be emitted (with equal probability) on the color side or the anticolor side (see discussion in section 4.3).A possible choice is shown in Fig. 33.Following the color flow, always starting and ending with a single arrow, we identify the following chains: 1-2-3-4, then 5-6-7, then 8-9-10-11-12, and finally 13-14.
One should note, as already mentioned several times, that even in the case of multiple scatterings, the remnants are color-white, as are the individual scatterings.In the multiple scattering procedure, where we determine the momenta of all the Pomerons, we also fix the flavors of the "Pomeron ends".Once these are known, the individual scattering can be treated independently of each other.This is why we discuss here individual scatterings.

Strings
In the previous sections, we discussed how to construct chains of partons, j 1 − j 2 − ... − j n , where the inner partons are gluons, and the two outer partons, in general, a quark and an antiquark (in any case a 3 and 3 color representation).These chains of partons are mapped (in a unique fashion) to kinky strings, where each parton corresponds to a kink, and the parton four-momentum defines the kink properties, as already done in earlier EPOS versions, for a detailed discussion see [11].Let us consider the example of Figs.32 and 34.Following the color flow, we have identified four chains, namely 1-2-3-4, 5-6-7, 8-9-10-11-12, and 13-14.Each of these chains is mapped to a kinky string, as indicated in Fig. 34.Each parton cor- responds to a kink.The inner kinks correspond to the inner gluons in the chains.The kinks carry the momenta of the corresponding partons.The chain 13-14 corresponds to a string with no inner kinks, so it is a flat (yo-yo) string.

Heavy quark issues
At each step in the SL cascade, there is the possibility of quark-antiquark production, and in the Born process as well.In the following, we discuss in particular the case of heavy flavor quarks, with the general notation Q for quarks and Q for antiquarks.Heavy flavor may be produced in different ways, as shown in Fig. 35.Starting from a gluon, a Q − Q pair may be produced in the SL cascade, as shown in Fig. 35(a), provided the virtuality is large enough.The number of allowed flavors is considered to be depending on the virtuality (variable flavor number scheme).It is also possible to create a Q − Q in the Born process, via g + g → Q + Q or q + q → Q + Q (for light flavor quarks q), as shown in Fig. 35(b), and finally Q − Q may be produced in the TL cascade, via g → Q + Q, as shown in Fig. 35(c).Let us first consider the Q − Q production in the SLC.We may have the situation as shown in Fig. 36( parton (here a Q) is emitted, and the corresponding antiparticle (here a Q) continues the SLC.But before reaching the Born process, it is emitted, and a gluon continues the SLC.The two heavy flavor partons have in general low transverse momenta.Another possibility is shown in Fig. 36(b), where a heavy flavor parton produced in the SLC "survives" till the Born process, and the latter has most likely the form Q + l → Q + l, with l being a light flavor parton.Other than the production during the SLC, heavy flavor may be produced in the Born process, via g + g → Q + Q or q + q → Q + Q (for light flavor quarks q), as shown in Fig. 36(c).Finally, heavy flavor may be produced during the timelike cascade, as shown in Fig. 37, initiated either from a TL parton in the SL cas- of partons: 1-2-3-4-5, 6-7-8, and 9-10-11.The initial TL partons (the horizontal blue lines with arrows) or most likely quarks and antiquarks (in any case 3 and 3 color representations).Let us assume that 3 is a quark, and 6 an antiquark (light flavor, both), then the two chains containing heavy flavor are of the form Q − g − q and q − g − Q, in both cases, the heavy flavor partons are "end partons" in the chains.These chains of partons are finally mapped (in a unique fashion) to kinky strings, where each parton corresponds to a kink, as shown in Fig. 39.
The general mapping procedure (chains of partons to kinky strings) as well as the string decay procedures are described in detail in [11].

Summary
We recently introduced new concepts [6], implemented in EPOS4, which allows to consistently accommodate fac- torization and saturation in high-energy proton-proton and nucleus-nucleus collisions, in a rigorous parallel scattering framework.EPOS4 has a modular structure: the multiple scattering contributions to the total cross section in pp or AA scattering are expressed in terms of (products of) cut Pomeron expressions G (each one representing a single scattering), and the latter ones are related to the "real QCD expressions" G QCD via some fundamental (new) equation.In other words, G QCD is the fundamental building block of the multiple scattering framework of EPOS4.In this paper, we provided detailed information about the precise structure of G QCD and its calculation, based on perturbative QCD, with special care concerning heavy flavor.We discussed the implementation (for the first time in the EPOS framework) of the "backward parton evolution method", which allows a much better control of the hard processes and comparisons with the "standard" pQCD calculations based on factorization.But compared to those, there are many technical subtleties in our approach, due to the fact that we do the parton evolutions individually for each of the scatterings in a multiple scattering configuration, and this has been discussed in great detail.It is actually the occurrence of singularities for cross section calculations as well as for parton generation which needed to be taken care of.We also discussed the way to transform the partonic multiple scattering structure into strings.
KW thanks Sergej Ostapchenko for many helpful discussions concerning the technical aspects of this paper.KW also thanks Tanguy Pierog for useful contributions related to the "deformation function".BG acknowledges support from ANID PIA/APOYO AFB220004.

APPENDIX A Evolution function A.1 Differential equation for evolution function
It will be useful to have a the differential equation, corresponding to Eq. ( 46).We may write then we compute the derivative ∂ ∂t and multiply by t, which gives

A.2 Evolution function theorem
For the proof of Eq. (56) we start from Eq. (169).We first define a matrix E N QCD (t a , t) with the components and we define a vector g i N (t a , t) to be the i th column of the matrix E N QCD (t a , t), divided by ∆ k (t).We further define a matrix γ N (t) with elements The Mellin transform of Eq. ( 169) is then which has the solution (the exponential of a matrix M is defined by its power series ∑ ∞ i=0 1 i! M i ).Using E ik QCD (x, t a , t a ) = δ(1 − x)δ ik , we have g i N (t a , t a ) = e i /∆ k (t) (with e i being a vector with the i th component unity and all others zero), and so we get for the matrix E N QCD : Considering t = t b , and splitting the integration interval, we get with X = 1 2 [A, B] + ..., with "..." representing higher commutators of X and Y.In leading log accuracy, we can neglect X, since it amounts to higher orders in α s , and so we get or, with the corresponding matrix elements, The inverse Mellin transformation gives QED.

B Parton-parton cross section formulas
Here we will show how to explicitly do some of the integrations, in order to compute cross sections, starting from Eqs. (45,44).

B.1 Rapidity integral
Here we simplify the expression to get integrated cross sections as rapidity integral.and δ mn accounting for properly counting identical partons.

B.2 LC momentum fraction integral
Here we want to show that σ ij hard in Eq. ( 225) can be written as an integral over light-cone momentum fractions.

B.3 Single jet parton-parton cross section
Here we show that d 3 σ ij hard /dy d 2 p t can be written as a single integral dx..., where y and p t refer to one outgoing parton (jet).
Again starting from Eq. ( 44), one may integrate over p 4 , to get where y and p t refer to particle 3.In the BornCMS, we have (with , not yet imposing energy conservation) and so 1 with {...}being the form the squared matrix elements are usually tabulated, with α s = g 2 /4π, and where the integration limits a and b are obtained from x 2 ≤ 1 and x 1 ≤ 1.This formula is useful to compute (without Monte Carlo) transverse momentum spectra of produced partons for tests and comparisons.
In the EPOS framework, we have evolution functions with an initial condition δ(1 − x), which forces us to distinguish four cases, namely both-sided, onesided(lower), one-sided(upper), or no emissions.In the case of both-sided emissions, we use Eq. ( 256), with Ẽab QCD instead of E ab QCD .In case of no emissions on either side, we have E ab QCD (x i , Q 2 i , µ 2 F ) = ∆ a QCD (Q 2 i , µ 2 F )δ(1 − x i )δ ab on both sides, so dx 1 dx 2 ... is trivial, and using δ(s + t + u) = δ( f (y)) = δ (y − y 0 ) f ′ (y 0 ) +δ (y + y 0 ) f ′ (−y 0 ) , with We finally mention that the relation (from s + t + u = 0) not only allows to express x 2 in terms of x 1 and y or x 1 in terms of x 2 and y (for given a = p 2 t /s lad ), as but also the other way round y in terms of x 1 and x 2 , as or useful formulas for later applications.

C Born scattering kinematics: Explicit formulas
Here we will provide explicit formulas for different kinematic quantities and relations related to the Born process, as being used in the numerical procedures.We consider two partons entering the Born process (let us note them 1 and 2), producing two partons (let us note them 3 and 4), so altogether we have the elementary process 1 + 2 → 3 + 4. All kinematic formulas depend on the four masses, m 1 , m 2 , m 3 , and m 4 .In order to present explicit formulas for the different cases, we will use the notation case m 1 m 2 m 3 m 4 , (271) so for example "case 0000" represents the case of massless incoming and outgoing partons, like g + g → g + g or u + ū → g + g.

Figure 7 :
Figure 7: The integration path in the k − plane.
8, with the corresponding flavors k, l, m

Figure 8 :
Figure 8:The differential parton-parton cross section

2 .
cl → cl : scattering of charmed quarks or antiquarks with light partons (case m0m0 with m = m c ) 3. bl → bl : scattering of bottom quarks or antiquarks with light partons (case m0m0 with m = m b ) 4. l l → c c : annihilation of light pairs into charmed pairs (case 00mm with m = m c ) 5. l l → b b : annihilation of light pairs into bottom pairs (case 00mm with m = m c ) 6. c c → l l : annihilation of charmed pairs into light pairs (case mm00 with m = m c ) 7. b b → l l : annihilation of bottom pairs into light pairs (case mm00 with m = m b ) 8. c c → c c : scattering of charmed pairs into charmed pairs (case mm m m with m = m = m c ) 9. c c → b b : scattering of charmed pairs into bottom pairs (case mm m m with m = m c , m = m b ) 10. b b → c c : annihilation of bottom pairs into charm pairs (case mm m m with m = m b , m = m c ) 11. b b → b b : annihilation of bottom pairs into bottom pairs (case mm m m with m = m = m b ) 12. cb → cb : scattering of charm and bottom partons (case m mm m with m = m c , m = m b )

Figure 12 :
Figure 12:Transverse momentum distributions of primary partons for Pomerons of type "val-val", "val-sea", "sea-val", and "sea-sea".We compare the Monte Carlo results (points) to the corresponding curves obtained with the help of explicit cross section formulas (see text).

Figure 13 :
Figure 13: 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 sub-plot .The red curves refer to EPOS PDFs, the green ones to CTEQ PDFs, and the black points are data from ZEUS and H1.

Figure 14 :
Figure 14: Parton yield dn/dp t dy for pp at 13 TeV.We 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 15 :
Figure 15: Transverse momentum distribution of charm quarks (c and c), based on EPOS PDFs and CTEQ PDFs.In both cases, we show results including flavor excitation (incl FlavEx) or not (w/o FlavEx).The "jumps" mark the activation of charm flavor excitation for µ F > 2 m c .

Figure 16 :Figure 17 :
Figure 16: Transverse momentum distribution of charm quarks (c and c), based on EPOS PDFs and CTEQ PDFs.In both cases, we show results including flavor excitation (incl FlavEx) or not (w/o FlavEx).The black circles represent FONLL results, and the gray squares NLO calculations.

Figure 18 :
Figure 18: Charmed partons and hadrons in pp collisions at 7 TeV.We show EPOS4 simulations, compared to FONLL results and ALICE data.

6 Figure 19 :
Figure 19: Partonic configuration of two colliding nuclei A and B, each one composed of two nucleons, with three scatterings (from three cut Pomerons).Dark blue lines mark active quarks, red dashed lines active antiquarks, and light blue thick lines projectile and target remnants.One of the target nucleons is just a spectator.

Figure 20 :Figure 21 :
Figure 20: Configuration colorwise equivalent to the one of Fig. 19.The outgoing antiquarks are drawn as incoming quarks (arrows towards vertices).

Figure 22 :
Figure 22: Simplified color flow diagram for the three scatterings of Fig. 21.The dot indicates the Born process.

Figure 23 :
Figure 23: timelike cascade, with Feynman diagram (a) and the corresponding color flow diagram (b).

Figure 26 :Figure 27 :
Figure 26: Defining color and anticolor side for the SL evolution.The red arrows represent the directions of the evolutions, starting from the projectile or the target side.

ηFigure 29 :
Figure 29: The two possibilities of color flow for a gluon emission: (a) emission on the color side, corresponding to η = 1 or (b) the anticolor side, corresponding to η = 2.

Figure 30 :
Figure 30: The two possibilities of color flow for the first emission.The l.h.s. of the figure corresponds to the color orientation η = 1, the r.h.s. to η = 2

Figure 32 :
Figure 32: A single parton-parton scattering diagram, with SL and TL cascades, including quarks.

Figure 33 :
Figure 33: A possible color flow diagram corresponding to Fig. 32.

Figure 35 :
Figure 35: Different possibilities to create heavy flavor, (a) in the spacelike cascade (SLC), (b) in the Born process, (c) in the timelike cascade (TLC).

Figure 36 :
Figure 36: Heavy flavor production (a,b) in the SL cascade and (c) in the Born process.The magenta point indicated the Born process.

Figure 37 :Figure 38 :
Figure 37: Heavy flavor production in the TL cascade.cade[Fig.37(a)], or from an outgoing parton of the Born process (Fig.37(b)).In the first case, the transverse momenta are in general small.The next step will be, for a given Feynman diagram, to construct the color flow diagram, as discussed in the previous sections.Let us take the graph of Fig.37(b), i.e., heavy flavor production during the TLC of an outgoing Born parton.As usual, the gluons are emitted to either side with equal probability, so a possible color flow diagram is the one shown in Fig.38.We identify three chains
) which has the structure e A+B with two matrices A and B. The Baker-Campbell-Hausdorff formula states e A e B = e A+B+X ,

Table 1 :
The symbols G, T, T, and σ.