# Parallel scattering, saturation, and generalized AGK theorem in the EPOS4 framework, with applications for heavy ion collisions at 5.02 ATeV and at 200 AGeV

Klaus Werner

► To cite this version:

Klaus Werner. Parallel scattering, saturation, and generalized AGK theorem in the EPOS4 framework, with applications for heavy ion collisions at 5.02 ATeV and at 200 AGeV. Physical Review C, 2024, 109 (3), pp.034918. 10.1103/PhysRevC.109.034918. hal-04258910

HAL Id: hal-04258910

https://hal.science/hal-04258910v1

Submitted on 16 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Parallel scattering, saturation, and generalized Abramovskii-Gribov-Kancheli (AGK) theorem in the EPOS4 framework, with applications for heavy-ion collisions at $\sqrt{s_{NN}}$ of 5.02 TeV and 200 GeV

K. Werner

SUBATECH, Nantes University, IN2P3/CNRS, IMT Atlantique, 44300 Nantes, France

<figure><img src="image_2.png" /></figure>

(Received 19 October 2023; revised 12 December 2023; accepted 23 January 2024; published 27 March 2024)

Ultrarelativistic heavy-ion collisions will first realize many nucleon-nucleon scatterings, happening instantaneously and therefore necessarily in parallel, due to the short collision time. An appropriate quantum mechanical tool to treat that problem is S-matrix theory, and it has been known for a long time how to derive a simple geometric probabilistic picture, still widely used, and here the Abramovskii-Gribov-Kancheli (AGK) theorem plays a crucial role. All this is done in a scenario where energy conservation is not taken care of, but this is needed, in particular for Monte Carlo simulations. When introducing energy-momentum sharing properly, the AGK theorem does not apply anymore, nor do simple geometric concepts such as binary scaling. I will discuss this (very serious) problem, and how it can be solved, in the EPOS4 framework. When connecting the multiple-Pomeron approach (for parallel scatterings) and perturbative QCD, one is actually forced to implement in a very particular way saturation scales in order to get an approach free of contradictions. One recovers a generalized AGK theorem (gAGK), valid at large $p_t$ (larger than the relevant saturation scales). I discuss how gAGK is related to factorization (in proton-proton scatterings) and binary scaling (in heavy-ion collisions). I will show some applications, using this new approach as an initial condition for hydrodynamical evolutions, for heavy-ion collisions at $\sqrt{s_{NN}}$ of 5.02 TeV and 200 GeV, to get some idea about the energy dependence.

DOI: 10.1103/PhysRevC.109.034918

## I. INTRODUCTION

It is of fundamental importance to realize that multiple nucleon-nucleon scatterings in heavy-ion collisions must happen in parallel, and not sequentially, based on very elementary considerations concerning timescales. The appropriate tool to take this into account has been known for a long time: an S-matrix approach, referred to as Gribov-Regge (GR) theory [1–3]. A key property is the Abramovskii-Gribov-Kancheli (AGK) theorem, also referred to as “AGK cancellations” [4], which allows to make the link between the multiple scattering approach and geometric properties such as binary scaling.

However, for a realistic scenario—in particular as a basis for event-by-event Monte Carlo procedures—it is mandatory to include energy-momentum conservation, which is not considered in the approach mentioned above. This sounds trivial, but actually implementing it is a highly complex operation, in a scenario where all NN collisions happen in parallel, avoiding any ordering of these collisions, which would make no sense. If a projectile nucleon interacts with several target nucleons, there is nothing like a first and subsequent collisions; they are all equal. Despite many technical difficulties, an unbiased energy-momentum sharing can be dealt with [5],

but unfortunately it destroys the “AGK cancellations” and, as a consequence, all the nice and simple geometric properties which follow, in particular binary scaling.

The aim of this paper is to first discuss and finally understand the serious problems which arise (unavoidably) in a scenario with appropriate energy-momentum sharing and at the same time unbiased parallel scattering, which at the end leads to an amazingly simple solution of all problems based on a particular implementation of saturation.

The present paper is the fourth in a series of publications [6–8] about the EPOS4 framework. One distinguishes between primary and secondary interactions. The former refers to the multiple parallel scatterings happening (at high energies) instantaneously and which result in complex configurations composed of many strings, whereas the latter refers to subsequent interactions of the string decay products, which amounts to first a core-corona separation based on the string segments, and then the fluid formation, evolution, and decay of the core part. It should also be mentioned that the discussion of primary interactions covers two topics, namely, the parallel scattering formalism, developed in terms of abstract objects for single scatterings called “Pomerons”, and the internal structure of the Pomerons, showing how the Pomeron is related to parton-parton cross sections expressed in terms of QCD diagrams. Reference [6] represents an overview of EPOS4; in Ref. [7], the internal Pomeron structure is treated, with very detailed discussions of the parton-parton scatterings based on perturbative QCD (which makes the link between the multiple scattering formalism and QCD). Reference [8] focuses on the secondary interactions, in particular on the

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP³.</footnote>
core-corona procedure and microcanonical hadronization of the fluid.

The present paper refers to the primary interactions, and it focuses on a very detailed and rigorous treatment of the multiple scattering aspect based on Pomerons, referring when necessary to the internal structure discussed in Ref. [7]. The aim is the presentation of a consistent multiple scattering formalism, which takes into account energy-momentum sharing (which is mandatory) and rigorous unbiased parallel scattering (which is mandatory as well), and which finally solves the problem of violating AGK (and simple geometrical pictures). The key element is an appropriate implementation of dynamical saturation scales.

In the following sections, I will first discuss the goal and the philosophy of the formalism developed in this paper, with precision concerning the assumptions and the limitations. Then I will discuss timescales and the corresponding applicability of parallel scattering as a function of the collision energy. Since “AGK cancellations” are the key element in the parallel scattering scheme, I will briefly review this concept in the classical approach [1–4], before generalizing towards a scenario including energy-momentum conservation. I will discuss all the difficulties that show up, and eventually the solution, based on saturation scales, leading to the generalized AGK theorem (gAGK), valid at large $p_t$, which allows to restore the above-mentioned geometrical properties like binary scaling. Finally, I show results for PbPb at $\sqrt{s_{NN}}$ of 5.02 TeV and AuAu at 200 GeV.

## II. THE PHILOSOPHY, BASIC ASSUMPTIONS, AND OBJECTIVES

Let me make some remarks about the philosophy behind the present work, the main goals, and the assumptions employed. The aim is not to construct an event generator, but to provide a theoretical framework for “multiple scatterings”, where one considers multiple partonic scatterings in proton-proton and multiple nucleon-nucleon scatterings in nucleus-nucleus collisions.

A crucial point is the fact that one is aiming at a formalism which allows access to event classes. Most theoretical work based on QCD is devoted to inclusive cross sections, and here factorization plays a dominant role. The approach in this paper does not contribute to these admittedly very important achievements, being very useful to study rare processes. But one of the highlights at the LHC and a major activity for (at least) the next decade are studies concerning particle production as a function of the “event activity” (measured in terms of multiplicity or transverse energy). Models are requested to compute particle yields or two-particle correlations, but with a very important additional condition: for given event classes. And this makes life complicated. Minimum bias results based on factorization are not able to provide answers—or provide wrong answers, citing, for example, the (experimentally well-known) increase of the mean transverse momentum with multiplicity in proton-proton scattering.

From the huge amount of experimental data and phenomenological attempts to interpret these, there is strong evidence that event activity is very strongly correlated

with multiple parton interactions (MPIs). Reference [9] discusses the role of multiple parton scatterings at the LHC: one distinguishes between “hard MPI” and “soft MPI”, where the former treats the scattering of two hard partons, based on multiple parton distribution functions, generalizing the factorization approach. But this does not address the questions related to studies of properties for given event classes, defined via “event activity”. This is discussed in Ref. [9] in the part “soft MPI”, and here the situation is not satisfactory; the “phenomenological approaches” have often no solid theoretical basis, and here I am not even talking about a first-principles basis.

The discussion of “event classes” is clearly also very relevant for scattering of nuclei with mass numbers A and B (A+B scattering). Nobody considers minimum bias results; everything is presented again for event classes, again defined via event activity. Also here many things are known from phenomenological studies, like the (measurable) event activity being very strongly correlated with the (nonmeasurable) impact parameter and also with the (nonmeasurable) number of nucleons being involved. There are also formulas concerning factorization in nucleus-nucleus scattering, but also here one needs a framework where one has access to information for given event classes, and one needs to identify the theoretical counterpart of these event activities.

It is completely out of reach to construct such a framework (with access to event classes) from first principles (QCD), and this is not at all what will be attempted in this paper. What can be done (and this will be developed in this paper) is to construct a scheme in the spirit of S-matrix theory, which was the main theoretical tool to treat scatterings before the QCD era. One assumes certain fundamental properties, like Lorentz invariance and analyticity of the T-matrix, although it is not possible to derive the particular T-matrix from first principles. But at least the assumptions are very clear and transparent. Here the same strategy is adopted. One starts with a certain assumption about the form of the T-matrix, which is only based on a qualitative (or phenomenological) understanding of experimental data. It cannot be proven from first principles, but it is well defined, and even very simple, and one respects perfectly the framework of quantum mechanics (this is worth noticing, because many “models” for heavy-ion scatterings are based on purely classical considerations).

Before discussing more in detail the assumptions of the present approach, some historical remarks are in order. In his famous lecture [10], Glauber discusses the scattering of high-energy protons with nuclei. There are several assumptions used, like a factorization of the wave function into a plane wave times a slowly (in space) changing function. The assumptions are clear and plausible, and one gets finally a simple result with a nice interpretation, of a particle moving on a straight line, accumulating phases (the so-called eikonal phase is an integral over the potential along the longitudinal axis). Another early (pre-QCD) approach is GR theory [1–3]. Also here very strong assumptions are made about the structure of the T-matrix, and the form of the sub-T-matrices. But also here, the assumptions cannot be proven from some fundamental theory, but they are plausible, clear, and transparent. Interesting enough, the Glauber and the GR approaches
<figure><img src="image_3.png" /><figcaption>FIG. 1. Factorization in pp scattering.</figcaption></figure>

(the latter being equal to the EPOS approach without energy sharing and without saturation) give similar results.

The scheme to be developed in this paper has to be understood in the spirit of these two (historical) approaches. It is not attempted to “derive” the formalism from a fundamental theory (QCD), but to construct a model, in the framework of quantum mechanics (with the possibility to treat elastic and inelastic scattering, related by the optical theorem), and which is finally compatible with QCD. There are assumptions, which are purely based on a phenomenological (and partly qualitative) understanding of experimental data. But, the assumptions are very clear and transparent (for this purpose I will try to summarize all the important points in the following).

But let me first have some qualitative discussion. The usual factorization picture where one attempts to separate “long range (soft)” and “short range (hard)” parts of the interaction may be graphically presented as shown in Fig. 1. In this plot and all the following ones, I show for simplicity only gluons; in reality all kinds of partons are present. The two light blue thick lines represent the projectile and the target protons. The proton structure and the so-called spacelike parton cascade are taken care of by using parton distribution functions (PDFs) $f$, which allow writing the jet cross section as a convolution of these PDFs and an elementary QCD cross section for the Born process in the middle. This approach provides excellent results concerning inclusive particle production. From studies of particle production for given event classes, for example, charm production versus charged particle multiplicity, it seems that the “real events” look more like that sketched in Fig. 2, where one has not only several (three in the figure) scatterings, but even several parallel “objects” (multiple scattering). The fact of having several of these objects not only increases (for example) the rate of heavy flavor production, it also increases the charged particle multiplicity, and each

<figure><img src="image_4.png" /><figcaption>FIG. 2. Multiple scattering.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 3. Modular structure of multiple scattering.</figcaption></figure>

object produces particles over the whole range in rapidity. This is why in the sketch of Fig. 2 the three objects are similar to the one shown in Fig. 1, with parton emissions, but per object. Some “evolution function” $E$ is indicated which should obey the same evolution equations as the $f$ in Fig. 1. The picture sketched in Fig. 2 suggests some modular structure, as indicated in Fig. 3, where one first simply considers “objects” indicated by the red boxes, where one might separate the multiple scattering aspect (expressed in terms of the objects) and the content of the objects, representing parton-parton scattering based on QCD diagrams.

The main objective of this work (with this paper being an important piece of) is to investigate (and understand) how a picture as sketched in Fig. 2 breaks down to the one sketched in Fig. 1 when it comes to inclusive particle production (in minimum bias scattering). Even more drastic, one has to show how in $A+B$ scattering a picture of $A \times B$ possible interactions, each one of the form of Fig. 2, breaks down as well to the one sketched in Fig. 1, simply multiplied by a factor $A \times B$ (which seems to be an experimental fact). These are crucial questions; any multiple scattering approach, which claims to describe observables “per event class”, must do this exercise.

The ideas sketched in the previous paragraph will be realized based on a model, generalizing the GR approach [to include energy-momentum sharing (GR$^+$), and to make it QCD compatible]. This is compatible with the picture sketched in Fig. 3 and the experimental data which suggest such a picture. There are serious problems coming up, for which solutions can finally be proposed. But they might be questioned. This is why I will present in the following all the important assumptions and hypotheses (marked H1a, H1b, etc.).

(H1a) The most important assumption says that the T-matrix (first for pp elastic scattering) has a modular structure: it is given in terms of products of sub-T-matrices $T_{Pom}$ referring to elementary scatterings between “proton constituents” by the exchange of a “Pomeron” (whatever this may be). I consider first the case without energy-momentum sharing. This is precisely what is done in Refs. [1–3]. One gets [see Eq. (32)]

$$iT = \sum_{n=0}^{\infty} \frac{1}{n!} \{iT_{Pom}\}^n. \quad (1)$$
(H1b) A fundamental assumption is the complete separation of the multiple scattering structure, expressed in terms of the sub-T-matrices $T_{Pom}$ and the underlying microscopic theory (QCD), which allows to first develop the multiple scattering theory based on $T_{Pom}$, and postpone the precise ‘internal structure’ of $T_{Pom}$ in terms of QCD diagrams.

The optical theorem (in impact parameter representation; see Sec. IV) allows to compute the inelastic cross section, as $\sigma_{in} = \int d^2b$ cut $T(b)$, with $b$ being the impact parameter. Cutting a complex diagram amounts to summing all possible cuts.

(H1c) Here one assumes that a Pomeron is cut completely, or not at all. Summing over all cuts then means summing over all possibilities of cut and uncut Pomeron. Assuming in addition $T_{Pom}$ to be purely imaginary, and defining the “cut Pomeron” $G = 2 \operatorname{Im} T_{Pom}$, one gets

$$ \sigma_{in} = \int d^2b \sum_{n=1}^{\infty} \frac{1}{n!} \sum_{m=0}^{n} \binom{n}{m} G^m(-G)^{n-m}. \quad (2) $$

In principle the T-matrices depend on the collision energy (expressed in terms of $s$) but for simplicity I will not write this dependence explicitly here (but it will be done properly later). In order to generalize towards $A+B$ scattering, one defines $b_i^A$ and $b_j^B$ being the transverse coordinates of the nucleons in the nuclei A and B, and an integration

$$ \int db_{AB} = \int d^2b \int \prod_{i=1}^{A} d^2b_i^A T_A(b_i^A) \int \prod_{j=1}^{B} d^2b_j^B T_B(b_j^B), \quad (3) $$

with the nuclear thickness functions $T_A$ and $T_B$ [see Eq. (41)] representing the nuclear geometry.

(H2a) It is assumed that the T-matrix expression for $A+B$ scattering is given as an integral $\int db_{AB} \dots$.

(H2b) The integrand is assumed to be a product of $A \times B$ nucleon-nucleon expressions as Eq. (1), where the impact parameter argument of the T-matrices corresponding to the $k$th nucleon-nucleon pair is $b_k = |b + b_{\pi(k)}^A - b_{\tau(k)}^B|$. The functions $\pi(k)$ and $\tau(k)$ refer to the nucleons associated to the pair $k$. Defining $G_k = G(b_k)$, one gets

$$ \sigma_{in}^{AB} = \int db_{AB} \sum_{m_1 l_1} \dots \sum_{m_{AB} l_{AB}} \prod_{k=1}^{AB} \frac{1}{m_k! l_k!} (G_k)^{m_k} (-G_k)^{l_k}. \quad (4) $$

What is discussed so far is essentially the GR approach, which does not address energy-momentum sharing. To incorporate that, one introduces the light-cone momentum fractions $x^\pm$ with respect to the initial light-cone momenta of the nucleons, with $0 \le x^+, x^- \le 1$. One defines $x_{k\mu}^\pm$ to be the light-cone momentum fractions of the external legs of the $\mu$th Pomeron of pair $k$, all of them connected to projectile nucleon $i = \pi(k)$ and target nucleon $j = \tau(k)$. The light-cone momenta of the

projectile and target remnants are named $x_{remn,i}^+$ and $x_{remn,j}^-$.

$$ x_{remn,i}^{+} = 1 - \sum_{k=1}^{AB} \sum_{\substack{\mu=1 \\ \pi(k)=i}}^{n_k} x_{k\mu}^{+}, \quad x_{remn,j}^{-} = 1 - \sum_{k=1}^{AB} \sum_{\substack{\mu=1 \\ \tau(k)=j}}^{n_k} x_{k\mu}^{-}, \quad (5) $$

which means that the initial values ($x^\pm = 1$) are shared among Pomeron and remnants.

(H3a) One assumes that energy-momentum conservation can be incorporated by simply using $x_{k\mu}^{\pm}$ as arguments of the sub-T-matrices and

(H3b) by adding Pomeron-nucleon vertices $V(x_{remn,i}^+)$ and $V(x_{remn,j}^-)$ representing the coupling of the Pomeron to the projectile and target remnants, which ensures energy-momentum conservation.

As a consequence, defining $G_{k\mu} = G(x_{k\mu}^+, x_{k\mu}^-, b_k)$, and an integration $\int dX_{AB}$ over all light-cone momentum fractions [see Eq. (53)], one gets (see Sec. VI C) $\sigma_{in}^{AB} = \int db_{AB} G^{AB}(\{b_{AB}\})$, with $\{b_{AB}\}$ being the multidimensional variable $\{b, \{b_i^A\}, \{b_j^B\}\}$, and with

$$ G^{AB}(\{b_{AB}\}) = \sum_{n_1=0}^{\infty} \dots \sum_{n_{AB}=0}^{\infty} \sum_{m_1 \le n_1} \dots \sum_{m_{AB} \le n_{AB}} \int dX_{AB} \times \prod_{k=1}^{AB} \frac{1}{n_k!} \binom{n_k}{m_k} \prod_{\mu=1}^{m_k} G_{k\mu} \prod_{\mu=m_k+1}^{n_k} -G_{k\mu} \times \prod_{i=1}^{A} V(x_{remn,i}^{+}) \prod_{j=1}^{B} V(x_{remn,j}^{-}), \quad (6) $$

with at least one $n_k$ being nonzero. This extension of the GR approach (by including energy-momentum sharing) will be referred to as GR$^+$.

One may define an expression $\bar{G}^{AB}(\{b_{AB}\})$ as in Eq. (6), but for all $n_k = 0$, which represents the case excluded in the summation of Eq. (6). One can prove [see Eq. (70)] the following relation:

$$ G^{AB}(\{b_{AB}\}) + \bar{G}^{AB}(\{b_{AB}\}) = 1. \quad (7) $$

In the case of $A = B = 1$ (pp scattering) one then gets

$$ \sigma_{in}^{pp} = \int db \{1 - \bar{G}^{11}(b)\}. \quad (8) $$

It is very tempting to interpret the expression $\{1 - \bar{G}^{11}(b)\}$ as the probability of an interaction at given impact parameter $b$. In that case, $\bar{G}^{11}(b)$ must be non-negative; otherwise one exceeds the “black disk limit”. It is also very plausible that this condition should hold for $A > 1$ and/or $B > 1$. So one requires the following:

(H4) Whatever approximation might be employed, it is mandatory that $\bar{G}^{AB}(\{b_{AB}\})$ is non-negative, for any value of A and B. This is a fundamental expectation which will guide the following hypotheses. It can be proven in the case without energy sharing, and is here expected to be true as well.
This is a crucial requirement; it has enormous consequences, as will be seen later. In order to investigate $G^{AB}$ and $\tilde{G}^{AB}$ (and in particular the sign of the latter), one has to be more specific. At this point the $G$ functions are just abstract objects, but eventually they should correspond to “real” expressions based on QCD diagrams representing parton-parton scattering. Studying such expressions, it could be shown [5] that one may obtain an almost perfect fit of the numerically computed functions $G_{QCD}$, with a parametrization being a sum of expressions of the form $\alpha_N(x^+x^-)^{\beta_N}$. This parametric form has been inspired by the asymptotic expressions for T-matrices (see Appendix A). So one postulates the following:

(H5a) With $\alpha_N$ and $\beta_N$ being coefficients depending on $b$ in terms of a few parameters, the functions $G(x^+, x^-, b)$ can be written as (see Appendix B)

$$G_{\text{QCDpar}}(x^+, x^-, b) = \sum_{N=1}^{4} \alpha_N (x^+ x^-)^{\beta_N}. \quad (9)$$

(H5b) Furthermore, the vertices can be parametrized as $V(x) = x^{\alpha_{\text{remn}}}$, again motivated by the asymptotic expressions for T-matrices.

These postulates are crucial; they allow one to do the integrals, which could not be done numerically. Defining $\{m_k\}$ to be the set of the $m_k$ variables, and $\{x_{k\mu}^{\pm}\}$ the set of the $x_{k\mu}^{\pm}$ variables, one finds (see Sec. VI C)

$$G^{AB}(\{b_{AB}\}) = \sum_{\{m_k\}} \int \prod_{k=1}^{AB} \prod_{\mu=1}^{m_k} dx_{k\mu}^{+} dx_{k\mu}^{-} P(\{m_k\}, \{x_{k\mu}^{\pm}\}), \quad (10)$$

and $\tilde{G}^{AB} = P(\{m_k\} = 0)$, with

$$P(\{m_k\}, \{x_{k\mu}^{\pm}\}) = \prod_{k=1}^{AB} \left[ \frac{1}{m_k!} \prod_{\mu=1}^{m_k} G_{k\mu} \right] \times W_{AB}, \quad (11)$$

where the expression $W_{AB}$ is a function of the variables

$$x_i^+ = 1 - \sum_{\substack{k=1 \\ \pi(k)=1}}^{AB} \sum_{\mu=1}^{m_k} x_{k\mu}^+, \quad x_j^- = 1 - \sum_{\substack{k=1 \\ \pi(k)=j}}^{AB} \sum_{\mu=1}^{m_k} x_{k\mu}^-, \quad (12)$$

with $1 \le i \le A$ and $1 \le j \le B$. One finds [see Eq. (88)]

$$W_{AB}(\{x_i^+\}, \{x_j^-\}) = \prod_{i=1}^{A} (x_i^+)^{\alpha_{\text{remn}}} \prod_{j=1}^{B} (x_j^-)^{\alpha_{\text{remn}}} \sum_{\{r_{Nk}\}} \left\{ \prod_{k=1}^{AB} \prod_{N=1}^{4} \frac{(-\alpha_N)^{r_{Nk}}}{r_{Nk}!} \right. \\ \times \prod_{i=1}^{A} \left[ \prod_{k=1}^{AB} \prod_{N=1}^{4} (\Gamma(\tilde{\beta}_N)(x_i^+)^{\tilde{\beta}_N})^{r_{Nk}} g\left(\sum_{k=1}^{AB} \sum_{N=1}^{4} r_{Nk} \tilde{\beta}_N\right) \right] \\ \left. \times \prod_{j=1}^{B} \left[ \prod_{k=1}^{AB} \prod_{N=1}^{4} (\Gamma(\tilde{\beta}_N)(x_j^-)^{\tilde{\beta}_N})^{r_{Nk}} g\left(\sum_{k=1}^{AB} \sum_{N=1}^{4} r_{Nk} \tilde{\beta}_N\right) \right] \right\}, \quad (13)$$

where $\sum_{\{r_{Nk}\}}$ means summing all the indices $r_{Nk}$, with $1 \le N \le 4$ and with $1 \le k \le AB$, from zero to infinity, where

<figure><img src="image_6.png" /><figcaption>FIG. 4. Nonlinear effects: ladders, which evolve first independently and in parallel, finally fuse.</figcaption></figure>

$r_{Nk}$ refers to the number of uncut “Pomerons of type N” of nucleon-nucleon pair $k$. It is useful for the discussion to consider “Pomeron types” N, although they are not physical objects, just coming from the parametrization in Eq. (84). I use $\tilde{\beta}_N = \beta_N + 1$, and a function $g$ defined as

$$g(z) = \frac{\Gamma(1 + \alpha_{\text{remn}})}{\Gamma(1 + \alpha_{\text{remn}} + z)}. \quad (14)$$

All the integrals could be done (see Appendix B). One gets $\tilde{G}^{AB} = W_{AB}(\{x_i^+ = 1\}, \{x_j^- = 1\})$, and since $\tilde{G}^{AB}(\{b_{AB}\})$ has to be non-negative, one has the requirement

$$W_{AB}(\{x_i^+ = 1\}, \{x_j^- = 1\}) \ge 0. \quad (15)$$

At least for $A = B = 1$ the numerical calculation of $W_{AB}$ can be done. Although not written explicitly, this quantity depends on $b$. For a large impact parameter $b$, one finds indeed $W_{11}(1, 1) > 0$, but below $b = 1$ fm, one gets negative values, up to roughly $-0.4$ at $b = 0$ (see Sec. VI C).

(i) So here the requirement of Eq. (15) is violated.

(ii) The result one finds is “unreasonable”; it contradicts common sense. Let me refer to this as the “sign problem”.

What does this mean? It seems that something is missing. Here one needs again some input from a qualitative understanding of high-energy scattering. It is known that with increasing energy, partons with very small momentum fractions $x \ll 1$ become increasingly important, the parton density becomes very large, and therefore the linear Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution scheme is not valid anymore. Nonlinear evolution takes over, considering explicitly gluon-gluon fusion. These phenomena are known as “saturation” [11–24]. The diagrams for each scattering actually look more like the one shown in Fig. 4. One expect “nonlinear effects”, which means that two ladders, which evolve first independently and in parallel, finally fuse. And only after that, the (linear) DGLAP evolution might be realized. So ideally one should generalize the framework; in addition to simple exchanges of Pomerons, one should include “Pomeron fusions”, i.e., triple Pomeron graphs, in an iterative fashion up to infinite order. This “more realistic scenario” should solve the sign problem. But to do that in a framework
<figure><img src="image_7.png" /><figcaption>FIG. 5. Nonlinear effects (inside the red ellipses) are “summarized” in the form of saturation scales.</figcaption></figure>

with energy-momentum sharing (considered to be crucial) seems impossible.

Here comes another important fact: nonlinear effects lead to strong destructive interference, which may be summarized in terms of a saturation scale [12,13]. So one might think of treating these saturation phenomena not explicitly, but by introducing saturation scales as the lower limit of the virtualities for the DGLAP evolutions, as sketched in Fig. 5. So the diagrams inside the red ellipses are replaced by two scales $Q^2_{sat,proj}$ and $Q^2_{sat,targ}$. Based on these considerations, one might think of keeping the modular structure as presented above, but introduce saturation via saturation scales, as sketched in Fig. 6, where one still uses objects called Pomerons as basic modules, but the intrinsic structure of these modules takes into account saturation scales.

So far I have discussed the idea; in the following I will discuss how to develop the corresponding formalism. But before doing so, there is another important element to be discussed, and this is related to binary scaling, which states that the inclusive cross section concerning the production of hard probes (high $p_t$ or heavy flavor particles) in A+B scattering is AB times the corresponding $pp$ cross section. This is shown experimentally, and anyway it is difficult to imagine how it could be different. Based on the formalism discussed so far (the one which leads to the sign problem), one can prove that binary scaling holds [see Eq. (65)]. One can as well prove the AGK theorem, saying that the inclusive A+B cross section is AB times the single Pomeron contribution, being an

<figure><img src="image_8.png" /><figcaption>FIG. 6. Modular structure of multiple scattering, including saturation scales, indicated by the big red dots.</figcaption></figure>

even stronger statement [see Eq. (64)]. These proofs are based on a probability interpretation (like $\{1 - W_{11}(1, 1)\}$ being the probability of an interaction at given impact parameter $b$). But this makes no sense if $W_{11}(1, 1)$ has negative values.

One postulates that the qualitative picture discussed above (and therefore a “realistic” description) can be achieved by the following:

(H6a) Solve the sign problem by modifying $g$ in Eq. (13) such that one gets rigorously non-negative results for whatever choice of parameters, without changing the large-$b$ behavior (this will be called “regularization”).

(H7a) Introduce saturation scales when connecting $G$ (used in the multiple scattering formalism) with $G_{QCD}$ (the corresponding QCD expression) such that for hard processes, binary scaling and the AGK theorem are recovered [their very significant violation is an unavoidable consequence of H6a]. The validity of AGK (proven for the case without energy sharing) is a strong expectation which significantly guides the proposed hypotheses.

Concerning H6a, one finds that the infinite sums in Eq. (13) are finally simply a product of exponentials [see Eq. (B64)] and therefore definitely always non-negative, leaving the large-$b$ behavior unchanged as required in H6a), if one postulates the following:

(H6b) For given coefficients $\tilde{\beta}_\lambda$ (arbitrary, but in practice of order unity) and three parameters $c_\mu$, the $g$ functions have to be modified as (referred to as “$g$-factorization”)

$$g\left(\sum_{\lambda} \tilde{\beta}_{\lambda}\right) \mapsto c_1 \prod_{\lambda} c_2 g(c_3 \tilde{\beta}_{\lambda}). \quad (16)$$

As shown in Figs. 23 and 24, this $g$-factorization looks like a reasonable approximation, and for large $b$ it does not change results (concerning $W_{AB}$, for example). But for small $b$ this modification is crucial; it changes results from “nonphysical” to “physical” ($W_{AB}$ always non-negative). It is a fundamental change, forcing the expression to become “physical” (to mimic the elements not treated explicitly).

This is clearly an *ad hoc* solution. But let me discuss some analog to better understand what is done: consider the power series $f(x) = \sum_{i=0}^{\infty} a_i \frac{1}{i!} (-x)^i$. Choosing $a_i = 1$, one gets the well-behaved exponential function shown in Fig. 7 as red squares. If one changes the values $a_i = 1$ randomly by 1% (let me call the corresponding function $f^{(1\%)}$), one gets the result shown as blue circles in Fig. 7: the values strongly fluctuate, and beyond $x = 5$ one gets also negative ones. So even small deviations of the coefficients from the “ideal ones” (namely, $\frac{1}{i!}$) lead to conceptual changes, namely, negative values where the ideal case gives rigorously positive values. Or, the other way around, if one takes the strongly fluctuating function $f^{(1\%)}$, and makes a small (1%) change of the coefficients (back to $a_i = 1$), then one gets the perfectly smooth and positive exponential function. One might see Eq. (13) similarly: one has an “almost exponential” function (with negative arguments as
<figure><img src="image_9.png" /><figcaption>FIG. 7. Exponential function (red squares) and “almost exponential” function (blue circles); see text.</figcaption></figure>

well), and—not surprisingly—one gets negative results. And a small correction (H6a) makes it well behaved. This is why I use the term “regularization”.

But anyway, this regularization is not unique, and this should be further investigated. Already in Eq. (16) one has three parameters. The present choice is one “which works”, but one might consider an optimized set of parameters which minimizes the difference between the factorized expression and the original $g$. One might also simply linearize $g$ as $g(z) \mapsto g(0) + g'(0) \times z$ and then use $g(0) \exp(g'(0)/g(0) \times z)$ which “optimizes” $g$ at small $z$. And then one should see to what extent the final results are affected by these different choices.

The regularization saves the probabilistic interpretation or, to be more precise, expressions, which must represent probabilities (as, for example, $\{1 - W_{11}(1, 1)\}$) behave properly (values between 0 and 1, correct normalization). One has probability laws allowing to generate all possible event configurations and compute all kinds of distributions, all based on the assumption “$G$ is equal to $G_{QCD}$”, the latter being the QCD expression, discussed in detail in Ref. [7] (this is a preliminary assumption, to be changed later). One finds quickly that binary scaling and more generally the AGK theorem are badly violated, and this can be traced back to the so-called deformation function $R_{deform}(x^+, x^-)$ being different from unity in the case of event classes of high event activity. Deformation functions are given as the inclusive distribution of Pomeron light-cone momentum fractions considering particular event classes (EC), $\{\frac{1}{\sigma_{\text{incl}}^{\text{EC}}} \frac{d\sigma_{\text{incl}}^{\text{EC}}}{dx^+ dx^-}\}$, divided by the corresponding expression for an event class with only isolated Pomerons. Having understood the origin of the problem, one parametrizes it; i.e., one parametrizes $R_{deform}(x^+, x^-)$ as a function of a variable called the “connection number”, $N_{\text{conn}} = \frac{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$, where $i$ and $j$ are the projectile and target nucleons a given Pomeron is connected to (see Sec. VII).

One understands that because of energy-momentum sharing, $R_{deform}(x^+, x^-)$ is unavoidably different from unity (actually $<1$ at large $x^+x^-$), and therefore the AGK theorem and binary scaling are unavoidably violated. There is only one way out: giving up $G = G_{QCD}$, which was thought to be the natural choice. In addition, one needs to take care of

H7a, namely, incorporation of saturation scales. Putting all the pieces together, in particular the quantitative understanding of why AGK and binary scaling fail, and how this is related to the deformation function, it is almost mandatory to postulate the following:

(H7b) The $G$ of the multiple scattering formalism and $G_{QCD}$ are related as

$$G(x^+, x^-) = \frac{n}{R_{\text{deform}}(x^+, x^-)} G_{\text{QCD}}(Q_{\text{sat}}^2, x^+, x^-), \quad (17)$$

such that $G$ itself does not depend on the environment (the latter represented via the connection number).

This means that the deformation does not affect $G$ but it changes (very much) the saturation scale $Q_{\text{sat}}^2$, and since the saturation scale affects eventually low-$p_t$ particle production, the high-$p_t$ behavior will not be affected, and one recovers binary scaling at high $p_t$, as it should be. Going through the proofs in Sec. VIII, the choice of H7b seems mandatory,

(H7c) based on the expectation that AGK must hold [leading to factorization (pp) and binary scaling (A+B) at large $p_t$].

But still, H7b and H7c are assumptions.

Let me close this section with an important remark. One should clearly distinguish between “multiple scattering formalism” and “Monte Carlo implementation”. A very important aspect of this paper is the attempt to provide a formalism which is 100% compatible with its Monte Carlo implementation. This sounds obvious, but this is in general not the case. If energy conservation is enforced in the code (via if statements), then this should have a correspondence in the theoretical formulas. It is easy to add an “if” in the code, but very difficult to have the corresponding feature in the theoretical formalism. Concerning all aspects of the current work, the model is completely defined via formulas (finally cross sections expressed in terms of multidimensional integrals), and the Monte Carlo implementation is only a numerical procedure to solve mathematical equations.

To summarize this section: The present work tries to put together theoretical knowledge on S-matrix theory, perturbative QCD, and saturation, together with a large amount of phenomenological (or even qualitative) understanding of experimental data (in particular at the RHIC and LHC), to construct a multiple scattering formalism, which goes beyond factorization, and which allows to address particle production per event classes (a major activity at the LHC for the next decade). Multiple scatterings must happen in parallel, there is no question about that, but the big assumption is the hypothesis that saturation can be incorporated via saturation scales, without any need to split Pomerons, so one can keep the picture of multiple parallel Pomerons. But the fact that saturation is not treated explicitly clearly requires some corrections at some stage. Doing cross-section calculations, one eventually finds expressions, usually interpreted as probabilities, which are negative at small impact parameters (one kind of exceeds
<figure><img src="image_10.png" /><figcaption>FIG. 8. Sketch of a nucleus-nucleus collision in space-time. Here, *t* is the time and *z* the longitudinal coordinate, and the blue lines represent the nucleons.</figcaption></figure>

the black disk limit). Being convinced that these “probability-like” expressions *must* be positive, one forces them to be so. One provides a method which sounds reasonable and which works, but the solution is not unique; one can even propose other possible solutions, and this should be explored further. The final step is the redefinition of the relation between *G* and *G*<sub>QCD</sub> by incorporating saturation scales. One proposes Eq. (17), and from very detailed discussions about the origin of the violation of binary scaling, and the discussion of why and how this proposal solves the problem, it looks like this choice is mandatory. At the end one has constructed a formalism that incorporates perturbative QCD (pQCD) and saturation, does parallel scattering in *pp* and *AA* scattering appropriately, can do analysis of particle production per event class, and is perfectly compatible with the so-called QCD Monte Carlo generators (with fewer hard processes implemented though).

### III. PARALLEL SCATTERING

One of the most important aspects of the formalism developed in this paper is the fact that collisions must happen in parallel—at high energies. In the following the energy dependence of this aspect will be studied. Let me consider a central nucleus-nucleus collision in space-time in its center-of-mass system, as sketched in Fig. 8. The points (intersections between projectile and target nucleon trajectories, assuming that they are close enough in the transverse direction) represent possible interactions. But what precisely happens depends on the relations between three crucial timescales:

$$ \tau_{\text{collision}} = \frac{2R}{\gamma v}, \qquad (18) $$

which is the duration of the *AA* collision,

$$ \tau_{\text{interaction}} = \frac{2R}{n\gamma v}, \qquad (19) $$

which is the time between two *NN* interactions, and

$$ \tau_{\text{form}} = \gamma_{\text{hadron}} \tau_{\text{form, CMS}} \qquad (20) $$

which is the particle (hadron) formation time, after an interaction of two nucleons, with *R* being the nuclear radius, *v* the (average) velocity of the nucleons, *γ* the corresponding Lorentz gamma factor, *n* the (average) number of nucleons in a row, *τ*<sub>form, CMS</sub> the (average) formation time of produced particles in their rest frames, and *γ*<sub>hadron</sub> their (average) gamma factor. The collision energy enters essentially because of the

<figure><img src="image_11.png" /><figcaption>FIG. 9. Sketch of a nucleus-nucleus collision at very low energies, with sequential scatterings.</figcaption></figure>

gamma factors in the denominators of $τ_{\text{collision}}$ and $τ_{\text{interaction}}$, which leads to strong Lorentz contractions at high energies. At very low energies, defined by

$$ \tau_{\text{form}} < \tau_{\text{interaction}}, \qquad (21) $$

after each collision, the particles are formed before the next interaction happens, so one has sequential scatterings (a cascade), as sketched in Fig. 9. The blue lines are the trajectories of projectile and target nucleons, and the red lines are produced hadrons. However, at very high energies, defined by

$$ \tau_{\text{collision}} < \tau_{\text{form}}, \qquad (22) $$

one has completely “parallel scatterings”; all nucleon-nucleon scatterings are realized before particle production starts, as sketched in Fig. 10. The blue lines are the trajectories of projectile and target nucleons; the red lines are produced hadrons. At TeV energies, the longitudinal dimensions of the nuclei are smaller than 0.01 fm/c, so the overlap area is essentially pointlike. In the intermediate energy range, defined by

$$ \tau_{\text{interaction}} < \tau_{\text{form}} < \tau_{\text{collision}}, \qquad (23) $$

one needs a “partially parallel approach”: several (but not all) *NN* scatterings are realized, before the particle production starts, as sketched in Fig. 11. The blue lines are the trajectories of projectile and target nucleons; the red lines are produced hadrons.

Let me define the low and high energy (per nucleon-nucleon pair) thresholds $E_{\text{HE}}$ by the identity $τ_{\text{collision}} = τ_{\text{form}}$ and $E_{\text{LE}}$ by the identity $τ_{\text{form}} = τ_{\text{interaction}}$. The high energy value $E_{\text{HE}}$ is then obtained from $\frac{2R}{\gamma v} = τ_{\text{form}}$. For a rough numerical estimate, I consider central rapidity hadrons ($γ_{\text{hadron}} = 1$) and $τ_{\text{form}} = 1$ fm/c, and a big nucleus with $R = 6.5$ fm,

<figure><img src="image_12.png" /><figcaption>FIG. 10. Sketch of a nucleus-nucleus collision at very high energies, with all nucleon-nucleon scatterings happening in parallel, before particle production starts.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 11. Sketch of a nucleus-nucleus collision at intermediate energies, with only some (but not all) NN scatterings being realized, before the particle production starts.</figcaption></figure>

which gives $\gamma_c^{\nu} = 13$, so one gets

$$E_{HE} \approx 24 \text{ GeV.} \qquad (24)$$

The low energy value $E_{LE}$ is given as $\tau_{form} = \frac{2R}{n\gamma\nu}$, and using $n = 7$, one gets $\gamma_c^{\nu} = \frac{2R}{nc\tau_{form}} \approx \frac{13}{7}$, which gives

$$E_{LE} \approx 4 \text{ GeV.} \qquad (25)$$

To conclude this part: For energies (in the sense of $E = \sqrt{s_{NN}}$) beyond $E_{HE} \approx 24 \text{ GeV}$ one has to employ a “parallel scattering approach”, for energies below $E_{LE} \approx 4 \text{ GeV}$ a pure hadron cascade is appropriate, and in between one needs a “partially parallel scattering approach”. The estimates for $E_{LE}$ and $E_{HE}$ are conservative in the sense that considering $\gamma_{hadron} > 1$ will lead to even smaller values.

In the case of nucleon-nucleon scattering, multiple parton scattering occurs with the corresponding weight increasing with energy. Here again, the $\gamma$ factor plays a crucial role, forcing multiple scatterings to happen in parallel. Looking at a single scattering, I sketch in Fig. 12 a typical example of a hard scattering in the middle of the diagram, preceded by successive parton emission from both sides. In particular, the first emissions lead to very fast partons, with large $\gamma$ factors, and therefore long lifetimes $\gamma\tau$, with $\tau$ of the order of 1 fm/c. This means that at high energies with $\gamma \gg 1$, the duration of the scattering is very large, which makes it impossible to have several such collisions one after the other. They must occur in parallel; one needs a structure as shown in Fig. 13. In principle, this is relevant for energies beyond several tens of GeV, but really important at TeV energies, with a multiple scattering weight being high.

From the above discussion, only based on elementary arguments based on timescales, I conclude that a parallel scattering scenario is mandatory, for pp and AA scattering, beyond center-of-mass energies of several tens of GeV (see Fig. 14).

<figure><img src="image_14.png" /><figcaption>FIG. 12. Sketch of a nucleon-nucleon scattering, with successive parton emissions before the hard process in the middle.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 13. Sketch of a parallel scattering in nucleon-nucleon scattering.</figcaption></figure>

So not only LHC energies are concerned, but even the energies of the RHIC Beam Energy Scan (BES) program.

The appropriate framework to treat strictly parallel scatterings is S-matrix theory, and before coming to the multiple scattering approach, I will first recall some basic facts about the S-matrix.

#### IV. S-MATRIX THEORY: BASIC DEFINITIONS, ASSUMPTIONS, FACTS

Let me recall some basic facts about S-matrix theory, needed to understand the paper. The S-matrix is the representation of the scattering operator $\hat{S}$, i.e., $S_{ab} = \langle a|\hat{S}|b \rangle$ for basis states $|a\rangle$ and $|b\rangle$. The T-matrix is defined as

$$S_{fi} = \delta_{fi} + i(2\pi)^4 \delta(p_f - p_i)T_{fi}, \qquad (26)$$

where $i$ and $f$ refer to initial and final state, respectively, and $p_i$ and $p_f$ are the corresponding momentum four-vectors. The operator $\hat{S}$ must be unitary, $\hat{S}^\dagger\hat{S} = 1$, which means that the scattering does not change the normalization of a state. One adds three “very plausible” hypotheses:

(i) $T_{ii}$ is Lorentz invariant ($\rightarrow$ use Mandelstam variables $s, t$).

(ii) $T_{ii}(s, t)$ is an analytic function of $s$, with $s$ considered to be a complex variable (Hermitian analyticity).

(iii) $T_{ii}(s, t)$ is real on some part of the real axis.

Using the Schwarz reflection principle (a theorem), $T_{ii}(s, t)$ first defined for $\text{Im}s \ge 0$ can be continued in a unique fashion via $T_{ii}(s^*, t) = T_{ii}(s, t)^*$. The cross section is given as a sum

<figure><img src="image_16.png" /><figcaption>FIG. 14. Validity of the parallel scattering scenario. The “hybrid” area refers to a partially parallel scenario.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIG. 15. (a) Elastic scattering diagram and (b) the corresponding cut diagram.</figcaption></figure>

over final states,

$$ \sigma_{\text{tot}} = \frac{1}{w} \sum_f |\mathbf{T}_{fi}|^2 (2\pi)^4 \delta^4(p_f - p_i), \quad (27) $$

with $w = 2s$ for large $s$. Using unitarity and the Schwarz reflection principle, and using here and in the following $\mathbf{T} = \mathbf{T}_{ii}$, one gets

$$ 2s \sigma_{\text{tot}} = \sum_f (2\pi)^4 \delta(p_f - p_i) |\mathbf{T}_{fi}|^2 = \frac{1}{i} \text{disc} \, \mathbf{T}, \quad (28) $$

with $\text{disc} \, \mathbf{T} = \mathbf{T}(s + ie, t) - \mathbf{T}(s - ie, t)$. Interpretation: $\frac{1}{i}\text{disc} \, \mathbf{T}$ can be seen as a so-called “cut diagram” (I will also use cut $\mathbf{T} = \frac{1}{i}\text{disc} \, \mathbf{T}$) with modified Feynman rules [complex conjugate expressions on the right side of the cut, and “intermediate particles” on mass shell via $2\pi\theta(p^0)\delta(p^2 - m^2)$]; see Fig. 15, where I show an elastic scattering diagram and the corresponding cut diagram, the latter being equal to an inelastic amplitude squared.

For applications such as relativistic proton-proton and heavy-ion scattering, one usually assumes purely transverse momentum transfer, i.e., $t = -p_{\perp}^2$, which allows to define a Fourier transform with respect to $p_{\perp}^2$. In addition, it is useful to divide by $2s$, so I define the so-called “impact parameter representation” as

$$ T(s, b) = \frac{1}{8\pi^2 s} \int d^2 p_{\perp} e^{-i\vec{p}_{\perp} \cdot \vec{b}} \mathbf{T}(s, t), \quad (29) $$

with $b = |\vec{b}|$. I also define

$$ G(s, b) = \frac{1}{i} \text{disc} T(s, b) = 2 \text{Im} T(s, b) = \text{cut} T(s, b). \quad (30) $$

Having $t = 0$ as the argument of $\mathbf{T}$ in Eq. (28), and using $\mathbf{T}(s, 0) = \int d^2b \, 2s T(s, b)$ obtained from Eq. (29), one finds

$$ \sigma_{\text{tot}} = \int d^2b \text{cut} T(s, b) = \int d^2b G(s, b), \quad (31) $$

which allows a simple geometrical interpretation:

(i) One may interpret $G$ as the probability of an interaction at impact parameter $b$, provided $G \le 1$.

Throughout the paper, I will systematically use the symbols $\mathbf{T}$, $\mathbf{T}$, and $G$, as defined above. For clarity, I recall the definitions in Table I.

Let me finally recall that cutting a complex diagram [corresponding to cut $T(s, b)$] amounts to summing over all possible

<table><caption>TABLE I. The symbols $\mathbf{T}$, $\mathbf{T}$, and $G$.</caption><tr><td><b>T</b></td><td>Diagonal element of the elastic scattering T-matrix as defined in standard quantum mechanics textbooks, where the asymptotic state is a system of two protons or two nuclei</td></tr><tr><td><i>T</i></td><td>Fourier transform with respect to the transverse momentum exchange of the elastic scattering T-matrix <b>T</b>, divided by 2s (formulas are simpler using this representation)</td></tr><tr><td><i>G</i></td><td>Defined as <i>G</i> = cut <i>T</i> = 2Im<i>T</i> = <sup>1</sup>&frasl;<sub><i>i</i></sub>disc<i>T</i> (where “disc” refers to the variable <i>s</i>), referring to the inelastic process associated with the cut of the elastic diagram corresponding to <i>T</i></td></tr></table>

cuts (cutting rules [25]). I will illustrate what this means in the next section.

## V. AGK CANCELLATIONS, FACTORIZATION, AND BINARY SCALING IN A SIMPLE SCENARIO WITHOUT ENERGY SHARING

To better understand the importance of AGK cancellations, and its relation with factorization and binary scaling, I first discuss a simple scenario, without energy sharing.

### A. A simple scenario

Let me consider the original Gribov-Regge approach for multiple scattering, following Refs. [1–4]. For pp scattering, the T-matrix is given in terms of products of sub-T-matrices,

$$ iT = \sum_{n=0}^{\infty} \frac{1}{n!} \{iT_{\text{Pom}}\}^n, \quad (32) $$

with $T_{\text{Pom}}$ referring to an elementary scattering between “proton constituents” by the exchange of a “Pomeron” (whatever this may be). Originally, lacking an underlying microscopic theory, one simply used asymptotic expressions (limit $s/t \to \infty$) of the form $\mathbf{T}_{\text{Pom}}(s, t) = A s^{B+Ct}$ with parameters A, B, and C (see Appendix A), which allows to easily compute $T_{\text{Pom}}(s, b)$ according to Eq. (29). The limit of infinite energy implies that energy sharing plays no role, and the arguments s, t of $T_{\text{Pom}}$ refer to pp scattering, and not to parton-parton scattering (I will come back to this later). The inelastic cross section is given as [similar to Eq. (31)]

$$ \sigma_{\text{in}} = \int d^2b \text{cut} T(s, b) \quad (33) $$

$$ = \int d^2b \sum_{n=0}^{\infty} \frac{1}{n!} \sum_{\text{cuts}} \{iT_{\text{Pom}}\}^n. \quad (34) $$

Summing all cuts simply means considering a graph with $n$ Pomeron exchanges, seen as a three-dimensional (3D) plot, and considering all possibilities to cut them, with at least one Pomeron being cut. The cuts are indicated by the dashed lines in Fig. 16. A cut should be imagined as an infinite plane. One usually assumes the sub-T-matrix to be purely imaginary, i.e., $T_{\text{Pom}} = i \frac{G}{2}$, with some real function $G = G(b)$, and a factor
<figure><img src="image_18.png" /><figcaption>FIG. 16. Five-Pomeron graph, with two Pomerons being cut.</figcaption></figure>

1/2 for convenience. Then one gets for the cut Pomeron cut $T_{Pom} = 2 \text{Im } T_{Pom} = G$. Concerning the uncut Pomerons, one sums up the contributions where the Pomeron is to the left or to the right of the cut, which gives $\{iT_{Pom}\} + \{iT_{Pom}\}^* = -G$. So cut and uncut Pomerons have opposite signs, and one gets

$$ \sigma_{\text{in}} = \int d^2 b \sum_{n=1}^{\infty} \frac{1}{n!} \sum_{m=0}^{n} \binom{n}{m} G^m(-G)^{n-m}, \quad (35) $$

where $m$ refers to the number of cut Pomerons. One may also write

$$ \sigma_{\text{in}} = \int d^2 b \sum_{m=1}^{\infty} \sum_{l=0}^{\infty} \frac{1}{m! l!} G^m(-G)^l. \quad (36) $$

### B. AGK cancellations in pp scattering

Let me consider inclusive cross sections, like jet cross sections, where $m$ cut-Pomeron events contribute $m$ times more than single Pomeron events, so one gets

$$ \sigma_{\text{incl}} = \int d^2 b \sum_{m=1}^{\infty} \sum_{l=0}^{\infty} m \frac{1}{m! l!} G^m(-G)^l, \quad (37) $$

which is equal to

$$ \sigma_{\text{incl}} = \int d^2 b G \underbrace{\sum_{m=0}^{\infty} \sum_{l=0}^{\infty} \frac{1}{m! l!} G^m(-G)^l}_{1}, \quad (38) $$

and so one finds

$$ \sigma_{\text{incl}} = \int d^2 b G = \int d^2 b \text{ cut } T_{\text{Pom}}, \quad (39) $$

which is an amazing result:

(i) Only a single cut Pomeron (cut $T_{Pom}$) contributes to the inclusive cross section.

(ii) All higher-order terms cancel each other (AGK cancellations).

(iii) This is the AGK theorem [4].

So although there is an infinite number of multiple Pomeron exchange diagrams, with an infinite number of possibilities to cut, they all cancel each other, with the exception of one diagram (see Fig. 17).

<figure><img src="image_19.png" /><figcaption>FIG. 17. The only diagram contributing to inclusive cross sections: a single cut Pomeron.</figcaption></figure>

### C. AGK cancellations in A+B scattering

Let me come to the scattering of two nuclei of mass numbers A and B. The inelastic cross section is a straightforward generalization of the pp result of Eq. (37); one simply needs to add the nuclear geometry in terms of an integration over the positions (in the transverse plane) of the nuclei, in addition to the integration over the nuclear impact parameter $b$, as

$$ \int db_{AB} = \int d^2 b \prod_{i=1}^{A} d^2 b_i^A T_A(b_i^A) \prod_{j=1}^{B} d^2 b_j^B T_B(b_j^B), \quad (40) $$

with the nuclear thickness function

$$ T_A(b) = \int dz \rho_A(\sqrt{b^2 + z^2}), \quad (41) $$

where $\rho_A$ is the nuclear density for nucleus A (and correspondingly for B). Considering all AB possible nucleon-nucleon pairs $k$, and summing over all possible numbers of cut ($m_k$) and uncut ($l_k$) Pomerons, one gets for the cross section $\sigma_{\text{in}}^{AB} = \int db_{AB} \text{cut } T^{AB}(s, b)$ the following expression:

$$ \sigma_{\text{in}}^{AB} = \int db_{AB} \sum_{m_1 l_1} \cdots \sum_{m_{AB} l_{AB}} \prod_{k=1}^{AB} \frac{1}{m_k! l_k!} (G_k)^{m_k} (-G_k)^{l_k}, \quad (42) $$

with at least one $m_k$ being nonzero, and with $G_k = G(b_k)$ and $b_k = |b + b_{\pi(k)}^A - b_{\tau(k)}^B|$ referring to the impact parameter of the nucleon-nucleon pair number $k$.

The inclusive cross section (again like jet cross sections), where pairs $k$ with $m_k$ cut Pomerons count with a factor $m_k$, is given as

$$ \sigma_{\text{incl}}^{AB} = \sum_{m_1 l_1} \cdots \sum_{m_{AB} l_{AB}} \left\{ \sum_{k'=1}^{AB} m_{k'} \right\} \times \int db_{AB} \prod_{k=1}^{AB} \frac{1}{m_k! l_k!} (G_k)^{m_k} (-G_k)^{l_k}, \quad (43) $$

where in $\{\cdots\}$ one sums up the cut-Pomeron numbers from the different nucleon-nucleon collisions. Using $m_{k'}/(m_k)! = 1/(m_{k'} - 1)$, and after renaming $m_{k'} - 1$ to $m_{k'}$ in the product
<figure><img src="image_20.png" /><figcaption>FIG. 18. The only diagram contributing to inclusive cross sections in A+B scattering: a single cut Pomeron (multiplied by AB).</figcaption></figure>

$\prod_{k=1}^{AB}$, one gets

$$ \sigma_{\text{incl}}^{AB} = \sum_{k'=1}^{AB} \int db_{AB} G_{k'} \\ \times \left\{ \sum_{m_1 l_1} \cdots \sum_{m_{AB} l_{AB}} \prod_{k=1}^{AB} \frac{1}{m_k! l_k!} G_k^{m_k} (-G_k)^{l_k} \right\}, \quad (44) $$

without any constraint concerning the $m_k$ summations. Therefore, the term $\{\cdots\}$ is unity, which means that the different multiple Pomeron terms cancel each other (AGK cancellations). The thickness functions are normalized as $\int d^2 b_i^A T_A(b_i^A) = 1$, and correspondingly for B, so all the integrations $\int d^2 b_i^A T_A(b_i^A)$ and $\int d^2 b_j^B T_B(b_j^B)$ in $\int db_{AB}$ give unity, except for the indices $i'$ and $j'$ corresponding to $k'$. So one finds

$$ \sigma_{\text{incl}}^{AB} = \sum_{k'=1}^{AB} \int d^2 b_i^A T_A(b_i^A) \int d^2 b_j^B T_B(b_j^B) \\ \times \int d^2 b G(|b + b_i^A - b_j^B|). \quad (45) $$

A variable change such as $\tilde{b}' = \tilde{b} + \tilde{b}_i^A - \tilde{b}_j^B$ allows one to write

$$ \sigma_{\text{incl}}^{AB} = \sum_{k'=1}^{AB} \int d^2 b_i^A T_A(b_i^A) \int d^2 b_j^B T_B(b_j^B) \int d^2 b' G(b'). \quad (46) $$

After using the normalization $\int d^2 b_i^A T_A(b_i^A) = 1$ and $\int d^2 b_j^B T_B(b_j^B) = 1$, one replaces $\sum_{k'=1}^{AB}$ by simply AB, and one gets

$$ \sigma_{\text{incl}}^{AB} = AB \times \int d^2 b G(b) = AB \times \int d^2 b \text{ cut } T_{\text{Pom}}, \quad (47) $$

which is again an amazing result:

(i) The inclusive A+B cross section is AB times the single cut-Pomeron contribution (cut $T_{\text{Pom}}$) (see Fig. 18).

(ii) This is the AGK theorem [4] for (minimum bias) scattering of nuclei with mass numbers A and B.

(iii) It also implies $\sigma_{\text{incl}}^{AB} = AB \times \sigma_{\text{incl}}^{pp}$, also referred to as binary scaling.

So although there is an infinite number of multiple Pomeron exchange diagrams, with an infinite number of possibilities to cut, they all cancel each other, with the exception of one diagram.

<figure><img src="image_21.png" /><figcaption>FIG. 19. (a) A single cut Pomeron based on pQCD and (b) the corresponding inelastic process.</figcaption></figure>

### D. Factorization

So I have shown that for pp scattering, and also for nuclear collisions, the inclusive cross sections are given in terms of one single cut Pomeron ($G = \text{cut } T_{\text{Pom}}$). What does this mean concerning factorization? To answer this question, one needs to specify the internal structure of the Pomeron. At the time of Gribov et al., it was unknown, and one could not do more than using parametrized expressions for the T-matrices, based on their asymptotic behavior. But nowadays one knows more. It is shown in Ref. [7] how to compute T-matrices and then $G = \text{cut } T$, corresponding to (a single) parton-parton scattering, based on perturbative QCD. The result is named $G_{\text{QCD}}$. It is a sum of several terms, the most important one shown in Fig. 19, representing an integration expressed in terms of several modules ($F_{\text{sea}}$, $E_{\text{psoft}}$, $E_{\text{QCD}}$, Born), with $E_{\text{QCD}}$ representing a spacelike parton evolution, and “Born” the QCD Born process. Assuming that the $G$ of the multiple scattering formalism is equal to $G_{\text{QCD}}$, i.e.,

$$ G = G_{\text{QCD}}, \quad (48) $$

one gets for the inclusive pp cross section (using the AGK theorem, and therefore considering only one single Pomeron) the following result (see Ref. [7]):

$$ \sigma_{\text{incl}}^{pp} = \int \frac{d^3 p_3 d^3 p_4}{E_3 E_4} \\ \times \left\{ \sum_{klmn} \iint d\xi_1 d\xi_2 f_{\text{PDF}}^k(\xi_1, \mu_F^2) f_{\text{PDF}}^l(\xi_2, \mu_F^2) \right. \\ \left. \times \frac{1}{32s\pi^2} \sum |\mathcal{M}^{kl \to mn}|^2 \delta^4(p_{\text{in}} - p_{\text{out}}) \frac{1}{1 + \delta_{mn}} \right\}, \quad (49) $$

with $f_{\text{PDF}}^k$ being a convolution of modules like $F_{\text{sea}}$, $E_{\text{psoft}}$, and $E_{\text{QCD}}$, with $\mathcal{M}^{kl \to mn}$ being the QCD matrix element for elementary parton-parton $2 \to 2$ scatterings (Born process), and with $\vec{p}_3$, $E_3$ and $\vec{p}_4$, $E_4$ being the momenta and energies of the outgoing partons. This amounts to factorization.
E. Conclusion

It is interesting that one can derive in a quantum mechan-
ical multiple scattering approach important features such as
binary scaling and factorization, as a consequence of AGK
cancellations. It is very useful to understand these phenom-
ena in a simple approach where calculations can be easily
done. But the main problem is the fact that energy sharing
is not present. The cut-Pomeron expressions G should depend
explicitly on the Pomeron energy s<sub>Pom</sub>, i.e., G = G(s<sub>Pom</sub>, b),
and when one has ten Pomerons in a pp scattering, then the
available energy has to be shared among these ten Pomerons.
It is also known that binary scaling fails at low p<sub>t</sub>, whereas
here binary scaling is always true. So the scenario discussed
in this section is a first step, but not the final solution.

But even in more realistic scenarios, to be discussed later,
one feature remains always correct: the fundamental property
is the “validity of AGK cancellations”, for *pp* scatterings, and
for AA scatterings, whereas binary scaling and factorization
are merely consequences.

To summarize this section:

(i) In a simplified parallel scattering S-matrix approach, without energy sharing (and consequently simple formulas), one can prove the “AGK theorem”, in *pp* and AA scattering (first shown in Ref. [4]).

(ii) “AGK theorem” means for *pp* scattering that inclusive cross sections $\sigma_{incl}^{pp}$ are given in terms of a single cut Pomeron, although the real events represent multiple scatterings.

(iii) “AGK theorem” means for $A+B$ scattering that inclusive cross sections $\sigma_{incl}^{AB}$ are given as $AB$ times the single cut-Pomeron contribution.

(iv) As a consequence, one gets $\sigma_{incl}^{AB} = AB \times \sigma_{incl}^{pp}$, which is called binary scaling. So it is a direct consequence of the AGK theorem.

(v) One can make statements about factorization, but
only when one specifies the internal structure of
the Pomeron, more precisely the microscopic pic-
ture underlying T<sub>Pom</sub> and G = cut T<sub>Pom</sub>. Assuming
G = G<sub>QCD</sub>, the latter being based on a QCD calcula-
tion of parton-parton scattering, including a DGLAP
evolution as defined in Ref. [7], one can deduce
"factorization".

(vi) The AGK theorem is the fundamental feature; others (like binary scaling and factorization) are just consequences.

(vii) The simple scenario of this section is useful for un-
derstanding the relation between the AGK theorem,
binary scaling, and factorization. However, a major
problem is the fact that energy sharing is not present,
but it should be.

In the following, I will discuss a scenario including energy
sharing.

<figure><img src="image_22.png" /><figcaption>FIG. 20. (a) Multiple scattering with energy sharing, for two col-
liding nuclei A and B, with Pomerons (cyan vertical lines), remnants
(pink horizontal lines), vertices (magenta dots), and cuts (black,
vertical). (b) Detailed view of a vertex connected to projectile i,
associated to pair k, with remnants light-cone momentum fraction
$x_{\text{remn},i}^{+}$ and Pomeron leg momenta $x_{k\mu}^{+}$.</figcaption></figure>

VI. MULTIPLE SCATTERING S-MATRIX APPROACH WITH ENERGY SHARING AND “SIGN PROBLEM”

In this section, I take over the S-matrix approach discussed
in Sec. V (GR approach), but I add energy-momentum sharing
(GR$^+$). I will discuss right away collisions of nuclei with mass
numbers A and B, where pp is just a special case (A = B = 1).

A. Scenario with energy-momentum sharing

I generalize Eq. (42) by taking into account energy-
momentum sharing, as sketched in Fig. 20(a). In general, one
considers all possible interactions between the *A* projectile
and the *B* target nucleons, and one therefore introduces a
“pair index” *k*, going from 1 to *AB*. I define *i* = π(*k*) and
*j* = τ(*k*) to be the associated nucleon indices on the projectile
and target sides. In Fig. 20(a), I consider *A* = *B* = 2 (for sim-
plicity), and I consider multiple (three) Pomeron interactions
between *i* = 1 and *j* = 1 and between *i* = 2 and *j* = 2, which
means zero Pomerons exchanged between *i* = 1 and *j* = 2
and between *i* = 2 and *j* = 1. In general, one sums over any
number of Pomerons for all possible pairs. In Fig. 20(b), I
show a detailed view of one of the projectile vertices con-
connected to projectile *i*, with the associated pair numbering *k*.
For the different components, one uses light-cone momentum
fractions *x*± with respect to the initial light-cone momenta
of the nucleons, with 0 ≤ *x*<sup>+</sup>, *x*<sup>-</sup> ≤ 1. I define *x*<sub>*kμ*</sub><sup>±</sup> to be the
light-cone momentum fractions of the external legs of the
μth Pomeron of pair *k*, all of them connected to projectile
nucleon *i* = π(*k*) and target nucleon *j* = τ(*k*). For each cut
Pomeron one has a term *G* (equal to cut *T*<sub>Pom</sub> = 2Im *T*<sub>Pom</sub>;
see Sec. IV for basic definitions), and for each uncut Pomeron
one has a term −*G*, both with arguments *x*<sub>*kμ*</sub><sup>±</sup>. The light-cone
momenta of the projectile and target remnants are named
*x*<sub>remn,*i*</sub><sup>+</sup> and *x*<sub>remn,*j*</sub><sup>−</sup>. In Fig. 20(b), I only show the projectile part.
Energy-momentum conservation amounts to

$$x_{\text{remn},i}^{+} = 1 - \sum_{k=1}^{AB} \sum_{\mu=1}^{n_k} x_{k\mu}^{+}, \quad x_{\text{remn},j}^{-} = 1 - \sum_{k=1}^{AB} \sum_{\mu=1}^{n_k} x_{k\mu}^{-} \qquad (50)$$

which means that the initial values ($x^{\pm} = 1$) are shared among Pomerons and remnants. The nuclear geometry will be taken care of by using $\int db_{AB}$, defined in Eqs. (40) and (41), which integrates over the impact parameter and the transverse co-ordinates $b_i^A$ and $b_j^B$ of projectile and target nucleons. With all these definitions, one may write the expression for the inelastic cross section being a sum over all possible numbers of Pomerons per pair $k$ and a sum over all possibilities to cut them, as

$$\sigma_{\text{in}}^{AB} = \sum_{n_1=0}^{\infty} \cdots \sum_{n_{AB}=0}^{\infty} \sum_{m_1 \le n_1} \cdots \sum_{m_{AB} \le n_{AB}} \int db_{AB} \\ \times \int dX_{AB} \prod_{k=1}^{AB} \frac{1}{n_k!} \left( \frac{n_k}{m_k} \right) \prod_{\mu=1}^{m_k} G_{k\mu} \prod_{\mu=m_k+1}^{n_k} -G_{k\mu} \\ \times \prod_{i=1}^{A} V(x_{\text{remn},i}^{+}) \prod_{j=1}^{B} V(x_{\text{remn},j}^{-}) \qquad (51)$$

[generalizing Eq. (42)], with at least one $n_k$ being nonzero, with $\prod_{\mu=1}^{m_k} G_{k\mu} = 1$ for $m_k = 0$, and with

$$G_{k\mu} = G(x_{k\mu}^{+}, x_{k\mu}^{-}, s, b_{k}), \qquad (52)$$

with $b_k = |b + b_{\pi(k)}^A - b_{\tau(k)}^B|$ referring to the impact parameter of the nucleon-nucleon pair number $k$. The indices $n_k$ refer to the numbers of Pomerons (cut and uncut) for pair $k$, whereas the indices $m_k$ refer to the number of cut Pomerons. The symbol $\int dX_{AB}$ represents the integration over light-cone momentum fractions:

$$\int dX_{AB} = \int \prod_{k=1}^{AB} \prod_{\nu=1}^{n_k} dx_{k\nu}^{+} dx_{k\nu}^{-} \qquad (53)$$

$$\frac{d^2 \sigma_{\text{incl}}^{AB}}{dx^+ dx^-} = \sum_{k'=1}^{AB} \sum_{n_1=0}^{\infty} \cdots \sum_{n_{AB}=0}^{\infty} \sum_{m_1 \le n_1} \cdots \sum_{m_{AB} \le n_{AB}} \sum_{\mu'=1}^{m_{k'}} \int db_{AB} \int dX_{AB} \\ \times \left\{ \prod_{k=1}^{AB} \left[ \frac{1}{m_k!(n_k - m_k)!} \prod_{\mu=1}^{m_k} G_{k\mu} \prod_{\mu=m_k+1}^{n_k} - G_{k\mu} \right] \prod_{i=1}^{A} V(x_{\text{remn},i}^{+}) \prod_{j=1}^{B} V(x_{\text{remn},j}^{-}) \delta(x^{+} - x_{k'\mu'}^{+}) \delta(x^{-} - x_{k'\mu'}^{-}) \right\}, \qquad (55)$$

with at least one $n_k$ being nonzero. The $\delta$ functions make some of the integrations of $\int dX_{AB}$ trivial. Integrating over $x_{k'\mu'}^+$ and $x_{k'\mu'}^-$ allows one to replace $G_{k',\mu'}$ by $G(x^+, x^-, s, b_{k'})$, which may be written in front of $\int dX_{AB}$. Then one may rename the integration variables such that the variables $x_{k'\mu'}^+$ disappear (to have cut-Pomeron variables $x_{k'\mu'}^+$ with $\mu \le m_k - 1$ and uncut-Pomeron variables $x_{k'\mu'}^-$ with $m_k' \le \mu \le n_k - 1$). In $x_{\text{remn},i}^+$ for $i = \pi(k')$ and in $x_{\text{remn},j}^-$ for $j = \tau(k')$ one replaces 1 by $1 - x^+$ and $1 - x^-$, respectively. The same procedure applies for all

Equation (50) ensures energy-momentum conservation. If one would remove the vertex part (replace $V$ by unity), the $\int dX_{AB}$ could be done, and one recovers Eq. (42), the case without energy sharing (I just organize the summations somewhat differently). With the vertex part, and full energy-momentum conservation, the up to 10 000 000 dimensional integrations are not separable, but this “technical problem” can be handled: see Ref. [5] for a detailed discussion about how to compute expressions like Eq. (51).

## B. AGK theorem

As discussed in Sec. V, a crucial feature of any multiple scattering approach is the validity of the AGK theorem, which is needed to get for inclusive cross sections binary scaling in AA collisions and factorization in pp. Therefore I am going to investigate this in the following.

Again I am considering collisions of nuclei with mass numbers $A$ and $B$, where pp is just a special case ($A = B = 1$). The inclusive cross section is a modification of the inelastic cross section in Eq. (51). In the case of “no energy sharing” as discussed in Sec. V, I simply added a factor $\{\sum_{k'=1}^{AB} m_{k'}\}$, which amounts to counting the number of cut Pomerons. Now, I introduce energy-momentum sharing, and each Pomeron is characterized by the light-cone momentum fractions $x^+$ and $x^-$, such that the squared (transverse) mass of the Pomeron is given by $x^+x^-s$, with $s$ being the Mandelstam $s$ of the nucleon-nucleon collisions. So one needs to count the Pomerons for given $x^+$ and $x^-$ intervals, which can be done by adding a factor

$$\sum_{k'=1}^{AB} \sum_{\mu'=1}^{m_k'} \delta(x^{+} - x_{k'\mu'}^{+}) \delta(x^{-} - x_{k'\mu'}^{-}) dx^{+} dx^{-} \qquad (54)$$

into Eq. (51), and after dividing by $dx^+dx^-$, one gets

values of $\mu'$, giving identical expressions, so the sum $\sum_{\mu'=1}^{m_{k'}}$ gives simply a factor $m_{k'}$, which one uses to replace $m_{k'}!$ by $(m_{k'} - 1)!$. The expression $[\cdots]$ in the second line of Eq. (55) for $k = k'$ may then be written as

$$\frac{1}{(m_k - 1)!(n_k - 1 - (m_k - 1))!} \prod_{\mu=1}^{m_k-1} G_{k\mu} \prod_{\mu=m_k-1+1}^{n_k-1} -G_{k\mu}. \qquad (56)$$
I then rename $(m_k - 1)$ to $m_k$ and $(n_k - 1)$ to $n_k$, which allows one to drop the condition “at least one $n_k$ nonzero”. In addition, one may now write

$$ \sum_{m_1 \le n_1} \cdots \sum_{m_{AB} \le n_{AB}} \cdots \prod_{k=1}^{AB} \prod_{k=1}^{AB} \sum_{m_k \le n_k}, \quad (57) $$

and so one gets

$$ \begin{aligned} \frac{d^2 \sigma_{\text{incl}}^{AB}}{dx^+ dx^-} = & \sum_{k'=1}^{AB} \sum_{n_1=0}^{\infty} \cdots \sum_{n_{AB}=0}^{\infty} \int db_{AB} G(x^+, x^-, s, b_{k'}) \\ & \times \int dX_{AB} \prod_{k=1}^{AB} \frac{1}{n_k!} \prod_{\mu=1}^{n_k} G_{k\mu} \sum_{m_k \le n_k} \binom{n_k}{m_k} (-1)^{n_k - m_k} \\ & \times \prod_{i=1}^{A} V(x_{\text{remn},i}^+) \prod_{j=1}^{B} V(x_{\text{remn},j}^-), \end{aligned} \quad (58) $$

without any constraints on the $n_k$. The relation

$$ \sum_{m_k \le n_k} \binom{n_k}{m_k} (-1)^{n_k - m_k} = (1-1)^{n_k} = \delta_{0n_k}, \quad (59) $$

for all values of $k$, means that only one term in the sum $\sum_{n_1} \cdots \sum_{n_{AB}}$ contributes, namely, the one with $n_1 = n_2 = \cdots = n_{AB} = 0$, and one gets

$$ \frac{d^2 \sigma_{\text{incl}}^{AB}}{dx^+ dx^-} = \sum_{k'=1}^{AB} \int db_{AB} G(x^+, x^-, s, b_{k'}) \times V(1-x^+)V(1-x^-). \quad (60) $$

Note that only one pair of vertices $V$ contributes, namely, that corresponding to the projectile $i = \pi(k')$ and to the target $j = \tau(k')$; for all others one has $V = V(1)$ which is per definition unity. All the integrations $\int d^2 b_i^A T_A(b_i^A)$ and $\int d^2 b_j^B T_B(b_j^B)$ in $\int db_{AB}$ give unity (normalization of the thickness functions), except for the indices $i'$ and $j'$ corresponding to $k'$, so one gets

$$ \frac{d^2 \sigma_{\text{incl}}^{AB}}{dx^+ dx^-} = \sum_{k'=1}^{AB} \int d^2 b_i^A T_A(b_i^A) \int d^2 b_j^B T_B(b_j^B) \times \int d^2 b G(x^+, x^-, s, b_{k'}) V(1-x^+)V(1-x^-), \quad (61) $$

with $\vec{b}_{k'} = \vec{b} + \vec{b}_i^A - \vec{b}_j^B$. A variable change such as $\vec{b}' = \vec{b} + \vec{b}_i^A - \vec{b}_j^B$ allows one to separate the three integrations, to use again the normalization of the thickness functions, and to replace $\sum_{k'=1}^{AB}$ by $AB$, to get

$$ \frac{d^2 \sigma_{\text{incl}}^{AB}}{dx^+ dx^-} = AB \int d^2 b G(x^+, x^-, s, b) V(1-x^+)V(1-x^-). \quad (62) $$

One recalls that $G(x^+, x^-, s, b)$ is a single cut Pomeron, so one may write

$$ \frac{d^2 \sigma_{\text{incl}}^{AB}}{dx^+ dx^-} = AB \times \int d^2 b \text{ cut } T_{\text{Pom}} \times V(1-x^+)V(1-x^-), \quad (63) $$

<figure><img src="image_23.png" /><figcaption>FIG. 21. The only diagram contributing to inclusive cross sections in A+B scattering: a single cut Pomeron (multiplied by AB).</figcaption></figure>

and again one finds this important result:

(i) The inclusive A+B cross section is $AB$ times the single cut-Pomeron contribution (cut $T_{\text{Pom}}$) including a vertex part (i.e., the single Pomeron inclusive cross section).

(ii) So one has

$$ \frac{d\sigma_{\text{incl}}^{AB}}{dx^+ dx^-} = AB \times \frac{d\sigma_{\text{incl}}^{\text{single Pom}}}{dx^+ dx^-} \quad (64) $$

(see Fig. 21), referred to as the AGK theorem (actually being an extension of the original AGK theorem, taking into account energy-momentum sharing).

(iii) As a direct consequence, one gets

$$ \frac{d\sigma_{\text{incl}}^{AB}}{dx^+ dx^-} = AB \times \frac{d\sigma_{\text{incl}}^{pp}}{dx^+ dx^-}, \quad (65) $$

referred to as binary scaling.

So I have proven in the scenario with energy sharing, for (minimum bias) scattering of nuclei with mass numbers A and B: (i) the validity of the AGK theorem, and (ii) binary scaling. So although there is an infinite number of multiple Pomeron exchange diagrams, with an infinite number of possibilities to cut, they all cancel each other, with the exception of one diagram.

I have been discussing the inclusive cross section $d\sigma_{\text{incl}}^{AB}/dx^+dx^-$, which counts the number of Pomeron with given values of $x^+$ and $x^-$. This is kind of a “master distribution”, which allows one to obtain the inclusive cross section for the production of measurable quantities like transverse momentum, provided the nature of the Pomeron is known, as discussed in Sec. V.

But before implementing the microscopic structure of the Pomeron, one has to confront a quite serious problem, related to probabilistic interpretations, as discussed in the following.

### C. The sign problem

Let me consider again the inelastic cross section for the scattering of two nuclei, which may be written as [see Eq. (51)]

$$ \sigma_{\text{in}}^{AB} = \int db_{AB} \text{ cut } T^{AB}(s, \{b_{AB}\}) = \int db_{AB} G^{AB}(s, \{b_{AB}\}), \quad (66) $$
<figure><img src="image_24.png" /><figcaption>FIG. 22. The multidimensional variable {$b, \{b_i^A\}, \{b_j^B\}$}, with b being the impact parameter, and with $b_i^A$ and $b_j^B$ being the transverse coordinates of the projectile and target nucleons, respectively (they are meant to be two-dimensional vectors).</figcaption></figure>

where one employs the usual relation $G = \text{cut } T$. Here, {$b_{AB}$} is the multidimensional variable {$b, \{b_i^A\}, \{b_j^B\}$} (see Fig. 22), and $\int db_{AB}$ the corresponding integration, defined in Eqs. (40) and (41). Equation (66) indicates that $G^{AB}$ may be interpreted as the probability of having an interaction for a given impact parameter b and given transverse coordinates {$b_i^A$} and {$b_j^B$} of the nucleons. Comparing with Eq. (51), one gets

$$ G^{AB}(s, \{b_{AB}\}) = \sum_{n_1=0}^{\infty} \cdots \sum_{n_{AB}=0}^{\infty} \sum_{m_1 \le n_1} \cdots \sum_{m_{AB} \le n_{AB}} \int dX_{AB} \\ \times \prod_{k=1}^{AB} \frac{1}{n_k!} (n_k) \prod_{\mu=1}^{m_k} G_{k\mu} \prod_{\mu=m_k+1}^{n_k} -G_{k\mu} \\ \times \prod_{i=1}^{A} V(x_{\text{remn},i}^{+}) \prod_{j=1}^{B} V(x_{\text{remn},j}^{-}), \quad (67) $$

with at least one $n_k$ being nonzero, and with the momentum fraction integrations $\int dX_{AB}$ being defined in Eq. (53). One may define an expression $\tilde{G}^{AB}(s, \{b_{AB}\})$ as Eq. (67), but for all $n_k = 0$, which represents the case excluded in the summation of Eq. (67), which might be interpreted as the probability of “no interaction”. One then gets

$$ G^{AB} + \tilde{G}^{AB} = \sum_{n_1=0}^{\infty} \cdots \sum_{n_{AB}=0}^{\infty} \sum_{m_1 \le n_1} \cdots \sum_{m_{AB} \le n_{AB}} \int dX_{AB} \\ \times \prod_{k=1}^{AB} \frac{1}{n_k!} (n_k) \prod_{\mu=1}^{m_k} G_{k\mu} \prod_{\mu=m_k+1}^{n_k} -G_{k\mu} \\ \times \prod_{i=1}^{A} V(x_{\text{remn},i}^{+}) \prod_{j=1}^{B} V(x_{\text{remn},j}^{-}), \quad (68) $$

without any constraint for the $n_k$. Using similar manipulations as for the $d\sigma_{\text{incl}}/dx^+dx^-$ calculation above, one gets

$$ G^{AB}(s, \{b_{AB}\}) + \tilde{G}^{AB}(s, \{b_{AB}\}) = \prod_{i=1}^{A} V(1) \prod_{j=1}^{B} V(1), \quad (69) $$

and with (by definition) $V(1) = 1$, one finds that the right-hand side (rhs) of this equation is unity.

So I have proven

$$ G^{AB}(s, \{b_{AB}\}) + \tilde{G}^{AB}(s, \{b_{AB}\}) = 1, \quad (70) $$

with $G^{AB}(s, \{b_{AB}\})$ and $\tilde{G}^{AB}(s, \{b_{AB}\})$ being interpreted as probabilities to interact for a given impact parameter b and given transverse coordinates {$b_i^A$} and {$b_j^B$} of the nucleons, and to not interact, respectively. This sounds promising. So far the interpretations seem to be consistent.

From Eqs. (70) and (68), one gets

$$ 1 = \sum_{n_1=0}^{\infty} \cdots \sum_{n_{AB}=0}^{\infty} \sum_{m_1 \le n_1} \cdots \sum_{m_{AB} \le n_{AB}} \int dX_{AB} \\ \times \prod_{k=1}^{AB} \frac{1}{n_k!} (n_k) \prod_{\mu=1}^{m_k} G_{k\mu} \prod_{\mu=m_k+1}^{n_k} -G_{k\mu} \\ \times \prod_{i=1}^{A} V(x_{\text{remn},i}^{+}) \prod_{j=1}^{B} V(x_{\text{remn},j}^{-}). \quad (71) $$

Let me rename the light-cone momentum variables: let $x_{k\mu}^{\pm}$ be the light-cone momenta of the $\mu$th cut Pomeron for pair k, and $\tilde{x}_{k\lambda}^{\pm}$ those of the $\lambda$th uncut Pomeron for pair k. Instead of $\int dX_{AB}$, one has $\int dX_{AB} \int d\tilde{X}_{AB}$, where now

$$ \int dX_{AB} = \int \prod_{k=1}^{AB} \left( \prod_{\mu=1}^{m_k} dx_{k\mu}^{+} dx_{k\mu}^{-} \right) \quad (72) $$

refers to cut Pomerons and

$$ \int d\tilde{X}_{AB} = \int \prod_{k=1}^{AB} \left( \prod_{\lambda=1}^{l_k} d\tilde{x}_{k\lambda}^{+} d\tilde{x}_{k\lambda}^{-} \right) \quad (73) $$

to uncut Pomerons. I use $G'_{k\lambda} = G(\tilde{x}_{k\lambda}^{+}, \tilde{x}_{k\lambda}^{-}, s, b_k)$. All this allows one to reorganize the sums in Eq. (71), and one may write

$$ 1 = \sum_{m_1=0}^{\infty} \cdots \sum_{m_{AB}=0}^{\infty} \sum_{l_1=0}^{\infty} \cdots \sum_{l_{AB}=0}^{\infty} \int dX_{AB} \int d\tilde{X}_{AB} \\ \times \prod_{k=1}^{AB} \left[ \frac{1}{m_k!} \prod_{\mu=1}^{m_k} G_{k\mu} \right] \prod_{k=1}^{AB} \left[ \frac{1}{l_k!} \prod_{\lambda=1}^{l_k} -G'_{k\lambda} \right] \\ \times \prod_{i=1}^{A} V(x_{\text{remn},i}^{+}) \prod_{j=1}^{B} V(x_{\text{remn},j}^{-}), \quad (74) $$

where $m_k$ refers to the number of cut Pomerons and $l_k$ to the number of uncut ones, for pair k. One recalls that $G_{k\mu}$ is defined as $G(x_{k\mu}^{+}, x_{k\mu}^{-}, s, b_k)$ with $b_k = |b + b_{\pi(k)}^{A} - b_{\tau(k)}^{B}|$. I define “cut Pomeron configurations”

$$ K = \{\{m_k\}, \{x_{k\mu}^{\pm}\}\}, \quad (75) $$

with $\{m_k\}$ being the set of the $m_k$ variables, and $\{x_{k\mu}^{\pm}\}$ the set of the $x_{k\mu}^{\pm}$ variables, for given values of the nuclear impact parameter b and the transverse coordinates $b_i^A$ and $b_j^B$ of all the nucleons in nuclei A and B. The sum $\sum_{m_1=0}^{\infty} \cdots \sum_{m_{AB}=0}^{\infty}$ is then simply the sum over all possible values of $\{m_k\}$, and Eq. (74) becomes

$$ 1 = \sum_{\{m_k\}} \int dX_{AB} P, \quad (76) $$
or explicitly

$$1 = \sum_{\{m_k\}} \int \prod_{k=1}^{AB} \prod_{\mu=1}^{m_k} dx_{k\mu}^{+} dx_{k\mu}^{-} P(\{m_k\}, \{x_{k\mu}^{\pm}\}), \quad (77)$$

with

$$P(\{m_k\}, \{x_{k\mu}^{\pm}\}) = \prod_{k=1}^{AB} \left[ \frac{1}{m_k!} \prod_{\mu=1}^{m_k} G_{k\mu} \right] \times W_{AB}, \quad (78)$$

with

$$W_{AB} = \sum_{l_1=0}^{\infty} \cdots \sum_{l_{AB}=0}^{\infty} \int d\tilde{X}_{AB} \prod_{k=1}^{AB} \left[ \frac{1}{l_k!} \prod_{\lambda=1}^{l_k} -G'_{k\lambda} \right] \\ \times \prod_{i=1}^{A} V \left( x_i^+ - \sum_{k=1}^{AB} \sum_{\lambda=1}^{l_k} \tilde{x}_{k\lambda}^+ \right) \prod_{j=1}^{B} V \left( x_j^- - \sum_{k=1}^{AB} \sum_{\lambda=1}^{l_k} \tilde{x}_{k\lambda}^- \right), \quad (79)$$

with

$$x_i^+ = 1 - \sum_{\substack{k=1 \\ \pi(k)=i}}^{AB} \sum_{\mu=1}^{m_k} x_{k\mu}^+, \quad x_j^- = 1 - \sum_{\substack{k=1 \\ \tau(k)=j}}^{AB} \sum_{\mu=1}^{m_k} x_{k\mu}^-. \quad (80)$$

The term $W_{AB}$ is a function of the variables $\{x_i^+\}$ with $1 \le i \le A$ and $\{x_j^-\}$ with $1 \le j \le B$, i.e., $W_{AB} = W_{AB}(\{x_i^+\}, \{x_j^-\})$. It also depends on $s$ and $b$ and the transverse coordinates $b_i^A$ and $b_j^B$ of all the nucleons, not written explicitly.

One concludes:

(i) the equation $1 = \sum_{\{m_k\}} \int dX_{AB} P$ allows the interpretation of

$$P = \prod_{k=1}^{AB} \left[ \frac{1}{m_k!} \prod_{\mu=1}^{m_k} G_{k\mu} \right] \times W_{AB}(\{x_i^+\}, \{x_j^-\}) \quad (81)$$

to be the probability of the “configuration” K,

(ii) where $K = \{\{m_k\}, \{x_{k\mu}^{\pm}\}\}$ represents $m_k$ cut Pomerons per pair $k$, with light-cone momentum fractions $x_{k\mu}^{\pm}$.

This would be a perfect basis for Monte Carlo applications: one “simply” needs to generate configurations according to that law. But what is still missing: one needs to prove that all the $P(\{m_k\}, \{x_{k\mu}^{\pm}\})$ are non-negative, and, to do so, one has to prove that $W_{AB}$ is non-negative.

As a side remark (just to understand the meaning of the quantity $W_{AB}$), one has

$$W_{AB}(\{x_i^+ = 1\}, \{x_j^- = 1\}) = P(\{m_k = 0\}), \quad (82)$$

corresponding to no cut Pomerons at all, which means “no interaction” in the sense of inelastic scattering. The expression $(1 - W_{AB})$, with all momentum fractions $x_i^+$ and $x_j^-$ being unity, corresponds to inelastic scattering, summed over all possible configurations, and the inelastic cross section is

obtained by integrating $\int db_{AB} \cdots$, as

$$\sigma_{\text{in}} = \int db_{AB} \{1 - W_{AB}(\{x_i^+ = 1\}, \{x_j^- = 1\})\}. \quad (83)$$

For example, the pp inelastic cross section is an integral over the impact parameter of $1 - W_{11}(1, 1)$. This again confirms the probabilistic interpretation, but still, under the condition $W_{AB} \ge 0$.

Let me now discuss how to calculate $W_{AB}$. To do so, one needs to specify $G(x^+, x^-, s, b)$. Using $G = G_{QCD}$ [see Eq. (48) and the discussion before], one can show [5] that one may obtain an almost perfect fit of the numerically computed functions $G_{QCD}$, with a parametrization of the form

$$G_{\text{QCDpar}}(x^+, x^-, s, b) = \sum_{N=1}^{N_{\text{par}}} \alpha_N (x^+ x^-)^{\beta_N}, \quad (84)$$

where $\alpha_N$ and $\beta_N$ depend on $s$ and $b$ given in terms of a few parameters, as shown in Appendix B. This parametric form has been inspired by the asymptotic expressions for T-matrices (see Appendix A). So I will use in the following

$$G = G_{\text{QCDpar}}. \quad (85)$$

Furthermore, the vertices are parametrized as

$$V(x) = x^{\alpha_{\text{remn}}}, \quad (86)$$

again motivated by the asymptotic expressions for T-matrices.

Now one has all the ingredients to compute $W_{AB}$. All the formulas will be given with the current choice of $N_{\text{par}} = 4$; it is straightforward to change this value. From Eq. (79), one gets

$$W_{AB}(\{x_i^+\}, \{x_j^-\}) = \sum_{\{l_k\}} \int \prod_{k=1}^{AB} \left( \prod_{\lambda=1}^{l_k} d\tilde{x}_{k\lambda}^+ d\tilde{x}_{k\lambda}^- \right) \\ \times \left\{ \prod_{k=1}^{AB} \left[ \frac{1}{l_k!} \prod_{\lambda=1}^{l_k} -G_{\text{QCDpar}}(\tilde{x}_{k\lambda}^+, \tilde{x}_{k\lambda}^-, s, b_k) \right] \right. \\ \times \prod_{i=1}^{A} \left( x_i^+ - \sum_{\substack{k=1 \\ \pi(k)=i}}^{AB} \sum_{\lambda=1}^{l_k} \tilde{x}_{k\lambda}^+ \right)^{\alpha_{\text{remn}}} \\ \left. \times \prod_{j=1}^{B} \left( x_j^- - \sum_{\substack{k=1 \\ \tau(k)=j}}^{AB} \sum_{\lambda=1}^{l_k} \tilde{x}_{k\lambda}^- \right)^{\alpha_{\text{remn}}} \right\}, \quad (87)$$

where $\sum_{\{l_k\}}$ means summing all the indices $l_k$, with $1 \le k \le AB$, from zero to infinity, where $l_k$ refers to the number of uncut Pomerons of nucleon-nucleon pair $k$. Remember that computing $W_{AB}$ amounts to summing over and integrating out all uncut Pomerons. As proven in Appendix B [see Eq. (B56)],
one finds

$$W_{AB}(\{x_i^+\}, \{x_j^-\}) = \prod_{i=1}^{A} (x_i^+)^{\alpha_{\text{remn}}} \prod_{j=1}^{B} (x_j^-)^{\alpha_{\text{remn}}} \sum_{\{r_{Nk}\}} \left\{ \prod_{k=1}^{AB} \prod_{N=1}^{4} \frac{(-\alpha_N)^{r_{Nk}}}{r_{Nk}!} \times \prod_{i=1}^{A} \left[ \prod_{k=1}^{AB} \prod_{N=1}^{4} (\Gamma(\tilde{\beta}_N)(x_i^+)^{\tilde{\beta}_N})^{r_{Nk}} g\left(\sum_{\pi(k)=i}^{AB} \sum_{\pi(k)=j}^{4} r_{Nk} \tilde{\beta}_N\right) \right] \times \prod_{j=1}^{B} \left[ \prod_{k=1}^{AB} \prod_{\tau(k)=j}^{4} (\Gamma(\tilde{\beta}_N)(x_j^-)^{\tilde{\beta}_N})^{r_{Nk}} g\left(\sum_{\pi(k)=j}^{AB} \sum_{\tau(k)=i}^{4} r_{Nk} \tilde{\beta}_N\right) \right] \right\}, \quad (88)$$

where $\sum_{\{r_{Nk}\}}$ means summing all the indices $r_{Nk}$, with $1 \le N \le 4$ and with $1 \le k \le AB$, from zero to infinity, where $r_{Nk}$ refers to the number of uncut Pomerons of type N of nucleon-nucleon pair k. It is useful for the discussion to consider “Pomeron types” N, although they are not physical objects, just coming from the parametrization in Eq. (84). I use $\tilde{\beta}_N = \beta_N + 1$, and a function g defined as

$$g(z) = \frac{\Gamma(1 + \alpha_{\text{remn}})}{\Gamma(1 + \alpha_{\text{remn}} + z)}. \quad (89)$$

The result in Eq. (88) is remarkable in the sense that all the integrations could be done analytically, expressed in terms of gamma functions $\Gamma$.

One cannot simplify Eq. (88) any further, without approximations. Looking at this expression, one sees that the two terms $g(\cdots)$ “make the problem”. Without them, the sums could be done, and one would get products of exponentials. But of course, one cannot simply drop these terms. However, if the g functions would have a particular property, namely, a factorization such as

$$g\left(\sum_{\lambda} \tilde{\beta}_{\lambda}\right) = c_1 \prod_{\lambda} c_2 g(c_3 \tilde{\beta}_{\lambda}), \quad (90)$$

for given coefficients $\tilde{\beta}_{\lambda}$ (arbitrary, but in practice of order unity) and three parameters $c_{\mu}$, then the expression in Eq. (88) can be factorized such that the infinite sums are finally just a product of exponentials, and as proven in Appendix B [see Eq. (B64)], and one gets

$$W_{AB} = \prod_{i=1}^{A} c_1(x_i^+)^{\alpha_{\text{remn}}} \prod_{j=1}^{B} c_1(x_j^-)^{\alpha_{\text{remn}}} \times \prod_{k=1}^{AB} \exp(-\tilde{G}(x_{\pi(k)}^+, x_{\tau(k)})), \quad (91)$$

for any choice of the parameters $c_{\mu}$, where $\tilde{G}$ is given as

$$\tilde{G}(x) = \sum_{N=1}^{4} \tilde{\alpha}_{N} x^{\tilde{\beta}_{N}}, \quad (92)$$

<figure><img src="image_25.png" /><figcaption>FIG. 23. Comparing $g(\sum_{\lambda} \tilde{\beta}_{\lambda})$ and $c_1 \prod_{\lambda} c_2 g(c_3 \tilde{\beta}_{\lambda})$ for randomly created sequences $\tilde{\beta}_{\lambda}$, for the parameter choice (A) $c_1 = 2$, $c_2 = 0.65$, and $c_3 = 1$.</figcaption></figure>

with

$$\tilde{\alpha}_N = \alpha_N \left( \frac{\Gamma(\tilde{\beta}_N) c_2 \Gamma(1 + \alpha_{\text{remn}})}{\Gamma(1 + \alpha_{\text{remn}} + c_3 \tilde{\beta}_N)} \right)^2, \quad (93)$$

$$\tilde{\beta}_N = \beta_N + 1. \quad (94)$$

What is very nice: $W_{AB}$ as given in Eq. (91) is strictly non-negative, and this was the missing piece which allows a probabilistic interpretation of the formulas [see the discussion around Eq. (81)], which is extremely important for any Monte Carlo application.

But the question arises: Is the property in Eq. (90) really true? To answer this question, I make some simulations, where I take $\alpha_{\text{remn}} = 1$ (currently used), so one has $g(z) = \Gamma(2)/\Gamma(2+z)$. In principle, Eq. (90) should be valid for any choice of sequences $\tilde{\beta}_{\lambda}$. So I generate randomly integer numbers $\lambda_{\text{MAX}}$ between zero and some upper limit (20), and then I generate a sequence $\tilde{\beta}_1, \tilde{\beta}_2, \dots, \tilde{\beta}_{\lambda_{\text{MAX}}}$ of $\lambda_{\text{MAX}}$ uniformly distributed random numbers. For each sequence, I compute $g(\sum_{\lambda} \tilde{\beta}_{\lambda})$ and $c_1 \prod_{\lambda} c_2 g(c_3 \tilde{\beta}_{\lambda})$, and then I compare the two. In Fig. 23, I plot the two quantities as a function of $\sum_{\lambda} \tilde{\beta}_{\lambda}$, for the parameter choice (A) $c_1 = 2$, $c_2 = 0.65$, and $c_3 = 1$. In Fig. 24, I show the comparison for a second parameter choice (B) $c_1 = 2$, $c_2 = 1.5$, and $c_3 = 2.8$. Neither of the two choices is perfect (one will understand later that a perfect choice cannot exist). Choice A is somewhat better for small $\sum_{\lambda} \tilde{\beta}_{\lambda}$, but then deviates for large values, whereas choice B gives an overall good description, but with bigger fluctuations. In the following, I will take choice A, mainly because B has only been discovered recently, and all EPOS4 simulations are based on A. A g-factorization similar to choice A was employed in Ref. [26].

The next obvious question is: How does the g-factorization [see Eq. (90)] compare to the exact result, concerning $W_{AB}$? By exact result I mean a numerical calculation of Eq. (88), which at least for pp scattering ($A = B = 1$) can be easily done, since the infinite sums converge fast. One recalls that $W_{11}$ depends also on the impact parameter, and on the collision energy.
<figure><img src="image_26.png" /><figcaption>FIG. 24. Same as Fig. 23, but for the parameter choice (B) $c_1 = 2$, $c_2 = 1.5$, and $c_3 = 2.8$.</figcaption></figure>

In Fig. 25, I show $W_{11}(\sqrt{x}, \sqrt{x})$ divided by the trivial factor $x^{\alpha_{\text{remn}}}$, as a function of $x = x^{+}x^{-}$ (so I use here $x^{+} = x^{-}$), for pp scattering at 13 TeV, for the impact parameters $b = 3.5, 2.5$, and $1.5$ fm. The solid curves refer to the exact results, and the dotted ones to the calculations based on g-factorization [see Eq. (90)]. The two methods give very similar results, and they look reasonable: at very large $b$, one gets for $x = 1$ unity, as it should be, since $W_{11}(1, 1)$ is interpreted as the probability of no interaction. With decreasing $b$, the probability of no interaction decreases.

The situation changes completely when one considers small values of the impact parameter $b$. In Fig. 26, I show $W_{11}(\sqrt{x}, \sqrt{x})$ divided by $x^{\alpha_{\text{remn}}}$ as a function of $x = x^{+}x^{-}$, for impact parameters $b = 0.5$ and $0$ fm. Again the solid curves refer to the exact results, and the dotted ones to the calculations based on g-factorization. The results are amazing: Not only are the curves from the two methods substantially

<figure><img src="image_27.png" /><figcaption>FIG. 25. $W_{11}(\sqrt{x}, \sqrt{x})/x^{\alpha_{\text{remn}}}$ as a function of $x = x^{+}x^{-}$, for pp scattering at 13 TeV, for the impact parameters $b = 3.5, 2.5$, and $1.5$ fm. The solid curves refer to the exact results, and the dotted ones to the calculations based on g-factorization [Eq. (90)].</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>FIG. 26. As Fig. 25, but for impact parameters $b = 0.5$ and $0$ fm.</figcaption></figure>

different, but $W_{11}$ even becomes negative at large $x$ for the exact calculations, whereas the g-factorized results converge towards zero, as it should be.

I will now focus on the large-$x$ behavior, and consider $1 - W_{11}$ for $x^{+} = x^{-} = 1$. As discussed earlier, $1 - W_{11}(1, 1)$ is interpreted as the probability of an inelastic interaction (summed over all possible numbers of exchanged Pomerons). I plot $1 - W_{11}(1, 1)$ as a function of the impact parameter $b$, as shown in Fig. 27, where the solid curve refers to the exact results, and the dotted one to the calculations based on g-factorization [see Eq. (90)]. The two curves are very close for large impact parameters but differ for small $b$ (compatible with Figs. 25 and 26). For small $b$, the curve based on g-factorization approaches unity, as it should, but the exact calculation (solid curve) exceeds unity.

D. Conclusion

So one observes some “unphysical behavior” in the sense that $1 - W_{11}(1, 1)$, interpreted as the “probability of an inelastic scattering”, exceeds unity at small impact parameter $b$, or the other way around, the quantity $W_{11}(1, 1)$, interpreted as “probability of no interaction”, becomes negative. This is

<figure><img src="image_29.png" /><figcaption>FIG. 27. $1 - W_{11}(1, 1)$ as a function of $b$, for pp scattering at 13 TeV. The solid curve refers to the exact result, and the dotted one to the calculations based on g-factorization [see Eq. (90)].</figcaption></figure>
what I refer to as the “sign problem”. Imposing g-factorization [see Eq. (90)] provides an excellent approximation to the exact solution at large impact parameter, but deviates at small b by correcting the sign problem. Therefore, it is a “regularization” (which fixes a problem) rather than a simple approximation.

Let me summarize this section (containing a large amount of information):

(i) One generalizes the Gribov-Regge (GR) S-matrix approach, to include energy-momentum sharing (GR$^+$).

(ii) One writes an expression for the inelastic nucleus-nucleus cross section $\sigma_{in}^{AB}$ summing over all possible numbers of cut and uncut Pomerons for each nucleon-nucleon pair, integrating over all Pomeron momenta, imposing strict energy-momentum sharing. The elementary object representing a single Pomeron is G (cut T-matrix element).

(iii) One proves the validity of the AGK theorem, which is an extremely important issue, since it allows one to deduce binary scaling in AA scattering and factorization in pp scattering.

(iv) From the expression for $\sigma_{in}^{AB}$, one deduces an equation $1 = \sum_{\{m_k\}} \int dX_{AB} P$, which allows the interpretation of

$$P = \prod_{k=1}^{AB} \left[ \frac{1}{m_k!} \prod_{\mu=1}^{m_k} G_{k\mu} \right] \times W_{AB}(\{x_i^+\}, \{x_j^-\}) \quad (95)$$

to be the probability of the “configuration” $K = \{\{m_k\}, \{x_{k\mu}^{\pm}\}\}$, representing $m_k$ cut Pomerons per pair k, with light-cone momentum fractions $x_{k\mu}^{\pm}$. This is very important since it allows Monte Carlo applications perfectly compatible with the theoretical S-matrix-based formulas. But one has to prove $W_{AB} \ge 0$.

(v) The computation of $W_{AB}$ is challenging due to a very high-dimensional nonseparable integration. Introducing a particular (and well-justified) parametrization of the $x^{\pm}$ dependence of G as $\sum_{N=1}^{N_{par}} \alpha_N (x^+ x^-)^{\beta_N}$, the integration can be done, and one deduces an analytic formula for $W_{AB}$.

(vi) $W_{AB}$ is still a complex expression, given as $\sum_{\{r_{Nk}\}} \dots$, which means summing indices $r_{Nk}$ (from 0 to $\infty$), with $1 \le N \le 4$ and with $1 \le k \le AB$. Doing a “small” manipulation (let me call it “regularization”), the sums can be done, one gets a simple product of exponentials, and, most importantly, $W_{AB} \ge 0$, as it should be.

(vii) For pp scattering, the regularized results and the exact results can be shown to be identical, for small values of the remnant momentum fractions $x^{\pm}$ (1 minus the sum of the Pomeron momentum fractions) and always for large impact parameters.

(viii) But for small impact parameter, for large $x^{\pm}$, and in particular for $x^{\pm} = 1$, the value of $W_{11}$ becomes negative. And $1 - W_{11}$, to be interpreted as the probability of an inelastic interaction, exceeds unity. Here, the regularized results deviate from the ex-

act ones, since for the former one gets strictly 1 – $W_{11} = 1$ for small impact parameters.

Imposing this regularization seems to solve the problem of negative probabilities, but also indicates that there still is something missing, and the theoretical framework is not yet complete. This will be discussed in the next section.

## VII. AGK VIOLATION AND DEFORMATION FUNCTIONS IN THE REGULARIZED THEORY

In the last section, I developed the parallel scattering S-matrix approach for heavy-ion scattering, with energy-momentum sharing, the latter being crucial for any realistic application. I could prove the validity of the AGK theorem, which is an essential property, the condition for binary scaling and factorization. An expression [see Eq. (81)]

$$P(K) = \prod_{k=1}^{AB} \left[ \frac{1}{m_k!} \prod_{\mu=1}^{m_k} G_{k\mu} \right] \times W_{AB}(\{x_i^+\}, \{x_j^-\}) \quad (96)$$

could be derived, to be interpreted as a probability law for particular multi-Pomeron configurations $K = \{\{m_k\}, \{x_{k\mu}^{\pm}\}\}$. The $W_{AB}$ term turned out to be negative, for certain values of the arguments, which is unphysical. A “regularization” allowed one to solve this problem, and all $W_{AB}$ (and the P’s) are finally non-negative. In addition, this regularization changes $\sigma_{in}^{AB}$ only little compared to the exact calculation (at least for $A=B=1$).

### A. AGK theorem in the regularized approach

One has the impression that after the regularization the problem is solved, but before coming to this conclusion, at least the (very important) AGK theorem has to be checked.

Based on the (reestablished) probability interpretation of the $P(K)$ expressions, one may generate corresponding configurations, and assuming

$$G = G_{QCD} \quad (97)$$

with $G_{QCD}$ being discussed in detail in Ref. [7] [see also Fig. 19 and Eq. (48)], one may generate partons, and employing the EPOS method to translate partons into strings and then hadrons (as well discussed in Ref. [7]), one finally produces a final state of hadrons. This allows the computation of inclusive cross sections with respect to the transverse momenta of partons or of hadrons. But I will first consider the inclusive cross section with respect to the light-cone momentum fractions $x^{\pm}$ of the Pomerons, obtained from $P(K)$ defined in Eq. (96) by summing over all $\{m_k\}$ and an integrating over $b_{AB}$ and $\{x_{k\mu}^{\pm}\}$, and by adding a factor

$$\sum_{k'=1}^{AB} \sum_{\mu'=1}^{m_{k'}} \delta(x^{+} - x_{k'\mu'}^{+}) \delta(x^{-} - x_{k'\mu'}^{-}) dx^{+} dx^{-} \quad (98)$$

[see Eqs. (54) and (55)]. One gets

$$\frac{d^2 \sigma_{\text{incl}}^{AB}}{dx^+ dx^-} = \sum_{k'=1}^{AB} \sum_{\{m_k\} \neq 0} \sum_{\mu'=1}^{m_{k'}} \int db_{AB} \int dX_{AB} \\ \times \{P(K) \delta(x^+ - x_{k'\mu'}^+) \delta(x^- - x_{k'\mu'}^-)\}, \quad (99)$$
<figure><img src="image_30.png" /><figcaption>FIG. 28. $R_{AGK}$ as a function of the $p_t$ of partons for minimum bias PbPb and pp scatterings at 5.02 TeV.</figcaption></figure>

where $\int dX_{AB}$ represents the integration over $\{x_{k\mu}^{\pm}\}$ [see Eq. (72)], and $\int db_{AB}$, defined in Eqs. (40) and (41), integrates over the impact parameter and the transverse coordinates $b_i^A$ and $b_j^B$ of projectile and target nucleons. The weight $P$ is given in Eq. (96), with $W_{AB}$ from Eq. (91), and with the arguments $\{x_i^+\}$, $\{x_j^-\}$ of $W_{AB}$ from Eq. (80). Equation (99) can be evaluated numerically, based on Monte Carlo simulations using Markov chains [5].

Let me make some checks. The “regularization” affects $\sigma_{in}^{AB}$ very little, but this is not necessarily true for other quantities, and most important is the validity of AGK theorem. Does it still hold in the regularized theory? Validity of the AGK theorem means that inclusive cross sections in the full nuclear $(A+B)$ scattering is $AB$ times the result for a single Pomeron,

$$ \frac{d\sigma_{\text{incl}}^{AB}}{dx^+dx^-} = AB \times \frac{d\sigma_{\text{incl}}^{\text{single Pom}}}{dx^+dx^-} \quad (100) $$

[see Eq. (64) and Fig. 21]. But since for a given $x^{\pm}$ distribution, and knowing the structure of a Pomeron ($G = G_{QCD}$), one obtains from $d\sigma_{incl}/dx^+dx^-$ in a unique fashion $d\sigma_{incl}/dq$ for any single-particle variable $q$, like for example the transverse momentum $p_t$ of partons. So validity of the AGK theorem implies

$$ \frac{d\sigma_{\text{incl}}^{AB}}{dp_t} = AB \times \frac{d\sigma_{\text{incl}}^{\text{single Pom}}}{dp_t}. \quad (101) $$

And this is easy to check: One computes $d\sigma_{incl}^{\text{single Pom}}/dp_t$ based on $G = G_{QCD}$ (see Ref. [7]). Then one runs the full Monte Carlo simulation and computes $d\sigma_{incl}^{AB}/dp_t$ for minimum bias scatterings. Finally, one computes the ratio

$$ R_{\text{AGK}}(p_t) = \frac{d\sigma_{\text{incl}}^{AB}}{dp_t} \Big/ \left\{ AB \times \frac{d\sigma_{\text{incl}}^{\text{single Pom}}}{dp_t} \right\}. \quad (102) $$

Validity of the AGK theorem then simply means $R_{AGK}(p_t) = 1$; everything else means violation of AGK. In Fig. 28, I plot $R_{AGK}$ for PbPb and for pp at 5.02 TeV, and one sees that in both cases one is far from unity, so AGK is badly violated, and the violation is much stronger for PbPb compared to pp, and it increases with $p_t$. One should keep in mind that one has for the moment just primary scatterings, no final state effects, and

<figure><img src="image_31.png" /><figcaption>FIG. 29. $R_{PbPb}$ as a function of the $p_t$ of partons for minimum bias PbPb scatterings at 5.02 TeV.</figcaption></figure>

one expects therefore at least at large $p_t$ a value of $R_{AGK}$ equal to unity. The commonly used “nuclear modification factor” $R_{PbPb}$ is given as the ratio of the two $R_{AGK}$ curves,

$$ R_{\text{PbPb}} = R_{\text{AGK}}^{\text{PbPb}}/R_{\text{AGK}}^{\text{PP}}. \quad (103) $$

shown in Fig. 29. The result is way below unity, around 0.4 at 50 GeV/c, and it decreases with $p_t$. As a side remark: Usually $R_{PbPb}$ is considered for given centrality classes, and the ratio is computed as $d\sigma_{incl}^{AB}/dp_t/\{N_{coll} \times d\sigma_{incl}^{\text{single Pom}}/dp_t\}$, but $N_{coll}$ is not a measurable quantity, whereas rescaling with $AB$ is well defined. The PbPb result at 5.02 TeV is just an example; similar results can be found for many different systems at many different energies. The effect increases with increasing nuclear mass and with energy.

## B. Deformation functions

So one has clearly a problem. $R_{PbPb}$ is not as it should be, and this is a consequence of the violation of the AGK theorem, with $R_{AGK} < 1$, actually substantially less than unity. As mentioned above, the $p_t$ distributions can be derived from $x^{\pm}$ distributions $d\sigma_{incl}^{AB}/dx^+dx^-$. Validity of AGK would mean

$$ R_{\text{AGK}}(x^+, x^-) = \frac{d\sigma_{\text{incl}}^{AB}}{dx^+dx^-} \Big/ \left\{ AB \times \frac{d\sigma_{\text{incl}}^{\text{single Pom}}}{dx^+dx^-} \right\} = 1, \quad (104) $$

and this is what one needs, to have $R_{AGK}(p_t) = 1$. The latter is expected for large $p_t$, which corresponds to large values of $x^+x^-$. In the following, I will investigate how the $x^{\pm}$ dependence of $d\sigma_{incl}^{AB}/dx^+dx^-$ changes (gets “deformed”) with respect to the reference curve for a single Pomeron, so I define the ratio of the normalized distributions,

$$ R_{\text{deform}}(x^+, x^-) = \frac{1}{\sigma_{\text{incl}}^{AB}} \frac{d\sigma_{\text{incl}}^{AB}}{dx^+ dx^-} \Big/ \left\{ \frac{1}{d\sigma_{\text{incl}}^{\text{single Pom}}} \frac{d\sigma_{\text{incl}}^{\text{single Pom}}}{dx^+ dx^-} \right\}, \quad (105) $$

referred to as the “deformation function”. It is useful to define the Pomeron squared energy fraction $x_{PE}$ and the Pomeron
<figure><img src="image_32.png" /><figcaption>FIG. 30. $R_{\text{deform}}(x_{\text{PE}})$ as a function of $x_{\text{PE}}$ for PbPb (central, 0-5%) and pp (event class 16-20 Pomerons) scatterings at 5.02 TeV.</figcaption></figure>

rapidity $y_{\text{PE}}$ as

$$x_{\text{PE}} = x^{+}x^{-} = \frac{M_{\text{Pom}}^{2}}{s}, \quad y_{\text{PE}} = 0.5 \ln \frac{x^{+}}{x^{-}}. \qquad (106)$$

The quantity $M_{\text{Pom}}$ is the transverse mass of the Pomeron. The corresponding distributions (for $\sigma_{\text{incl}}^{\text{AB}}$ and $\sigma_{\text{incl}}^{\text{single Pom}}$) are

$$\frac{d\sigma_{\text{incl}}}{dx_{\text{PE}}dy_{\text{PE}}} = J \times \frac{d\sigma_{\text{incl}}}{dx^{+}dx^{-}}, \qquad (107)$$

with $J$ being the corresponding Jacobian determinant. This defines $R_{\text{deform}}(x_{\text{PE}}, y_{\text{PE}})$ as the corresponding ratio. The Pomeron rapidity $y_{\text{PE}}$ turns out to be close to zero, so one usually integrates over $y_{\text{PE}}$, which gives (again for $\sigma_{\text{incl}}^{\text{AB}}$ and $\sigma_{\text{incl}}^{\text{single Pom}}$) the distribution $d\sigma_{\text{incl}}/dx_{\text{PE}}$, and the ratio of the $A+B$ cross section over the single Pomeron one provides, $R_{\text{deform}}(x_{\text{PE}})$.

In Fig. 30, I show results for $R_{\text{deform}}(x_{\text{PE}})$ in PbPb (central, 0-5%) and pp (event class 16-20 Pomerons) scatterings at 5.02 TeV. These curves should be unity, but they are far from that. At large $x_{\text{PE}}$, the ratios are much smaller than unity, whereas they exceed unity at small $x_{\text{PE}}$. Many simulations were done to compute $R_{\text{deform}}(x_{\text{PE}})$ for different energies and systems; see Fig. 31, where I show the results for PbPb at 5.02 TeV, XeXe at 5.44 TeV, AuAu at 0.2 TeV (centralities 0-5%, 5-10%, 10-20%, 30-40%, 50-60%, 70-80%, 80-90%, 90-100%), and pp at 13 TeV (number of Pomerons: 51-60, 31-40, 21-25, 16-20, 10-15, 5-9, 2-4, 1). In all cases, one sees a suppression at large $x_{\text{PE}}$, and the effect gets bigger with increasing event activity. The effect is actually biggest in pp, for rare events with very large Pomeron numbers (the average number is around two).

The current status is as follows: Whereas originally one could prove the validity of AGK, one needed to introduce a “regularization” to avoid negative probabilities, which in turn ruined this very important AGK property: for all systems, energies, and centralities, one gets deformed $x_{\text{PE}}$ distributions with respect to the reference distribution of a single Pomeron, which leads to AGK violation with respect to $x_{\text{PE}}$, which as a consequence violates AGK with respect to $p_{\text{t}}$, which finally leads to a violation of binary scaling in AA scattering. It also

<figure><img src="image_33.png" /><figcaption>FIG. 31. $R_{\text{deform}}(x_{\text{PE}})$ as a function of $x_{\text{PE}}$ for different systems, different energies (in TeV), different event classes “Mi” (M1 = highest and M8 = lowest multiplicity). See text.</figcaption></figure>

leads to a violation of factorization in pp, but this requires in addition to consider the internal structure of G.

So at the heart of all these problems are the “deformed” $x^{\pm}$ or $x_{\text{PE}}$ distributions. In the following, I will try to better understand and eventually parametrize these deformations, as a first step towards a solution. Let me consider in AA collisions (including pp as a special case) a particular multi-Pomeron configurations $K = \{\{m_k\}, \{x_{k\mu}^{\pm}\}\}$, with $m_k$ cut Pomerons per nucleon-nucleon pair $k$, with Pomeron light-cone momentum fractions $x_{k\mu}^{\pm}$. And let me consider a particular Pomeron, connected to projectile nucleon $i$ and target nucleon $j$ (see Fig. 32). In the configuration $K$, there might be other Pomerons, connected to one (or both) of these nucleons. The corresponding Pomeron-nucleon connections are marked as

<figure><img src="image_34.png" /><figcaption>FIG. 32. A Pomeron connected to projectile nucleon $i$ and target nucleon $j$, together with other Pomerons connected to one (or both) of these nucleons.</figcaption></figure>
red and blue dots. It is obvious that the additional Pomerons connected to the same nucleons i and j compete with each other; they have to share the initial energy-momentum of the two nucleons. The more Pomerons are connected, the less energy is available for one particular Pomeron. To quantify this statement, I define the “connection number”

$$N_{\text{conn}} = \frac{N_{\text{P}} + N_{\text{T}}}{2}, \qquad (108)$$

with $N_{\text{P}}$ being the number of Pomerons connected to i, and with $N_{\text{T}}$ being the number of Pomerons connected to j (the variable $N_{\text{conn}}$ corresponds to half of the number of red and blue points in Fig. 32).

In the following, I will discuss the effect of energy sharing related to the connection number. One wants to understand the dependence of the $x^{\pm}$ distributions on the connection number $N_{\text{conn}}$, so rather than Eq. (99), I will compute

$$\frac{d^2\sigma_{\text{incl}}^{\text{AB}}(N_{\text{conn}})}{dx^+dx^-} = \sum_{k'=1}^{AB} \sum_{\{m_k\}\neq0} \sum_{\mu'=1}^{m_{k'}} \delta_{N_{\text{conn}}(k',\mu')} \int db_{\text{AB}} \int dX_{\text{AB}} \\ \times \{P(K)\delta(x^+ - x_{k'\mu'}^+)\delta(x^- - x_{k'\mu'}^-)\}, \quad (109)$$

where $\delta_a^b$ is the Kronecker delta. So I only consider Pomerons $k'$, $\mu'$ with connection number $N_{\text{conn}}(k', \mu')$ equal to $N_{\text{conn}}$. Equation (109) can be easily evaluated numerically, based on Monte Carlo simulations.

Although the above method to compute the inclusive $x^{\pm}$ cross sections for given $N_{\text{conn}}$ is perfectly doable, I proceed somewhat differently. I define event classes (centrality in AA scattering, number of Pomerons in pp), and then compute the average values of the cross sections $d^2\sigma_{\text{incl}}^{\text{AB}}/dx_{\text{PE}}$ and $N_{\text{conn}}$ per event class. This allows to compute and then investigate $R_{\text{deform}}(x_{\text{PE}})$ for particular event classes, with the associated mean $N_{\text{conn}}$, in different collision systems at various energies. Studying the corresponding plots, one realizes that all curves $R_{\text{deform}}(x_{\text{PE}})$ can be parametrized in the following way. One first defines

$$\begin{aligned} w &= \ln(R_{\text{deform}})/20, \\ u &= -\ln(x_{\text{PE}})/20, \end{aligned} \qquad (110)$$

and then uses

$$\begin{aligned} v &= |u - a_1|/0.5, \\ w &= a_2 - v^{a_3} a_4, \\ \text{if } (u > a_1) \quad w &= \max(0, w), \\ w &= \max(a_5, w). \end{aligned} \qquad (111)$$

In Fig. 33, I plot $R_{\text{deform}}(x_{\text{PE}})$ as a function of $x_{\text{PE}}$, for PbPb at 5.02 TeV, for the centralities M1 to M8, i.e., 0–5%, 5–10%, 10–20%, 30–40%, 50–60%, 70–80%, 80–90%, 90–100%. The red solid curves refer to the simulated results, the black dotted ones to the parameterized curves.

In Table II, I report the values of $\langle N_{\text{conn}} \rangle$ for the different centrality classes, as well as the values for the parameters. Somewhat unexpectedly, the deformations do not change much with centrality in the range 0–20% (M1–M3), explained by a correspondingly little change in $\langle N_{\text{conn}} \rangle$. As a side remark, the average number of Pomerons or the average number of

<figure><img src="image_35.png" /><figcaption>FIG. 33. $R_{\text{deform}}(x_{\text{PE}})$ as a function of $x_{\text{PE}}$ for PbPb at 5.02 TeV, for the centralities M1 to M8. See text.</figcaption></figure>

collisions do change considerably from M3 to M1, but not the connection number, and it is the latter that counts concerning energy sharing.

Tables like Table II are considered to define the $N_{\text{conn}}$ dependence of $R_{\text{deform}}(x_{\text{PE}})$ for a given system at a given energy, first for the $N_{\text{conn}}$ values in the table, but also for arbitrary values via interpolation or extrapolation.

At this point one should discuss some important details: In order to compute $R_{\text{deform}}(x_{\text{PE}})$, one needs to do simulations based on the probability law in Eq. (96), and to do so one needs to know the single Pomeron function $G(x^+, x^-)$, which Eq. (96) is based upon. As discussed in Sec. VI, and actually needed to derive the final formulas for the laws, one

<table><caption>TABLE II. $R_{\text{deform}}$ parameters for PbPb at 5.02 TeV. See text.</caption><thead><tr><th>Class</th><th>⟨N<sub>conn</sub>⟩</th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th><th>a<sub>4</sub></th><th>a<sub>5</sub></th></tr></thead><tbody><tr><td>M1</td><td>7.69</td><td>0.49</td><td>0.076</td><td>2.28</td><td>0.38</td><td>-0.15</td></tr><tr><td>M2</td><td>7.68</td><td>0.49</td><td>0.076</td><td>2.28</td><td>0.38</td><td>-0.15</td></tr><tr><td>M3</td><td>7.63</td><td>0.49</td><td>0.076</td><td>2.28</td><td>0.38</td><td>-0.14</td></tr><tr><td>M4</td><td>7.25</td><td>0.49</td><td>0.075</td><td>2.26</td><td>0.36</td><td>-0.13</td></tr><tr><td>M5</td><td>6.12</td><td>0.47</td><td>0.067</td><td>2.12</td><td>0.33</td><td>-0.12</td></tr><tr><td>M6</td><td>4.18</td><td>0.45</td><td>0.054</td><td>2.08</td><td>0.26</td><td>-0.11</td></tr><tr><td>M7</td><td>3.23</td><td>0.44</td><td>0.046</td><td>2.00</td><td>0.22</td><td>-0.10</td></tr><tr><td>M8</td><td>2.21</td><td>0.42</td><td>0.040</td><td>1.90</td><td>0.20</td><td>-0.07</td></tr></tbody></table>
<figure><img src="image_36.png" /><figcaption>FIG. 34. $R_{deform}(x_{PE})$ as a function of $x_{PE}$ for XeXe at 5.44 TeV, for the centralities M1 to M8. See text.</figcaption></figure>

parametrizes the $x^{\pm}$ dependence of $G$ as [see Eq. (84)]

$$G(x^{+}, x^{-}, s, b) = \sum_{N=1}^{N_{\text{par}}} \alpha_{N} (x^{+}x^{-})^{\beta_{N}}, \quad (112)$$

inspired by asymptotic form of T-matrices (see Appendix A). I argued that this form provides an excellent fit to numerically computed expressions of $G_{QCD}$, based on pQCD, discussed in detail in Ref. [7]. However, in the present work, one does not require $G = G_{QCD}$ anymore—it is just a starting point—and one is free to change parameters to get an appropriate behavior of elementary quantities as total and elastic cross sections. At this point, one does not need (yet) to specify the internal structure of $G$ (the relation between $G$ and $G_{QCD}$). One only needs to know the parametric form [Eq. (112)] of $G$; this is enough to compute the deformations.

I now continue to show the results concerning the simulation and parametrization of $R_{deform}(x_{PE})$ for other systems.

In Fig. 34, I present results for XeXe at 5.44 TeV, for the centralities M1 to M8, i.e., 0–5%, 5–10%, 10–20%, 30–40%, 50–60%, 70–80%, 80–90%, 90–100%. The red solid curves refer to the simulated results, the black dotted ones to the parametrized curves. These results are very similar to the 5.02 TeV PbPb ones, and also the table corresponding to Table II is almost identical. Also here, the deformations (and the $N_{\text{conn}}$ values) vary little with centrality.

In Fig. 35, I show results for AuAu at 200 GeV, for the centralities M1 to M8, i.e., 0–5%, 5–10%, 10–20%, 30–40%, 50–60%, 70–80%, 80–90%, 90–100%. As already observed at higher energies, the deformations are as well varying only

<figure><img src="image_37.png" /><figcaption>FIG. 35. $R_{deform}(x_{PE})$ as a function of $x_{PE}$ for AuAu at 200 GeV, for the centralities M1 to M8. See text.</figcaption></figure>

slowly with centrality. Here, the values for $\langle N_{\text{conn}} \rangle$ evolve between 1 (M8) and 5 (M1). Obviously the covered $x_{PE}$ range is smaller because the collision energy is smaller.

In Fig. 36, I finally show results for $pp$ at 13 TeV, where the different event classes M1–M8 refer to ranges in the number of Pomerons: 51–60, 31–40, 21–25, 16–20, 10–15, 5–9, 2–4, 1. Contrary to the heavy-ion results, the values for $\langle N_{\text{conn}} \rangle$ vary strongly: they go from unity up to 28.5. And correspondingly the deformations are huge. This will have important consequences, concerning measurable observables.

### C. Conclusions

Let me summarize this section.

(i) **Reminder:** Validity of the AGK theorem means *in-clusive cross sections* of collisions of nuclei of mass numbers A and B are equal to AB times the cross section of a *single* Pomeron exchange. The most “basic” inclusive cross section is the one with respect to the variables $x^{\pm}$ (or $x_{PE}$ and $y_{PE}$); everything else can be derived from this.

(ii) As a first test of the “regularized theory”, I check the validity of AGK. Doing a PbPb simulation at 5.02 TeV, one finds that AGK fails badly, and indeed also binary scaling does not work at all.

(iii) To understand the problem, I study the ratio of inclusive cross sections with respect to $x_{PE}$ for A+B scattering divided by the single Pomerons result, the
<figure><img src="image_38.png" /><figcaption>FIG. 36. $R_{\text{deform}}(x_{\text{PE}})$ as a function of $x_{\text{PE}}$ for pp at 13 TeV, for the event classes M1 to M8. See text.</figcaption></figure>

latter being the reference curve. Both cross sections are normalized. The ratio is called the “deformation function” $R_{\text{deform}}(x_{\text{PE}})$.

(iv) Many simulations are done, for different systems and energies. In all cases, $R_{\text{deform}}(x_{\text{PE}})$ deviates from unity; one sees always a suppression at large $x_{\text{PE}}$, and the effect gets bigger with increasing event activity.

(v) I define for a given Pomeron a quantity called connection number $N_{\text{conn}}$, counting the number of Pomerons connected to the same projectile or target, strongly correlated with the deformations: large $N_{\text{conn}}$ means large deformation (due to energy sharing). The variable $N_{\text{conn}}$ represents the “environment” of a given Pomeron; it states if the Pomeron is isolated or competes with others with respect to energy-momentum sharing.

(vi) I study the dependence of the deformations on $\langle N_{\text{conn}} \rangle$ for different event classes, and one finds simple parametrizations for $R_{\text{deform}}(x_{\text{PE}})$ in terms of five $N_{\text{conn}}$-dependent parameters. This allows one to define $R_{\text{deform}}(x_{\text{PE}})$ for arbitrary $N_{\text{conn}}$ values via interpolation and extrapolation.

So deformations $R_{\text{deform}}(x_{\text{PE}}) \neq 1$ represent the problem (AGK violation), but at least one understands the origin of the problem (energy-momentum sharing), and one is able to parametrize it. This will be useful for the next step: the solution of the problem.

<figure><img src="image_39.png" /><figcaption>FIG. 37. A single cut Pomeron.</figcaption></figure>

## VIII. THE SOLUTION VIA SATURATION AND GENERALIZED AGK THEOREM

One has seen in Sec. VII that the validity of the AGK theorem is badly violated, with respect to the variable $x_{\text{PE}}$. There is nothing one can do about it; it is an unavoidable consequence of energy-momentum sharing, and the effect can be quantified in terms deformation functions $R_{\text{deform}}(x^+, x^-)$ [in practice $R_{\text{deform}}(x_{\text{PE}})$], it being the ratio of $A+B$ inclusive cross sections over the single Pomeron one, using normalized distributions. Here validity of AGK would mean $R_{\text{deform}}(x^+, x^-) = 1$, which is not at all observed; one has always a suppression at large $x_{\text{PE}} = x^+x^-$. This invalidity of AGK has important consequences, like the violation of binary scaling in nuclear collisions.

Let me recall that by AGK theorem I mean some extension of the original one, in the sense that an inclusive cross section with respect to some variable $q$ in $A+B$ scattering is equal to $AB$ times the corresponding single Pomeron cross section.

### A. The role of saturation

But does one really need AGK for inclusive cross sections with respect to the variables $x^\pm$ (or $x_{\text{PE}}$)? Not necessarily: one expects binary scaling to be valid not always, but only for high-$p_t$ processes, and the variable $p_t$ is related to $x^\pm$, but how precisely depends on the internal structure of $G$. One remembers that the single cut Pomeron $G$ is the fundamental building block of the multiple scattering formalism; in the graphs I use simply a vertical line (in cyan) as shown in Fig. 37, where the vertical black dashed line represents the cut. All the multiple scattering formulas depend on $G$, but not on the internal structure. Only when it comes to statements concerning, for example, $p_t$ of produced partons does one need to specify how $G$ is expressed in terms of QCD diagrams. So far I am assuming

$$G = G_{\text{QCD}}, \qquad (113)$$

where $G_{\text{QCD}}$ is essentially a DGLAP parton ladder, discussed in detail in Ref. [7] [see also Fig. 19 and the discussion before Eq. (48)]. But the assumption in Eq. (113) is obviously wrong, because it leads to a strong violation of binary scaling at large $p_t$, as shown in Sec. VII.

There is another serious problem with Eq. (113): As discussed in detail in Ref. [7], the essential part of $G_{\text{QCD}}$ is a cut parton ladder, based on DGLAP parton evolutions. But this is certainly not the full story: With increasing energy, partons with very small momentum fractions $x \ll 1$ become increas-
<figure><img src="image_40.png" /><figcaption>FIG. 38. Cut diagram (a) with nonlinear effects via triple Pomeron vertices and (b) with the nonlinear effects (inside the red circles) “summarized” in the form of saturation scales, which replace these nonlinear parts.</figcaption></figure>

ingly important, since the parton density becomes large, and therefore the linear DGLAP evolution scheme is not valid anymore, and nonlinear evolution takes over, considering explicitly gluon-gluon fusion. These phenomena are known as “saturation” [11–24].

At least for scatterings carrying a large value of $x^+x^-$, one expects “nonlinear effects”, which means that two ladders, which evolve first independently and in parallel, finally fuse. And only after that, the (linear) DGLAP evolution is realized. In the “Pomeron language”, this corresponds to diagrams with triple (and more) Pomeron vertices, as sketched in Fig. 38(a). Such nonlinear effects lead to strong destructive interference at low transverse momentum ($p_t$), which may be summarized in terms of a saturation scale [12,13]. This suggests treating these “saturation phenomena” not explicitly, but by introducing saturation scales as the lower limits of the virtualities for the DGLAP evolutions, as sketched in Fig. 38(b). Also in Refs. [27,28] saturation scales “summarize” nonlinear effects, but based on the assumption that the partons exactly cover the transverse size of the nucleon (or nucleus). In the present paper, the criterion is quite different, as discussed below.

So one has two problems—a wrong association $G = G_{QCD}$ and a missing treatment of saturation—but fortunately, the two problems are connected, and there is an amazingly simple solution based on saturation scales that solves both problems.

## B. The solution

Instead of the “naive” assumption $G = G_{QCD}$, one postulates

$$G(x^+, x^-) = \frac{n}{R_{\text{deform}}(x^+, x^-)} G_{\text{QCD}}(Q_{\text{sat}}^2, x^+, x^-), \quad (114)$$

such that $G$ itself does not depend on the environment. By environment I mean, for a given Pomeron in a given configuration of multiple Pomeron exchanges, the connections of other Pomerons to the same projectile and target nucleons as the given Pomeron. The environment is here quantified in terms of $N_{\text{conn}}$. The quantity $R_{\text{deform}}$ is the $N_{\text{conn}}$-dependent deformation function discussed in Sec. VII, and $n$ is a constant, not

<figure><img src="image_41.png" /><figcaption>FIG. 39. The saturation scale $Q_{\text{sat}}^2$ as a function of $x_{\text{PE}} = x^+x^-$, for several $N_{\text{Pom}}$ event classes (from top to bottom, 10–15, 5–9, and 2–4 Pomerons).</figcaption></figure>

depending on $x_{\text{PE}}$. As discussed earlier, one first parametrizes $G$ as in Eq. (84), with the parameters being fixed by comparing simulation results to elementary experimental data, and then uses Eq. (114) to determine $Q_{\text{sat}}^2$. In this way, $Q_{\text{sat}}^2$ depends on $N_{\text{conn}}$ and on $x^\pm$, as does $R_{\text{deform}}$:

$$R_{\text{deform}} = R_{\text{deform}}(N_{\text{conn}}, x^+, x^-), \quad (115)$$

$$Q_{\text{sat}}^2 = Q_{\text{sat}}^2(N_{\text{conn}}, x^+, x^-), \quad (116)$$

$$\text{but } G \text{ independent of } N_{\text{conn}}. \quad (117)$$

The $N_{\text{conn}}$ dependence of $Q_{\text{sat}}^2$ means that the low virtuality cutoff for the DGLAP evolutions in $G_{\text{QCD}}$ is not a constant, but its value depends on the environment in terms of $N_{\text{conn}}$ and on the energy of the Pomeron. I will refer to this as “dynamical saturation scales”.

The fundamental relation in Eq. (114), together with Eqs. (115)–(117) and in particular the dynamical saturation scales, will solve the problems related to the AGK theorem, as will be discussed in the following sections.

To get some idea about the $x^\pm$ dependence of $Q_{\text{sat}}^2$, obtained from Eq. (114), I show in Fig. 39 results for $pp$ scattering at 7 TeV, for several event classes defined via the number of Pomerons, $N_{\text{Pom}}$. I also indicate the average values of $N_{\text{conn}}$ for the three event classes. One sees that $\langle N_{\text{conn}} \rangle$ varies considerably between the different classes, and so does $Q_{\text{sat}}^2$.

In Fig. 40, I show results for PbPb at 2.76 TeV. Here, the variation of the average $N_{\text{conn}}$ values is quite moderate, and towards central collisions $Q_{\text{sat}}^2$ even “saturates” (no variation anymore). Correspondingly, $Q_{\text{sat}}^2$ varies little from semiperial towards central events, in particular compared to the strong variation of $Q_{\text{sat}}^2$ with the event activity in $pp$ scattering.

## C. Dynamical saturation scales

In order to understand how these dynamical saturation scales work, let me investigate the inclusive cross
<figure><img src="image_42.png" /><figcaption>FIG. 40. The saturation scale $Q_{\text{sat}}^2$ as a function of $x_{\text{PE}}$, for several event classes (from top to bottom 5–10%, 10–20%, 30–40%, 50–60% centrality).</figcaption></figure>

section

$$ \frac{d^2 \sigma_{\text{incl}}^{\text{AB}}(N_{\text{conn}})}{dx^+ dx^-} = \sum_{k'=1}^{AB} \sum_{\{m_k\}\neq0} \sum_{\mu'=1}^{m_{k'}} \delta_{N_{\text{conn}}}^{N_{\text{conn}}(k',\mu')} \int db_{AB} \int dx_{AB} \\ \times \{P(K)\delta(x^+ - x_{k'\mu'}^+) \delta(x^- - x_{k'\mu'}^-)\} \quad (118) $$

[see Eqs. (109) and (40), (41), (72), (81), and (91)]. Here, K refers to a multi-Pomeron configuration $K = \{\{m_k\}, \{x_{k\mu}^{\pm}\}\}$, with $m_k$ cut Pomerons per nucleon-nucleon pair k and with Pomeron light-cone momentum fractions $x_{k\mu}^{\pm}$. In the case of a single Pomeron in pp scattering, the expression simplifies enormously, and one gets

$$ \frac{d\sigma_{\text{incl}}^{\text{single Pom}}}{dx^+ dx^-} = \int d^2 b \, G(x^+, x^-, b) W_{11}((1-x^+)(1-x^-)). \quad (119) $$

Employing Eqs. (114)–(117), using $N_{\text{conn}} = 1$ and $R_{\text{deform}} = 1$, this expression may be written as

$$ \int d^2 b \, G_{\text{QCD}}(Q_{\text{sat}}^2(1, x^+, x^-), x^+, x^-) \\ \times W_{11}((1 - x^+)(1 - x^-)). \quad (120) $$

Since the saturation scale is playing a crucial role in the formalism, I write the above result for a single Pomeron pp scattering as

$$ \frac{d\sigma_{\text{incl}}^{\text{single Pom}}}{dx^+ dx^-} = \frac{d\sigma_{\text{incl}}^{\text{single Pom}}}{dx^+ dx^-} [Q_{\text{sat}}^2(1, x^+, x^-)], \quad (121) $$

with the general definition (for any f)

$$ \frac{d\sigma_{\text{incl}}^{\text{single Pom}}}{dx^+ dx^-} [f(x^+, x^-)] \\ = \int d^2 b \, G_{\text{QCD}} \left( \underbrace{f(x^+, x^-)}_{\text{sat. scale}}, x^+, x^- \right) \\ \times W_{11}((1 - x^+)(1 - x^-)). \quad (122) $$

Let me consider the general case $d^2\sigma_{\text{incl}}^{AB(1)}/dx^+dx^-$, i.e., an A+B scattering with $N_{\text{conn}} = 1$, where one has only

isolated Pomerons, and there is never more than one Pomeron connected to a projectile and a target nucleon, for all nonzero contributions. Concerning $\int dX_{AB}\{\cdots\}$, for given $k'$ and $\mu'$, the $x_{k'\mu'}^\pm$ integration is trivial thanks to the $\delta$ functions, and one gets a factor

$$ G(x^+, x^-, b_{k'}) W_{11}((1 - x^+)(1 - x^-)), \quad (123) $$

multiplied by $P(K')$, where $K'$ is simply the configuration K minus the $\mu'$th Pomeron of pair $k'$ (with actually $\mu' = 1$, always). The quantity W is defined in Eq. (91). Although usually not written explicitly, G depends also on the impact parameter, in this case $b_{k'} = |b + b_i^A - b_{j'}^B|$, where $i' = \pi(k')$ and $j' = \tau(k')$ are the projectile and target nucleons, respectively, associated to pair $k'$. The term $P(K')$ does not depend on $x^+$ and $x^-$, because there are no other Pomerons connected to $i'$ and $j'$. The term $P(K')$ does also not depend on $b_i^A$ and $b_{j'}^B$. Let me define $\int dX_{AB}'$ the integration except for the $x_{k'\mu'}^\pm$ integrations (already done), and $\int db_{AB}'$ the impact parameter integrations except for $d^2b_i^A T_A(b_i^A) \int d^2b_{j'}^B T_B(b_{j'}^B)$. Then the integrations $\int db_{AB}' \int dX_{AB}' P(K')$ provide a factor independent of $x^+, x^-, b_i^A$, and $b_{j'}^B$. So finally, for the case $N_{\text{conn}} = 1$, one finds

$$ \frac{d^2 \sigma_{\text{incl}}^{AB(1)}}{dx^+ dx^-} \propto \sum_{k'=1}^{AB} \sum_{\{m_k\}\neq0} \delta_{N_{\text{conn}}}^{N_{\text{conn}}(k',1)} \\ \times \int d^2 b_i^A T_A(b_i^A) \int d^2 b_{j'}^B T_B(b_{j'}^B) \\ \times \int d^2 b \, G(x^+, x^-, |b + b_i^A - b_{j'}^B|) \\ \times W_{11}((1 - x^+)(1 - x^-)), \quad (124) $$

where the symbol “$\propto$” means proportional. A variable change as $\bar{b}' = \bar{b} + \bar{b}_i^A - \bar{b}_{j'}^B$, and then replacing $b'$ by $b$, and using $\int d^2b_i^A T_A(b_i^A) = 1$ and $\int d^2b_{j'}^B T_B(b_{j'}^B) = 1$, allows one to write

$$ \frac{d^2 \sigma_{\text{incl}}^{AB(1)}}{dx^+ dx^-} \propto \sum_{k'=1}^{AB} \sum_{\{m_k\}\neq0} \delta_{N_{\text{conn}}}^{N_{\text{conn}}(k',1)} \int d^2 b \, G(x^+, x^-, b) \\ \times W_{11}((1 - x^+)(1 - x^-)), \quad (125) $$

which is up to a factor equal to

$$ \int d^2 b \, G(x^+, x^-, b) W_{11}((1 - x^+)(1 - x^-)). \quad (126) $$

Employing Eqs. (114)–(117), using $N_{\text{conn}} = 1$ and $R_{\text{deform}} = 1$, this expression becomes identical to what one found for $d\sigma_{\text{incl}}^{\text{single Pom}}/dx^+dx^-$ [see Eqs. (121) and (122)], and one may write

$$ \frac{d^2 \sigma_{\text{incl}}^{AB(1)}}{dx^+ dx^-} \propto \frac{d\sigma_{\text{incl}}^{\text{single Pom}}}{dx^+ dx^-} [Q_{\text{sat}}^2(1, x^+, x^-)]. \quad (127) $$

So even in nuclear collisions, when one restricts oneself to $N_{\text{conn}} = 1$, one recovers the cross section of a single isolated Pomeron.
In the case of $N_{\text{conn}} > 1$, one has [see Eqs. (105) and (119)]

$$ \frac{d^2 \sigma_{\text{incl}}^{\text{AB}}(N_{\text{conn}})}{dx^+ dx^-} \propto R_{\text{deform}}(N_{\text{conn}}, x^+, x^-) \frac{d \sigma_{\text{incl}}^{\text{single Pom}}}{dx^+ dx^-} \quad (128) $$

$$ = R_{\text{deform}}(N_{\text{conn}}, x^+, x^-) \int d^2 b \, G(x^+, x^-, b) \times W_{11}((1-x^+)(1-x^-)). \quad (129) $$

One may now use Eqs. (114)–(117), so the rhs of Eq. (129) is given as

$$ R_{\text{deform}}(N_{\text{conn}}, x^+, x^-) \int d^2 b \frac{n}{R_{\text{deform}}(N_{\text{conn}}, x^+, x^-)} \times G_{\text{QCD}}(Q_{\text{sat}}^2(N_{\text{conn}}, x^+, x^-), x^+, x^-) \times W_{11}((1-x^+)(1-x^-)). \quad (130) $$

Here, $R_{\text{deform}}(N_{\text{conn}}, x^+, x^-)$ cancels out. Note that this quantity does not depend on $b$. So one finds

$$ \frac{d^2 \sigma_{\text{incl}}^{\text{AB}}(N_{\text{conn}})}{dx^+ dx^-} \propto \int d^2 b \, G_{\text{QCD}}(Q_{\text{sat}}^2(N_{\text{conn}}, x^+, x^-), x^+, x^-) \times W_{11}((1-x^+)(1-x^-)), \quad (131) $$

and using Eq. (122), one gets

$$ \frac{d^2 \sigma_{\text{incl}}^{\text{AB}}(N_{\text{conn}})}{dx^+ dx^-} \propto \frac{d \sigma_{\text{incl}}^{\text{single Pom}}}{dx^+ dx^-} [Q_{\text{sat}}^2(N_{\text{conn}}, x^+, x^-)], \quad (132) $$

where the rhs is again the single Pomeron (pp scattering) expression, just with $Q_{\text{sat}}^2(N_{\text{coll}}, x^+, x^-)$ as saturation scale. The crucial point of Eq. (132) is the fact that, thanks to Eq. (114) and since $G$ does not depend on $N_{\text{conn}}$, the $R_{\text{deform}}$ expressions disappear. So I have shown the following:

(i) For $A+B$ scattering as well as for pp scattering, even with large $N_{\text{conn}}$, the inclusive cross sections with respect to $x^{\pm}$ are always expressed in terms of single Pomeron cross sections, depending on an $N_{\text{conn}}$-dependent saturation scale (“dynamical saturation scale”).

The $N_{\text{conn}}$ dependence of $x^{\pm}$ distributions is guided by the saturation scale, and nothing else.

### D. Generalized AGK theorem

Equations (132) and (122) tell us that also in $A+B$ scatterings, the partonic structure is given by $G_{\text{QCD}}$, and therefore also the $p_t$ distribution of the outgoing partons is encoded in the single Pomeron expression $G_{\text{QCD}}$, for any $N_{\text{conn}}$. Only the saturation scales $Q_{\text{sat}}^2$ depend on $N_{\text{conn}}$, and these saturation scales suppress small-$p_t$ particle production, but should not affect high-$p_t$ results, as sketched in Fig. 41. To be a bit more quantitative: As discussed in detail in Sec. 3 of Ref. [7], in the case of pp scattering with a single Pomeron (of the form $G_{\text{QCD}}$) involved, one can deduce an inclusive dijet cross

<figure><img src="image_43.png" /><figcaption>FIG. 41. Sketch of the suppression of low-$p_t$ partons with increasing $Q_{\text{sat}}^2$ for a single Pomeron, where the red line corresponds to some minimum value $Q_{\text{sat min}}^2$.</figcaption></figure>

section of the form

$$ E_3 E_4 \frac{d^6 \sigma_{\text{incl}}^{\text{single Pom}}}{d^3 p_3 d^3 p_4} [Q_0^2] = \sum_{klmn} \iint dx^+ dx^- f_{\text{PDF}}^k(x^+, Q_0^2, \mu_F^2) f_{\text{PDF}}^l(x^-, Q_0^2, \mu_F^2) \times \frac{1}{32s\pi^2} \sum_{\mathcal{M}^{kl \to mn}} |\mathcal{M}^{kl \to mn}|^2 \delta^4(p_1 + p_2 - p_3 - p_4) \frac{1}{1 + \delta_{mn}}, \quad (133) $$

with explicit expressions for the PDFs given in terms of elements (or modules) of $G_{\text{QCD}}$, like vertices and evolution functions. I indicate explicitly the dependence on the low virtuality cutoff $Q_0^2$. The differential cross section in Eq. (133) and the one in Eq. (132) refer to the same cross section $\sigma_{\text{incl}}^{\text{single Pom}}$, just in Eq. (133) one integrates over $x^{\pm}$, whereas in Eq. (132) the momenta $p_3$ and $p_4$ are integrated out. These momenta refer to the outgoing partons, from the Born process, and $\mu_F^2$ is the factorization scale, being of the same order as the outgoing momenta. Assuming $f_{\text{PDF}}^k(x^+, Q_0^2, \mu_F^2)$ to be independent of $Q_0^2$ for $\mu_F^2 > Q_0^2$, then the only effect of $Q_0^2$ is a reduction of particles with $p_t^2$ below $Q_0^2$. In the case of $A+B$ scattering for given $N_{\text{conn}}$ [see Eq. (132)], one has (up to a factor) the same formula, just $Q_0^2$ has to be replaced by $Q_{\text{sat}}^2(N_{\text{coll}}, x^+, x^-)$, with correspondingly a reduction of particle production below this scale. The minimum bias (MB) inclusive cross section may be written as a superposition of the different contribution for given values of $N_{\text{conn}}$ with weights $w^{(N_{\text{conn}})}$,

$$ E_3 E_4 \frac{d^6 \sigma_{\text{incl}}^{\text{AB (MB)}}}{d^3 p_3 d^3 p_4} \propto \sum_{N_{\text{conn}}=1}^{\infty} w^{(N_{\text{conn}})} E_3 E_4 \times \frac{d^6 \sigma_{\text{incl}}^{\text{single Pom}}}{d^3 p_3 d^3 p_4} [Q_{\text{sat}}^2(N_{\text{conn}}, x^+, x^-)]. \quad (134) $$

If one is only interested in $p_t^2$ values bigger than the saturation scales, then one may replace the satura-
<figure><img src="image_44.png" /><figcaption>FIG. 42. The ratio $R_{AGK}$ (see text) for minimum bias PbPb collisions at 5.02 TeV.</figcaption></figure>

tion scales with some constant value, say $Q_0^2$, and using $\sum_{N_{\text{conn}}=1}^{\infty} w^{(N_{\text{conn}})} = 1$, one finds the “generalized AGK theorem”,

$$E_3 E_4 \frac{d^6 \sigma_{\text{incl}}^{\text{AB (MB)}}}{d^3 p_3 d^3 p_4} = A B \times E_3 E_4 \frac{d^6 \sigma_{\text{incl}}^{\text{single Pom}}}{d^3 p_3 d^3 p_4} [Q_0^2] \quad (\text{at large } p_t), \quad (135)$$

where the “normalization constant” $n$ in Eq. (114) has been used to ensure the factor $AB$. The term “large $p_t$” means $p_t^2$ values bigger than all the “relevant” $Q_{\text{sat}}^2$ values (with negligible weights $w^{(N_{\text{conn}})}$).

In the following, it will be shown that AGK really works in practice. One defines $p_t$ to be the transverse momentum of particle 3 in Eq. (135), and one integrates over the longitudinal momentum of particle 3 and the three momentum components of particle 4. Then one computes the ratio of the full Monte Carlo simulation over $AB$ times the single Pomeron distribution,

$$R_{\text{AGK}} = \left. \frac{d\sigma_{\text{incl}}^{\text{AB (MB)}}}{dp_t} \right|_{\text{Full MC}} / \left\{ AB \frac{d\sigma_{\text{incl}}^{\text{single Pom}}}{dp_t} [Q_0^2] \right\}, \quad (136)$$

showing the result in Fig. 42. The single Pomeron result is calculated numerically (without Monte Carlo). The simulations are done for minimum bias PbPb collisions at 5.02 TeV. One can see that the ratio is close to unity for large values of $p_t$, whereas low $p_t$ values are suppressed. In other words, AGK holds at high $p_t$. Not shown here, but also in $pp$, the full simulation over the single Pomeron reference curve is close to unity at large $p_t$ and one gets therefore $R_{AA} \approx 1$.

## E. Conclusion

Let me summarize this (crucial) section:

(i) It is first recalled that the validity of the AGK theorem (needed to get binary scaling) is badly violated, which

manifests itself by $R_{\text{deform}}(x^+, x^-) \neq 1$, and one understands that this is unavoidable.

(ii) It is also recalled that binary scaling is not always expected, but only at high $p_t$ (concerning the outgoing partons), and that the behavior of $p_t$ distributions depends on the internal structure of the single Pomeron $G$. So far, I was using $G = G_{\text{QCD}}$, with $G_{\text{QCD}}$ (explained in Ref. [7]) essentially based on DGLAP evolutions preceding a hard QCD scattering. This association seems to be wrong.

(iii) It is pointed out that a well-known feature is completely missing—saturation—and this problem seems to be related to the before-mentioned problem of a wrong association $G = G_{\text{QCD}}$.

(iv) It is proposed to solve both problems, by postulating the association $G = k \times G_{\text{QCD}}(Q_{\text{sat}}^2, x^+, x^-)$, with $k$ being inversely proportional to the deformation function, which defines a saturation scale $Q_{\text{sat}}^2$, depending on the “environment” in terms of the Pomeron connection number $N_{\text{conn}}$, which replaces the virtuality cutoff $Q_0^2$ usually used in DGLAP evolutions. In this way, one incorporates saturation.

(v) One can prove that for $A+B$ scattering, inclusive cross sections—with respect to $x^\pm$—are always expressed in terms of the single Pomeron cross section, but depending on an $N_{\text{conn}}$-dependent saturation scale.

(vi) One can prove that for minimum bias $A+B$ scattering, the inclusive cross section—with respect to the transverse momenta of the outgoing partons, for large transverse momenta—is equal to $AB$ times the one of the corresponding single Pomeron cross section. One refers to this as the “generalized AGK theorem”, valid at high $p_t$, in a scenario with energy sharing.

As a final remark: Within a rigorous parallel scattering scenario (which seems mandatory), and respecting energy conservation (which seems mandatory as well), the only way to not get in contradiction with factorization and binary scaling seems to be the consideration of saturation via $G = k \times G_{\text{QCD}}(Q_{\text{sat}}^2)$ with $k$ being inversely proportional to the deformation function. In this sense, parallel scattering, energy conservation, saturation, and factorization (and binary scaling in $A+B$) are deeply connected.

## IX. HIGH- AND LOW-$p_t$ DOMAIN

In Sec. VIII, I have discussed in great detail the “asymptotic behavior”, i.e., the fact that one recovers in the new formalism the validity of the AGK theorem at high transverse momenta, saying that for inclusive cross sections everything breaks down to the single Pomeron case (see Fig. 43), although one has in reality many scatterings happening in parallel.

To compute the inclusive cross section versus $p_t$ (of the outgoing partons), one needs to consider the internal structure of $G_{\text{QCD}}$, which is discussed in detail in Secs. 2 and 3 of Ref. [7], where I show that $G_{\text{QCD}}$ is given as a sum of several contributions: There is the “sea-sea” contribution, where one has a “pseudosoft block” preceding the first perturbative par-
<figure><img src="image_45.png" /><figcaption>FIG. 43. The only diagram contributing to inclusive cross sections in A+B scattering: a single cut Pomeron (multiplied by AB).</figcaption></figure>

ton, as indicated in Fig. 44(a). Here, a sea quark or a gluon is the first parton entering the parton ladder. The vertices $F_{\text{sea}}$ represent couplings to the projectile and target nucleons. Then there is the “val-val” contribution where on both sides a valence quark is the first parton entering the parton ladder, as shown in Fig. 44(b). The vertices $F_{\text{val}}$ represent couplings to the projectile and target nucleons. In both cases, the central part is a parton ladder based on DGLAP parton evolutions [3,29,30]. The inelastic processes corresponding to the two cut diagrams of Fig. 44 are shown in Fig. 45. In addition to sea-sea and val-val, one has also the mixed contributions “sea-val” and “val-sea”. In Ref. [7], one finds explicit expressions for the four cases, expressed in terms of the modules $F_{\text{sea}}$, $F_{\text{val}}$, $E_{\text{psoft}}$, $E_{\text{QCD}}$, and “Born” (the QCD matrix elements of the Born process). As shown in Ref. [7], one may rearrange the integrations such that one may define parton distribution functions $f_{\text{PDF}}$, corresponding to the sum of two diagrams as shown in Fig. 46, which allows one to write jet cross sections as an integral over PDFs and a QCD matrix element [see Eq. (133)], which amounts to factorization.

<figure><img src="image_46.png" /><figcaption>FIG. 44. The contributions (a) $G_{\text{QCD}}^{\text{sea-sea}}$ and (b) $G_{\text{QCD}}^{\text{val-val}}$.</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>FIG. 45. The two inelastic processes corresponding to the two cut diagrams of Fig. 44.</figcaption></figure>

I use these PDFs to compute the jet (parton) cross section for pp at 13 TeV, using Eq. (133), integrating out the momentum of the second parton and the azimuthal angle of the first parton, which finally gives (see Ref. [7])

$$ \frac{d^2\sigma}{dy dp_t^2} = \sum_{klmn} \int dx f_{\text{PDF}}^k(x_1, \mu_F^2) f_{\text{PDF}}^l(x_2, \mu_F^2) \\ \times \frac{\pi \alpha_s^2}{s^2} \left\{ \frac{1}{g^4} \sum_{\vec{l}} |M^{kl \to mn}(s, t)|^2 \right\} \frac{1}{1 + \delta_{mn}} x_1 x_2 \frac{1}{x}, \qquad (137) $$

$$ \text{with } x_1 = x + \frac{p_t}{\sqrt{s_{\text{pp}}} e^y}, \quad x_2 = \frac{x_1}{x} \frac{p_t}{\sqrt{s_{\text{pp}}} e^{-y}}, \\ s = x_1 x_2 \sqrt{s_{\text{pp}}}, \quad t = -p_t x_1 \sqrt{s_{\text{pp}} e^{-y}}, \qquad (138) $$

with $\{\cdots\}$ being the form in which the squared matrix elements are usually tabulated, with $\alpha_s = g^2/4\pi$. I define the parton yield $dn/dp_t dy$ as the cross section $d\sigma/dy dp_t^2$, divided by the inelastic pp cross section, times $2p_t$, showing the result in Fig. 47. In addition to the results based on EPOS PDFs (red solid line), I show the corresponding curves based on “The Coordinated Theoretical-Experimental Project on QCD” (CTEQ) PDFs [31] (green dashed line), the full EPOS simulation (blue circles), and experimental data from ATLAS [32] (black triangles). At large values of $p_t$, all the different distributions 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.

The possibility of using factorization is extremely useful, when one is interested in rare processes such as particle production at large transverse momentum, and it is very important that in EPOS4 one recovers factorization at large $p_t$, since it is a crucial element and its violation would simply disqualify

<figure><img src="image_48.png" /><figcaption>FIG. 46. The two contributions of the parton distribution function (a common remnant vertex factor V is not shown).</figcaption></figure>
<figure><img src="image_49.png" /><figcaption>FIG. 47. Parton yield $dn/dp_t dy$ for $pp$ at 13 TeV. I show results based on EPOS PDFs (red solid line), CTEQ PDFs (green dashed line), the full EPOS simulation (blue circles), and experimental data from ATLAS (black triangles).</figcaption></figure>

the model. However, the majority of all particles are produced at low $p_t$, and even when one is interested in high-$p_t$ particles, one needs to worry about possible interactions with the “bulk” of low-$p_t$ particles. And here, one needs to employ the full multiple scattering machinery, and this concerns $pp$ and nuclear scatterings. Multiple Pomeron configurations are generated using the probability law [see Eqs. (81) and (91)]

$$P(K) = \prod_{k=1}^{AB} \left[ \frac{1}{m_k!} \prod_{\mu=1}^{m_k} G_{k\mu} \right] \times W_{AB}(\{x_i^+\}, \{x_j^-\}), \quad (139)$$

with $G_{k\mu} = G(x_{k\mu}^+, x_{k\mu}^-, s, b_k)$ [see Eq. (52)] being expressed in terms of $G_{QCD}$ via the fundamental equation (114), by introducing dynamical saturation scales, as discussed in Sec. VIII. In Sec. 3 of Ref. [7], it is explained in detail how to generate partons, based on the explicit expressions for $G_{QCD}$, which will allow to obtain “partonic configurations”, as shown in Fig. 48 for the example of a collision of two nuclei A and B, where (for simplicity) each nucleus is 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]. There are two scatterings of sea-sea type, and one of val-sea type. One considers each incident nucleon as a reservoir of three valence quarks plus quark-antiquark pairs. The “objects” which represent the “external legs” of the individual

<figure><img src="image_50.png" /><figcaption>FIG. 48. Partonic configuration of two colliding nuclei, 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.</figcaption></figure>

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 $\bar{3}$ color representation.

The transition from partonic configurations (as in Fig. 48) to strings, based on color flow diagrams, is discussed in detail in Sec. 4 of Ref. [7]. A short summary will be given in the following.

Let me for simplicity consider the quark-antiquark option for the 3 and the $\bar{3}$ color representations, and first of all look at the “sea” cases (on the projectile or the target side) of Fig. 48. In each case, a quark-antiquark pair is emitted as final state timelike (TL) parton pair (marked 1, 2, 3, 4, and 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 one does not treat explicitly. Then emerging from this soft Pomeron, one sees a first perturbative SL gluon (another possibility is the emission of a quark), which initiates the partonic cascade. In the case of “val”, one also has a quark-antiquark pair as an external leg, but here first an antiquark is emitted as a TL final particle (marked 6), plus an SL quark starting the partonic cascade.

In the case of multiple scattering as in Fig. 48, 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-\bar{s}$, then the remnant for an incident proton has flavor uds. In addition, the remnants get massive, much more than simply resonance excitation. One may have remnants with masses of $10 \text{ GeV/c}^2$ or more, which contribute significantly to particle production (at large rapidities).
<figure><img src="image_51.png" /><figcaption>FIG. 49. Configuration colorwise equivalent to the one of Fig. 48. The outgoing antiquarks are drawn as incoming quarks (arrows towards vertices).</figcaption></figure>

In the following, I discuss the color flow for a given configuration, for example, the one in Fig. 48. Since the remnants are by construction white, one does not need to worry about them; one just considers the rest of the diagram. In addition, colorwise, the “soft Pomeron” part behaves as a gluon. Finally, I use the following convention for the SL partons which are immediately emitted: one uses for the quarks an array away from the vertex, and for the antiquarks an array towards the vertex. The diagram equivalent to Fig. 48 is then the one shown in Fig. 49. Based on Fig. 49, considering the fact that, in the parton evolution and the Born process, the gluons are emitted randomly to the right or to the left, I show in Fig. 50 a possible color flow diagram for the three scatterings. Horizontal lines refer to TL partons, which later undergo a timelike cascade, while the vertical lines refer to spacelike intermediate partons. I added integers just to mark the different TL partons. For the leftmost scattering, starting from one end, say “1”, one follows the color flow to “5”, and then starting from 6 to 10, so one gets two chains: 1-2-3-4-5 and 6-7-8-9-10. The end partons of each chain are always quarks or antiquarks; 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.

In the above example, I considered (for simplicity) only gluons, and no timelike cascade. In general, the situation is a bit more complicated, as shown in Fig. 51. Here, one of the gluons emitted in the Born process splits into a gluon (3), a quark (4), and an antiquark (5). Two gluons emitted in the SL cascade, split into the partons 9,10 and 12,13. One

<figure><img src="image_52.png" /><figcaption>FIG. 50. A possible color flow diagram for the three scatterings of Fig. 49.</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>FIG. 51. A color flow diagram with SL and TL cascades.</figcaption></figure>

again follows the color flow, always starting and ending with a single arrow, and one identifies the following chains: 1-2-3-4, then 5-6-7, then 8-9-10-11-12, and finally 13-14.

All 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 as described in Ref. [5], where one also discusses the string decay.

## X. SECONDARY INTERACTIONS: THE ROLE OF CORE, CORONA, AND REMNANTS, AT RHIC AND LHC ENERGIES

The S-matrix part discussed so far concerns “primary scatterings”, happening instantaneously at $t=0$. As a result, one obtains multiple Pomeron configurations, which translates into complex partonic configurations, and eventually into kinky strings and remnants, as discussed in Sec. IX (for details see Ref. [7]). String decay traditionally produces string segments which correspond to hadrons. But one considers the possibility of having a dense environment, and here the string segments cannot “evolve” into hadrons. So one uses the term “prehadrons” for these segments, and they either “fuse” to produce the core, or become hadrons if they escape the core. A similar argument is used for excited remnants, which may decay into hadrons (see Fig. 48), but again this may happen in a dense area, so one names them prehadrons as well. Based on these prehadrons from string or remnant decay, a core-corona procedure will be employed, which allows identifying the core, which will then be treated as a fluid that evolves and eventually decays into hadrons, which still may collide with each other. In the present paper, I focus on PbPb collisions at 5.02 TeV and AuAu scattering at 200 GeV. In Sec. 3 of Ref. [8], one shows explicitly partonic configurations (similar as Fig. 48) but for both LHC and at RHIC energies. They are similar, but with decreasing energy, it becomes simply more and more likely that the Pomerons are replaced by purely soft ones. Also, the Pomerons get less energetic, producing fewer particles. A dedicated paper on lower RHIC energies is in preparation.

Let me discuss some more details about core and corona. Based on the prehadrons from strings and remnants, one employs a so-called core-corona procedure (introduced in Ref. [33], updated in Ref. [34]), at some given (early) proper time $\tau_0$, to separate “core prehadrons” from “corona prehadrons”. The former constitute bulk matter and will be
<figure><img src="image_54.png" /><figcaption>FIG. 52. The prehadron yield as a function of space-time rapidity, for different centralities in PbPb collisions at 5.02 TeV. The curves refer to all prehadrons (red solid curve), all core prehadrons (red dotted curve), prehadrons from remnant decay (blue solid curve), and core prehadrons from remnant decay (blue dotted curve).</figcaption></figure>

treated via hydrodynamics; the latter become simply hadrons and propagate with reduced energy (due to the energy loss). For details, see Sec. 3 of Ref. [8].

In the following, it will be tried to understand the relative importance of the core part, and of the fraction coming from remnant decay. In Fig. 52, I show results for different centralities in PbPb collisions at 5.02 TeV, namely (from top to bottom), 0–10%, 30–40%, 70–80%, and 80–100% (based on the distribution of the impact parameter). I plot four different curves: all prehadrons (red solid curve), all core prehadrons (red dotted curve), prehadrons from remnant decay (blue solid curve), and core prehadrons from remnant decay (blue dotted curve). The remnant contributions show up preferentially at large rapidities, and in all cases, they do contribute to the core. Comparing the red solid and dotted curves, one sees that the core fraction (ratio of $dn/d\eta_s$ of the core contribution over all) is in all cases substantial: 0.97 for 0–10%, 0.95 for 30–40%, 0.85 for 70–80%, and 0.78 for 80–100%. One also sees that this “core dominance” extends over a wide rapidity range.

In Fig. 53, I show results for different centralities in AuAu collisions at 200 GeV. One sees that the core still dominates, but the core fractions are significantly smaller compared to PbPb at 5.02 TeV. Since the overall yields decrease with decreasing energy, the relative importance of the remnant

<figure><img src="image_55.png" /><figcaption>FIG. 53. Same as Fig. 52, but for AuAu at 200 GeV.</figcaption></figure>

contributions (remnant over all) increases. But the remnant contribution at central rapidities remains small.

Having identified core prehadrons, one computes the corresponding energy-momentum tensor $T^{\mu\nu}$ and the flavor flow vector at some space-time position $x$ at initial proper time $\tau = \tau_0$ as

$$T^{\mu\nu}(x) = \sum_i \frac{p_i^\mu p_i^\nu}{p_i^0} g(x - x_i) \quad (140)$$

and

$$N_q^\mu(x) = \sum_i \frac{p_i^\mu}{p_i^0} q_i g(x - x_i), \quad (141)$$

with $q_i \in u, d, s$ being the net flavor content and $p_i$ the four-momentum of prehadron $i$. The function $g$ is some Gaussian smoothing kernel (see Ref. [8] for more details). The Lorentz transformation into the comoving frame provides the energy density $\epsilon$ and the flow velocity components $v^i$, which will be used as the initial condition for a hydrodynamical evolution [34,35]. This is based on the hypothesis that equilibration happens rapidly and affects essentially the space components of the energy-momentum tensor.

In Fig. 54, I show the energy density at the initial proper time $\tau_0$ as a function of the transverse coordinate $r$ for different centralities (defined via impact parameter) in PbPb collisions at 5.02 TeV (upper plot) and in AuAu collisions
<figure><img src="image_56.png" /></figure>

<figure><img src="image_57.png" /><figcaption>FIG. 54. Energy density at the initial proper time $\tau_0$ as a function of the transverse coordinate $r$, for an azimuthal angle $\phi = 0$ (red solid curve) and $\phi = \pi/2$ (red dotted curve). The blue dashed lines represent the freeze-out energy density. I show results for different centralities (defined via impact parameter) in PbPb collisions at 5.02 TeV (upper plot) and in AuAu collisions at 200 GeV (lower plot).</figcaption></figure>

at 200 GeV (lower plot). The blue dashed lines represent the freeze-out energy density. For each event, one determines (based on the energy density distribution) the event plane angle $\psi$ and rotates the system accordingly (to have after rotation event plane angles zero). The plots in Fig. 54 represent averages over such rotated events; the solid lines correspond to azimuthal angles $\phi = 0$, and the dotted lines to $\phi = \pi/2$. The difference between the two lines reflects the azimuthal asymmetry. Even in peripheral PbPb collisions, there is some core production, and one gets actually an energy density of about 4–5 GeV/fm$^3$ for 70–80% centrality, but the radial extension is small, and the lifetime as well. The numerical values in fm/c used for $\tau_0$ are between 1 (peripheral) and 1.5 (central) for PbPb and always 1 for AuAu.

It follows a viscous hydrodynamic expansion. Starting from the initial proper time $\tau_0$, the core part of the system evolves according to the equations of relativistic viscous hydrodynamics [34,35], where one uses presently $\eta/s = 0.08$. The “core-matter” hadronizes on some hypersurface defined by a constant energy density $\epsilon_H$ (presently 0.57 GeV/fm$^3$). In earlier versions [36], one used a so-called Cooper-Frye procedure. This is problematic in particular for small systems: not only do energy and flavor conservation become important, but one also encounters problems due to the fact that one gets small “droplets” with huge baryon chemical potential, with strange results for heavy baryons. In EPOS4, one uses systematically microcanonical hadronization, as discussed in Ref. [8]. After the hadronization of the fluid, the created hadrons as well as the corona prehadrons (having been promoted to hadrons) may still interact via hadronic scatterings, and here one uses UrQMD [37,38].

In the following, I will study core and corona contributions to hadron production. I will distinguish the following:

(A) The core+corona contribution: primary interactions (S-matrix approach for parallel scatterings), plus core-corona separation, hydrodynamic evolution, and microcanonical hadronization, but without hadronic rescattering.

(B) The core contribution: as case A, but considering only core particles.

(C) The corona contribution: as case A, but considering only corona particles.

(D) The full EPOS4 scheme: as case A, but in addition hadronic rescattering.

In cases A, B, and C, one needs to exclude the hadronic afterburner, because the latter affects both core and corona particles, so in the full approach, the core and corona contributions are not visible anymore.

In Fig. 55 (upper plot), I show ratios X/“core+corona” versus $p_t$, with X being the “corona” contribution (blue), the “core” (green), and the “full” contribution (red), for PbPb collisions at 5.02 TeV, for (from top to bottom) pions ($\pi^\pm$), kaons ($K^\pm$), protons ($p$ and $\bar{p}$), and lambdas ($\Lambda$ and $\bar{\Lambda}$). The four columns represent four different centrality classes, namely, 0–5%, 20–40%, 60–80%, and 80–100%. Looking at the green (core) and blue (corona) curves, one observes that the core contribution increases with centrality, but it also increases with the hadron mass (from top to bottom).
<figure><img src="image_58.png" /><figcaption>FIG. 55. The X/"core+corona" ratio as a function of $p_t$ (for $|\eta| < 1$), with X being the "corona" contribution (blue), the "core" (green), and the "full" contribution (red), for four centrality classes and four different particle species, for PbPb at 5.02 TeV (upper plot) and AuAu at 200 GeV (lower plot).</figcaption></figure>

Concerning the $p_t$ dependence, one observes a maximum of the green core curves around 1–2 GeV/c; at very low $p_t$ the core contribution goes down, so even at very small $p_t$ values the corona contributes. The crossing of the green core

and the blue corona curves (core = corona) occurs between around 2 GeV/c (mesons, peripheral) and 5 GeV/c (baryons, central). The red curve, "full" over "core + corona", represents the effect of the hadronic cascade in the "full" case. The pions are not much affected, but for kaons and even more for protons and lambdas, rescattering makes the spectra harder. One should keep in mind that rescattering involves particles from fluid hadronization, but also corona particles from hard processes. Concerning the baryons, rescattering reduces (considerably) low $p_t$ yields, due to baryon-antibaryon annihilation.

In Fig. 55 (lower plot), I show the corresponding results for AuAu collisions at 200 GeV, being similar compared to PbPb at 5.02 TeV, but the "core" contributions are weaker.

XI. RESULTS

In this section, I show simulation results compared to data. I will not add too many comments to each curve; the main purpose is to check if the concepts discussed in the previous sections give a coherent picture (and reproduce the data) or not.

Although the $p_t$ spectra cover usually many orders of magnitude, I have chosen the dimensions such that differences of 10% between data and simulations are always visible.

I will show many results for AuAu collisions at 200 GeV and some results for PbPb at 5.02 TeV. From the theory point of view, the high-energy case contains in principle everything (as discussed in the preceding sections); one does not need to add "features" at lower energies, simply certain phenomena "die out" when reducing the energy: The number of Pomerons per nucleon-nucleon collision gets smaller, the Pomerons get less energetic, and the remnant contributions get relatively more important.

A. Spectra for AuAu at 200 GeV

In this section, I will show particle spectra for AuAu collisions at 200 GeV. From Fig. 53, one knows that at this energy the remnant contribution is important, in particular beyond space-time rapidities $|\eta_s|$ of 2 ($\eta_s$ is numerically similar to the pseudorapidity $\eta$). One also knows from Fig. 55 that at central rapidity the core largely dominates, up to $p_t$ values of 3–5 GeV/c.

In Fig. 56, I show pseudorapidity distributions of charged particles in AuAu collisions at 200 GeV, comparing EPOS4 simulations (lines) to BRAHMS data [39] (points). I show results for different centralities, from top to bottom: 0–5%, 5–10%, ..., 40–50%.

In Figs. 57 and 58, I show transverse mass distributions of $\pi^+, \pi^-, K^+, K^-$, $p$, and $\bar{p}$ for different centrality classes (as indicated in the figures), comparing EPOS4 simulations (lines) to data from STAR [40]. From top to bottom (for each subplot), the curves are multiplied by $3^{-i}$ ($i = 0, 1, 2, 3, ...$).

In Fig. 59, I show transverse momentum distributions of $\pi^+, \pi^-, K^+$, and $K^-$ in central (0–5%) AuAu collisions at 200 GeV for different rapidities (as indicated in the plots). EPOS4 simulations (lines) are compared to data from BRAHMS [41]. From top to bottom, the curves are multiplied by $3^{-i}$ ($i = 0, 1, 2, 3, ...$).
<figure><img src="image_59.png" /><figcaption>FIG. 56. Pseudorapidity distributions of charged particles in AuAu collisions at 200 GeV, comparing EPOS4 (lines) to BRAHMS data (points).</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>FIG. 57. Transverse mass distributions of $\pi^{+}$, $\pi^{-}$ in AuAu collisions at 200 GeV for different centrality classes. EPOS4 simulation (lines) are compared to data from STAR. From top to bottom, one multiplies the curves by $3^{-i}$ ($i = 0, 1, 2, 3, ...$).</figcaption></figure>

<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /><figcaption>FIG. 58. Same as Fig. 57, but for $K^{+}$, $K^{-}$, $p$, and $\bar{p}$.</figcaption></figure>
<figure><img src="image_63.png" /></figure>

<figure><img src="image_64.png" /><figcaption>FIG. 59. Transverse momentum distributions of $\pi^{+}$, $\pi^{-}$, $K^{+}$, and $K^{-}$ in central (0–5%) AuAu collisions at 200 GeV for different rapidities. EPOS4 simulations (lines) are compared to data from BRAHMS. From top to bottom, the curves are multiplied by $3^{-i}$ ($i = 0, 1, 2, 3, ...$).</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>FIG. 60. Transverse momentum distributions of $\pi^{+}$, $\pi^{-}$ in AuAu collisions at 200 GeV at central rapidity for different centralities. EPOS4 simulations (lines) are compared to data from PHENIX [42]. From top to bottom, the curves are multiplied by $3^{-i}$ ($i = 0, 1, 2, 3, 4$).</figcaption></figure>

In Figs. 60 and 61, I plot transverse momentum distributions of $\pi^{+}$, $\pi^{-}$, $K^{+}$, $K^{-}$, $p$, and $\bar{p}$ at central rapidity for different centralities. EPOS4 simulations (lines) are compared to data from PHENIX [42]. From top to bottom, the curves are multiplied by $3^{-i}$ ($i = 0, 1, 2, 3, 4$).

In Fig. 62, I show transverse momentum distributions of $K_{0}$, $\Lambda$, $\bar{\Lambda}$, $\Xi^{-}$, $\bar{\Xi}^{+}$, and $\Omega$ at central rapidity for different centralities. EPOS4 simulations (lines) are compared to data from STAR [43] (for $K_{0}$, $\Lambda$, $\bar{\Lambda}$) and [44]. From top to bottom, the curves are multiplied by $3^{-i}$ ($i = 0, 1, 2, 3, 4$).

In general, the simulation results are close to the data, concerning identified particles as pions, kaons, and protons, as well as hyperons.

Concerning protons and antiprotons in Fig. 61, the data drop remarkably at low $p_{t}$ which is not seen in the simulations, and as well not in other proton and antiproton data (see, for example, Fig. 58).

In all cases, one suppresses (or not) feed-down from decays, to be consistent with the corresponding experimental data.

## B. Flow harmonics for AuAu at 200 GeV

The so-called flow harmonics $v_{n}$ are important observables, characterizing anisotropic azimuthal flow, defined as

$$v_n = \langle \cos(n(\phi - \psi_n)) \rangle, \quad (142)$$
<figure><img src="image_66.png" /><figcaption>FIG. 61. Same as Fig. 60, but for $K^{+}$, $K^{-}$, $p$, and $\bar{p}$.</figcaption></figure>

with $n$ being the order of the flow harmonic, $\phi$ the azimuthal angle, and $\psi_n$ the event plane angle of harmonic $n$. Several methods have been developed over the years to determine $v_n$ from the momentum vectors of the observed particles. In the following, I will compare EPOS4 $v_n$ results with experimental data, employing exactly the same methods as used in the

<figure><img src="image_67.png" /><figcaption>FIG. 62. Transverse momentum distributions of $K_0$, $\Lambda$, $\bar{\Lambda}$, $\Xi^{-}$, $\Xi^{+}$, and $\Omega$ in AuAu collisions at 200 GeV at central rapidity for different centralities. EPOS4 simulations (lines) are compared to data from STAR. From top to bottom, one multiplies the curves by $3^{-i}$ ($i=0, 1, 2, 3, 4$).</figcaption></figure>
<figure><img src="image_68.png" /><figcaption>FIG. 63. Pseudorapidity dependence of $v_2$ of charged particles for different multiplicity classes in AuAu collisions at 200 GeV. I compare the simulations (red lines) with data from PHOBOS (black points).</figcaption></figure>

experiments, with details being found in the corresponding citations.

In Fig. 63, I plot the pseudorapidity dependence of $v_2$ for different multiplicity classes in AuAu collisions at 200 GeV. I compare the simulations (red lines) with data from PHOBOS [45] (black points). In Fig. 64, I plot the centrality dependence (using $N_{part}$) of $v_2$ in AuAu collisions at 200 GeV. I compare the simulations (red lines) with data from PHOBOS [45] (black points). In Fig. 65, I plot the transverse momentum dependence of $v_2$ for pions (left column), kaons (middle column), and protons (right column) in AuAu collisions at 200 GeV, for different centralities (from top to bottom): 0–10%, 10–20%, 20–30%, 30–40%, 40–50%, and 50–60%. I compare the full simulations (red lines) and those without hadronic cascade (yellow lines) with data from PHENIX [46] (black points). In Fig. 66, I show results for $v_3$ (upper plot) and $v_4$ (lower plot), comparing full simulations (red lines) with data from PHENIX [46] (black points).

### C. Spectra for PbPb at 5.02 TeV

I will show some $p_t$ spectra of identified particles, in PbPb collisions at 5.02 TeV. In Figs. 67 and 68, I show transverse momentum distributions at central rapidity of pions, kaons, and protons in PbPb collisions at 5.02 TeV for different centrality classes (from top to bottom: 0–5%, 5–10%, ..., 80–90%). EPOS4 simulations (lines) are compared to data from

<figure><img src="image_69.png" /><figcaption>FIG. 64. Centrality dependence (using $N_{part}$) of $v_2$ of charged particles in AuAu collisions at 200 GeV. I compare the simulations (red lines) with data from PHOBOS (black points).</figcaption></figure>

<figure><img src="image_70.png" /><figcaption>FIG. 65. Transverse momentum dependence of $v_2$ for pions (left column), kaons (middle column), and protons (right column) in AuAu collisions at 200 GeV, for different centralities (from top to bottom): 0–10%, 10–20%, 20–30%, 30–40%, 40–50%, and 50–60%. I compare the full simulations (red lines) and those without hadronic cascade (yellow lines) with data from PHENIX (black points).</figcaption></figure>
<figure><img src="image_71.png" /></figure>

<figure><img src="image_72.png" /><figcaption>FIG. 66. Transverse momentum dependence of $v_3$ (upper plot) and $v_4$ (lower plot) for pions (left column), kaons (middle column), and protons (right column) in AuAu collisions at 200 GeV, for different centralities. I compare the full simulations (red lines) with data from PHENIX (black points).</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>FIG. 67. Transverse momentum distributions of pions in PbPb collisions at 5.02 TeV for different centrality classes. EPOS4 simulation data (lines) are compared to data from ALICE. From top to bottom, one multiplies the curves by $2^{-i}$ ($i = 0, 1, 2, 3, 4, ...$).</figcaption></figure>

ALICE [47]. From Fig. 55, one knows that corona particles dominate beyond 2–3 GeV/c for peripheral collisions and beyond 4–5 GeV/c for central ones. They are affected by the energy loss parameter in the core-corona procedure, and allow one to determine it. Concerning the kaons and protons, one clearly sees the “flow effect”, i.e., an increase at intermediate $p_t$, compared to pions.

The picture looks consistent, although there is a systematic excess of protons at intermediate $p_t$. It should be emphasized that EPOS4.0.0 has been “tuned” to thousands of experimental results (with only a small fraction shown in this paper), and the aim is to see to what extent one gets a consistent overall picture, and not to optimize a particular curve.

### D. Flow harmonics for PbPb at 2.76 and 5.02 TeV

In the following, I show results for the flow harmonics $v_2$ and $v_3$ in PbPb collisions at LHC energies. In Fig. 69, I plot the pseudorapidity dependence of $v_2\{2\}$ and $v_2\{4\}$ (see Ref. [48] for the definitions) for different multiplicity classes (using percentiles) in PbPb collisions at 2.76 GeV. I compare the EPOS4 simulations (red lines) with data from ALICE [48] (black points). In Fig. 70, I plot the $p_t$ dependence of $v_2$ for different multiplicity classes (using percentiles) in PbPb collisions at 5.02 TeV, for (from left to right) pions,
<figure><img src="image_74.png" /><figcaption>FIG. 68. Same as Fig. 67, but for kaons (upper plot) and protons (lower plot).</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>FIG. 69. Pseudorapidity dependence of $v_2\{2\}$ and $v_2\{4\}$ for different multiplicity classes (using percentiles) in PbPb collisions at 2.76 GeV. I compare the EPOS4 simulations (red lines) with data from ALICE (black points).</figcaption></figure>

protons, kaons ($K_s$), $\phi$ mesons, and lambdas. I compare the simulations (red lines) with data from ALICE [49] (black points). In Fig. 71, I show the corresponding plots for $v_3$. The simulations give in general a decent description of the data, but the simulation results for $v_2$ at high $p_t$ are significantly below the data, pointing to some possible problem related to parton energy loss, which is presently realized as energy loss of prehadrons during the core-corona procedure. Some real parton energy loss procedure in the EPOS4 framework is presently developed.

## XII. SUMMARY AND CONCLUSION

It was first recalled that it is of fundamental importance to realize that multiple nucleon-nucleon scatterings in high-energy nucleus-nucleus ($A+B$) collisions must happen in parallel, and not sequentially. There is nothing like a first, second, third, etc., collision; they are all equal and simultaneous. In this context, I presented a discussion about timescales, and the corresponding applicability of this “parallel scattering scenario” as a function of the collision energy. It was estimated that beyond a collision energy of 24 GeV the parallel scattering scenario is mandatory. The term “parallel scattering” actually refers to both nucleon-nucleon scatterings and parton-parton scatterings in each nucleon-nucleon collision. I then reviewed early work on parallel scattering based on S-matrix theory, using the concept of subprocesses referred
energy limit”, so energy sharing among the Pomerons is not an issue. However, in a realistic scenario—in particular when it is used as a basis for event-by-event Monte Carlo procedures—it is mandatory to include energy-momentum sharing.

I showed how to implement energy-momentum sharing in the S-matrix approach, and how to prove the validity of the AGK theorem, which means that the inclusive cross section (with respect to the Pomeron energy) for A+B scattering is AB times the inclusive cross section of a single Pomeron. As a direct consequence, one gets binary scaling, which means that the inclusive cross section for A+B scattering is AB times the inclusive cross section for pp scattering. Starting from the expression for the inelastic A+B scattering cross section, I showed how to derive a probability law $\sum_{\{m_k\}} \int dX_{AB}P(K) = 1$, with a “configuration” $K = \{\{m_k\}, \{x_{k\mu}^{\pm}\}\}$ representing $m_k$ cut Pomerons per pair $k$, with light-cone momentum fractions $x_{k\mu}^{\pm}$, where $P$ is given as a product of single Pomeron expressions $G$ times some known function $W_{AB}$ (being the result of integrating out all elastic scatterings, i.e., uncut Pomerons). The symbol $G$ refers to a cut diagram representing a single scattering. The calculation $W_{AB}$ (multidimensional integral) could be done, showing that it is not strictly non-negative, which ruins the probability interpretation of $P$.

I showed that it is possible to regularize $W_{AB}$ to make it strictly non-negative, which allows recovering the probability interpretation of $P$, and which allows generating multiple scattering (or multiple Pomeron) configurations. In order to generate partons, one needs to specify the internal structure of the Pomeron $G$. I use the simplest choice, namely, $G = G_{QCD}$, where $G_{QCD}$ refers to a pQCD diagram, with the main element being a parton-parton scattering graph with parton evolutions on both sides and QCD matrix element in the middle. Computing the inclusive cross section (with respect to the variable $p_t$) in A+B scattering, one finds that the AGK theorem is badly violated (which as a consequence ruins binary scaling). So the model is not usable as it is.

In order to understand the problem, I defined a quantity $R_{deform}(x^+, x^-)$, referred to as the “deformation function”, being the ratio of the normalized inclusive A+B cross section over the single Pomeron one, considering the Pomeron light-cone momentum fractions $x^{\pm}$ as variables. One needs $R_{deform}(x^+, x^-) = 1$ in order to validate the AGK theorem, at least in the case $G = G_{QCD}$. But $R_{deform}$ is far from unity, which is the origin of the violation of AGK. I investigated $R_{deform}$ more in detail, and it was understood that a quantity called “connection number”, $N_{conn}$, which counts the number of Pomerons connected to a given pair $i, j$ of projectile and target nucleons, plays an important role. The bigger $N_{conn}$, the bigger the number of Pomerons which need to share the given initial energy of the pair, and therefore, with increasing $N_{conn}$, the deformation function $R_{deform}$ decreases more and more below unity, at large values of $x^{\pm}$. This is unavoidable, and therefore a violation of AGK is unavoidable. At least, one could “quantify the problem”, being able to find a simple parametrization of $R_{deform}$, for given values of $N_{conn}$, for all systems and all energies, with tabulated parameters. So in the following, $R_{deform}$ could be considered to be known.

At this point, I realized that there were two problems: (a) the assumption $G = G_{QCD}$, making the link between the “Pomeron approach” and QCD, seemed to be wrong; (b) saturation was not implemented, which has been known since a long time to play an important role. But fortunately, the two problems are connected, and there is an amazingly simple solution based on saturation scales that solves both problems. Instead of the “naive” assumption $G = G_{QCD}$, I postulated

$$G(x^+, x^-) = \frac{n}{R_{\text{deform}}(x^+, x^-)} G_{\text{QCD}}(Q_{\text{sat}}^2, x^+, x^-), \quad (143)$$

such that $G$ itself does not depend on the environment, in terms of $N_{conn}$, which means that the $N_{conn}$ dependence of $R_{deform}$ is “absorbed” by the saturation scale $Q_{sat}^2$, which in this way depends on $N_{conn}$ and also on $x^{\pm}$. Based on Eq. (143), I could prove the validity of the AGK theorem, saying that the inclusive cross section (with respect to $p_t$) for A+B scattering is AB times the inclusive cross section of a single Pomeron— but only at large $p_t$ (bigger than the relevant saturation scales). In this way, I was finally able to connect the multiple Pomeron approach (for parallel scatterings) and pQCD, by introducing saturation scales. The latter seem indispensable for getting a consistent picture.

I discussed consequences of having, on one hand, a multiple scattering approach for instantaneous parallel scatterings, and on the other hand, the validity of the AGK theorem. The latter allows one to define parton distribution functions and do simple calculations of inclusive cross sections at very high $p_t$ (the same way as models based on factorization do). The former does allow one to do much more, treating problems far beyond the factorization scheme, by using the full multiple scattering machinery.

An important application is the study of collective effects in A+B or pp scattering. Form the multiple (parallel) scattering approach (referred to as primary interactions) one gets first of all a (more or less) large number of prehadrons, which are then subject to secondary interactions, in the form of core-corona separation, hydrodynamic evolution of the core, and microcanonical decay of the plasma. All this has been discussed briefly, not being the main subject of this paper. Finally, I showed some results on heavy-ion collisions at LHC and RHIC energies, covering $p_t$ spectra and flow harmonics of identified hadrons, in order to show what extent the model works.

## ACKNOWLEDGMENTS

K.W. thanks Tanguy Pierog for many contributions during the past decade and in particular the important proposal (and first implementation) in 2015 of using a Pomeron definition as $G \propto G_{QCD}(Q_{sat}^2)$ with a parametrized $G$. He also made a first attempt to take into account the “deformation” by using rescaled longitudinal momentum arguments.

## APPENDIX A: ASYMPTOTIC BEHAVIOR OF T MATRICES FOR pp SCATTERING

Consider a reaction (all particles having mass $m$)

$$1 + 2 \rightarrow 3 + 4 \qquad (A1)$$
in the x-z plane. In the c.m. system, one has ($P = |\vec{p}|$)

$$P_1 = P_2 = P_3 = P_4, \qquad (A2)$$

and

$$E_1 = E_2 = E_3 = E_4. \qquad (A3)$$

The vectors ($E, p_x, p_z$) are given as (for particles 1,2,3,4)

$$(E, p_x, p_z) = \begin{cases} (\sqrt{P^2+m^2}, 0, P) \\ (\sqrt{P^2+m^2}, 0, -P) \\ (\sqrt{P^2+m^2}, P \sin\theta, P \cos\theta) \\ (\sqrt{P^2+m^2}, -P \sin\theta, -P \cos\theta), \end{cases} \qquad (A4)$$

which gives

$$s = 4(P^2 + m^2), \qquad (A5)$$

$$t = 0 - P^2 \sin^2 \theta - P^2 (1 - \cos \theta)^2 = -2P^2(1 - \cos \theta), \qquad (A6)$$

$$u = 0 - P^2 \sin^2 \theta - P^2 (1 + \cos \theta)^2 = -2P^2(1 + \cos \theta). \qquad (A7)$$

Defining $z = \cos \theta$, one gets

$$t = 2\left(\frac{s}{4} - m^2\right)(z - 1), \qquad (A8)$$

$$u = 2\left(\frac{s}{4} - m^2\right)(-z - 1), \qquad (A9)$$

which gives for zero masses the relation

$$z = 1 + \frac{2t}{s}. \qquad (A10)$$

One may always expand the T-matrix (partial wave expansion)

$$\mathbf{T}(s, t) = \sum_{j=0}^{\infty} (2j + 1) \mathcal{T}(j, s) P_j(z) \qquad (A11)$$

with

$$\mathcal{T}(j, s) = \frac{1}{2} \int_{-1}^{1} dz \, \mathbf{T}(s, t) P_j(z) \qquad (A12)$$

due to the orthogonality property of Legendre polynomials,

$$\frac{2n+1}{2} \int_{-1}^{1} dx P_n(x) P_m(x) = \delta_{mn}. \qquad (A13)$$

One analytically continues to the unphysical region of the $1+2 \to 3+4$ process, with

$$\mathbf{T}_{1+2 \to 3+4}(s < 0, t > 0), \qquad (A14)$$

corresponding to the physical region of the $1+\bar{3} \to \bar{2}+4$ ($t$-channel) process, with

$$\mathbf{T}_{1+\bar{3} \to \bar{2}+4}(t > 0, s < 0), \qquad (A15)$$

and after exchanging $s$ and $t$, one gets

$$\mathbf{T}_{1+\bar{3} \to \bar{2}+4}(s > 0, t < 0). \qquad (A16)$$

<figure><img src="image_76.png" /><figcaption>FIG. 72. Contour integrations in $j$ plane.</figcaption></figure>

With this new definition of $\mathbf{T}$, the partial wave expansion reads

$$\mathbf{T}(s, t) = \sum_{j=0}^{\infty} (2j + 1) \mathcal{T}(j, t) P_j(z), \qquad (A17)$$

with

$$\mathcal{T}(j, t) = \frac{1}{2} \int_{-1}^{1} dz \, \mathbf{T}(s, t) P_j(z), \qquad (A18)$$

and with

$$z = 1 + \frac{2s}{t}. \qquad (A19)$$

The Watson-Sommerfeld transform [50,51] amounts to writing the partial wave expansion as

$$\mathbf{T}(s, t) = \frac{1}{2i} \int_C d\bar{j} \frac{1}{\sin \pi \bar{j}} (2j + 1) \mathcal{T}(j, t) P_j(z), \qquad (A20)$$

with a contour integration in the complex $j$ plane, as shown in Fig. 72(a). Opening the contour [see Fig. 72(b)] to integrate along the imaginary axis (with $\text{Re}j = -\frac{1}{2}$), one picks up poles of $\mathcal{T}$ at $j = \alpha_n(t)$, with residues $\beta'_n(t)$, so one gets

$$\mathbf{T}(s, t) = \frac{1}{2i} \int_{C'} d\bar{j} \frac{2j + 1}{\sin \pi \bar{j}} \mathcal{T}(j, t) P_j(z) + \sum \underbrace{\frac{\pi}{\sin \pi \alpha_n(t)} \beta'_n(t) P_{\alpha_n(t)}(z)}_{\beta_n(t)} \qquad (A21)$$

As one will see, the partial wave amplitudes have contributions which alternate in sign, so there is a factor

$$(-1)^j = \exp(i\pi j) \qquad (A22)$$

which diverges on the imaginary axis. The problem will be solved by separating even and odd terms. To see this problem, let me first investigate $\mathcal{T}(j, t)$. One considers $\mathbf{T}(s, t)$ as a function of $z$ and $t$, and writes (Cauchy)

$$\mathbf{T}(z, t) = \frac{1}{2\pi i} \int_C \frac{\mathbf{T}(z', t)}{z' - z} dz', \qquad (A23)$$

where one chooses $C$ as in Fig. 73. One finds (assuming the semicircles do not contribute)

$$\mathbf{T}(z, t) = \frac{1}{2\pi i} \int_{z_0^+}^{\infty} \frac{\mathrm{disc} \, T(z', t)}{z' - z} dz' + \frac{1}{2\pi i} \int_{-\infty}^{-z_0} \frac{\mathrm{disc} \, T(z', t)}{z' - z} dz'. \qquad (A24)$$
<figure><img src="image_77.png" /><figcaption>FIG. 73. Contour integration in z plane.</figcaption></figure>

This so-called dispersion relation may be inserted into Eq. (A18):

$$ \mathcal{T}(j, t) = \frac{1}{2\pi i} \int_{z_0}^{\infty} \mathrm{disc} \, \mathbf{T}(z', t) \frac{1}{2} \int_{-1}^{1} \frac{P_j(z)}{z' - z} dz dz' \\ + \frac{1}{2\pi i} \int_{-\infty}^{-z_0} \mathrm{disc} \, \mathbf{T}(z', t) \frac{1}{2} \int_{-1}^{1} \frac{P_j(z)}{z' - z} dz dz'. \qquad (A25) $$

One may use Neumann's formula

$$ \frac{1}{2} \int_{-1}^{1} \frac{P_j(z)}{z' - z} dz = Q_j(z'), \qquad (A26) $$

with Legendre functions $Q_j$ of the second kind, and get

$$ \mathcal{T}(j, t) = \frac{1}{2\pi i} \int_{z_0}^{\infty} \mathrm{disc} \, \mathbf{T}(z', t) Q_j(z') dz' \\ + \frac{1}{2\pi i} \int_{z_0}^{\infty} \mathrm{disc} \, \mathbf{T}(-z', t) Q_j(-z') dz'. \qquad (A27) $$

Using the symmetry property $Q_j(-z) = (-1)^{j+1}Q_j(z)$ and the definition $D_j^{\pm} = \frac{1}{2\pi i} \int_{z_0}^{\infty} \mathrm{disc} \, \mathbf{T}(\pm z', t) Q_j(z') dz'$, one gets

$$ \mathcal{T}(j, t) = \{D_j^+ + (-1)^{j+1} D_j^-\}, \qquad (A28) $$

which shows the above-mentioned $(-1)^j$ problem. Writing

$$ \mathcal{T}(j, t) = \sum_{\eta \in \{-1, 1\}} \frac{1 + \eta(-1)^j}{2} \{D_j^+ + (-1)^{j+1} D_j^-\} \qquad (A29) $$

allows one to separate even and odd terms,

$$ \mathcal{T}^1(j, t) = D_j^+ - D_j^-, \qquad \mathcal{T}^{-1}(j, t) = D_j^+ + D_j^-, \qquad (A30) $$

as

$$ \mathcal{T}(j, t) = \sum_{\eta \in \{-1, 1\}} \frac{1 + \eta(-1)^j}{2} \mathcal{T}^{\eta}(j, t) \qquad (A31) $$

with a so-called signature $\eta$. The Watson-Sommerfeld transform then reads

$$ \mathbf{T}(s, t) = \sum_{\eta} \frac{1}{2i} \int_{C'} dj \frac{1+\eta(-1)^j}{2} \frac{1}{\sin \pi j} (2j+1) \mathcal{T}^{\eta}(j, t) P_j(z) \\ + \sum_{\eta} \sum_{n_{\eta}} \frac{1+\eta(-1)^{\alpha_{n_{\eta}}(t)}}{2} \beta_{n_{\eta}}(t) P_{\alpha_{n_{\eta}}(t)}(z), \qquad (A32) $$

and there is no $(-1)^j$ problem any more. The poles $\alpha_{n}(t)$ are called (i) even-signature Regge poles ($\eta = +1$) or (ii) odd-signature Regge poles ($\eta = -1$). Asymptotically, for $s \gg |t|$, i.e., $|z - 1| = |\frac{2\eta}{t}| \gg 1$, the Legendre polynomial

$$ P_j(z) = \sum_{k=0}^{j} \binom{j}{k} \binom{j+k}{k} \left(\frac{z-1}{2}\right)^k \qquad (A33) $$

is dominated by its leading term,

$$ P_j(z) \rightarrow \binom{2j}{j} \left(\frac{z-1}{2}\right)^j, \qquad (A34) $$

so one has in this limit

$$ P_j(z) = P_j\left(1 + \frac{2s}{t}\right) \rightarrow \frac{\Gamma(2j + 1)}{\Gamma^2(j + 1)} \left(\frac{s}{t}\right)^j. \qquad (A35) $$

The contour integration (along $j = -\frac{1}{2} + iy$) for $s \gg t$ vanishes, because of factor $(s/t)^{-1/2}$, and one gets (absorbing the $\Gamma$ functions and $t^{-j}$ into $\beta$), for $s \gg t$:

$$ \mathbf{T}(s, t) = \sum_{\eta} \sum_{n_{\eta}} \frac{1 + \eta(-1)^{\alpha_{n_{\eta}}(t)}}{2} \beta_{n_{\eta}}(t) s^{\alpha_{n_{\eta}}(t)}. \qquad (A36) $$

With $\alpha(t)$ being the rightmost (leading) pole, one gets

$$ \mathbf{T}(s, t) = \frac{1 + \eta e^{-i\pi\alpha(t)}}{2} \beta(t) s^{\alpha(t)} \qquad (A37) $$

viewed at the time to correspond to the exchange of some hypothetical “particle” called Reggeon or Pomeron. Absorbing $(1 + \eta e^{-i\pi\alpha(t)})/2$ into $\beta$, and (because of $t \ll s$)

$$ \alpha(t) = \alpha(0) + \alpha't, $$

and neglecting the $t$ dependence of $\beta$, one gets finally

$$ \mathbf{T}(s, t) = \beta s^{\alpha(0)+\alpha'}, \qquad (A38) $$

with $\alpha(0)$ called intercept, and $\alpha'$ called slope. This is referred to as Regge pole expression for the T-matrix [52].

## APPENDIX B: EXPLICIT FORMULAS FOR $W^{AB}$

The $W_{AB}$ expression is given as [see Eqs. (79) and (80)]

$$ W_{AB} = \sum_{l_1=0}^{\infty} \cdots \sum_{l_{AB}=0}^{\infty} \int d\tilde{X}_{AB} \prod_{k=1}^{AB} \left[ \frac{1}{l_k!} \prod_{\lambda=1}^{l_k} -G'_{k\lambda} \right] \\ \times \prod_{i=1}^{A} V \left( x_i^+ - \sum_{k=1}^{AB} \sum_{\lambda=1}^{l_k} \tilde{x}_{k\lambda}^+ \right) \prod_{j=1}^{B} V \left( x_j^- - \sum_{k=1}^{AB} \sum_{\lambda=1}^{l_k} \tilde{x}_{k\lambda}^- \right), \qquad (B1) $$

with

$$ x_i^+ = 1 - \sum_{k=1}^{AB} \sum_{\mu=1}^{m_k} x_{k\mu}^+, \qquad x_j^- = 1 - \sum_{k=1}^{AB} \sum_{\tau(k)=j}^{m_k} x_{k\tau\mu}^- \qquad (B2) $$

Using $G = G_{QCD}$ [see Eq. (48) and the discussion before], one can show [5] that one may obtain an almost perfect fit of the
numerically computed functions $G_{QCD}$, of the form

$$G_{QCDpar}(x^+, x^-, s, b) = \sum_{N=1}^{N_{par}} \alpha_N (x^+ x^-)^{\beta_N}, \quad (B3)$$

$$\alpha_N = (\alpha_{D_N} + \alpha_{D_N}^*) s^{(\beta_{D_N} + \gamma_{D_N} b^2)} e^{-\frac{b^2}{8D_N}}, \quad (B4)$$

$$\beta_N = \beta_{D_N} + \beta_{D_N}^* + \gamma_{D_N} b^2 - \alpha_{\text{part}}, \quad (B5)$$

with $\alpha_{D_N}^* \neq 0$ and $\beta_{D_N}^* \neq 0$ only if $\alpha_{D_N} = 0$, and with $N_{par}$ being a small number (4). This parametric form has been inspired by the asymptotic expressions for T-matrices (see Appendix A). So I will use $G = G_{QCDpar}$ in the following. Furthermore, the vertices are parametrized as

$$V(x) = x^{\alpha_{\text{remn}}}. \quad (B6)$$

1. Doing the integrals

In the following discussion, I use explicitly the (current) value $N_{par} = 4$, for simplicity, but the results can be easily adapted to other values. Using $G = G_{QCDpar}$ and $V(x) = x^{\alpha_{\text{remn}}}$, and defining $G_{Nk\lambda} = \alpha_N(\tilde{x}_{k\lambda}^+ \tilde{x}_{k\lambda}^-)^{\beta_N}$, one gets

$$W_{AB} = \sum_{l_1=0}^{\infty} \cdots \sum_{l_{AB}=0}^{\infty} \int d\tilde{X}_{AB} \prod_{k=1}^{AB} \left[ \frac{1}{l_k!} \prod_{\lambda=1}^{l_k} \sum_{N=1}^{4} -G_{Nk\lambda} \right] \\ \times \prod_{i=1}^{A} \left( x_i^+ - \sum_{k=1 \atop \pi(k)=i}^{AB} \sum_{\lambda=1}^{l_k} \tilde{x}_{k\lambda}^+ \right)^{\alpha_{\text{remn}}} \\ \times \prod_{j=1}^{B} \left( x_j^- - \sum_{k=1 \atop \tau(k)=j}^{AB} \sum_{\lambda=1}^{l_k} \tilde{x}_{k\lambda}^- \right)^{\alpha_{\text{remn}}}. \quad (B7)$$

One has for any $k$, using simply $l = l_k$,

$$\sum_{l=0}^{\infty} \cdots \frac{1}{l!} \prod_{\lambda=1}^{l} \sum_{N=1}^{4} -G_{Nk\lambda} = \sum_{l=0}^{\infty} \cdots \frac{1}{l!} \prod_{N_1=1}^{4} \cdots \prod_{N_l=1}^{4} \prod_{\lambda=1}^{l} -G_{N_l k \lambda}. \quad (B8)$$

Keeping in mind that this is an integrand, and one may exchange integration variables $\tilde{x}_{k\lambda}^{\pm}$ and $\tilde{x}_{k\lambda}'$, one may exchange for a given value of $N$ the terms $G_{Nk\lambda}$ and $G_{Nk\lambda}'$ without changing the result, which means that in the sum over $\prod_{\lambda=1}^l -G_{N_k k\lambda}$ there are $(r_1! r_2! r_3! r_4!)$ identical expressions with $r_N$ times $N_\lambda = N$, for $N = 1, \dots, 4$, which one may write as

$$\prod_{\lambda=1}^{r_1} -G_{1k\lambda} \prod_{\lambda=r_1+1}^{r_1+r_2} -G_{2k\lambda} \prod_{\lambda=r_1+r_2+1}^{r_1+r_2+r_3} -G_{3k\lambda} \prod_{\lambda=r_1+r_2+r_3+r_4}^{-} -G_{4k\lambda}, \quad (B9)$$

with the constraint $r_1 + r_2 + r_3 + r_4 = l$, which will be released when one does the sum $\sum_{l=0}^{\infty}$. Putting it all together,

one gets (using $l_k = l$ and $r_{Nk}$ instead of $r_N$)

$$\sum_{l_k=0}^{\infty} \cdots \frac{1}{l_k!} \prod_{\lambda=1}^{l_k} \sum_{N=1}^{4} -G_{Nk\lambda} = \sum_{r_{1k}=0}^{\infty} \cdots \sum_{r_{4k}=0}^{\infty} \cdots \frac{1}{r_{1k}! \cdots r_{4k}!} \prod_{\lambda=1}^{r_{1k}} -G_{1k\lambda} \cdots \prod_{\lambda=r_{1k}+r_{2k}+r_{3k}+r_{4k}}^{r_{1k}+r_{2k}+r_{3k}+r_{4k}} -G_{Nk\lambda}. \quad (B10)$$

Inserted into Eq. (B7), using $G_{Nk\lambda} = \alpha_N(\tilde{x}_{k\lambda}^+ \tilde{x}_{k\lambda}^-)^{\beta_N}$, one gets

$$W_{AB} = \sum_{r_{11}=0}^{\infty} \cdots \sum_{r_{41}=0}^{\infty} \cdots \sum_{r_{1AB}=0}^{\infty} \cdots \sum_{r_{4AB}=0}^{\infty} \\ \times \int d\tilde{X}_{AB} \prod_{k=1}^{AB} \left[ \frac{1}{r_{1k}! \cdots r_{4k}!} \prod_{\lambda=1}^{r_{1k}} - \alpha_1 (\tilde{x}_{k\lambda}^+ \tilde{x}_{k\lambda}^-)^{\beta_1} \right. \\ \left. \cdots \prod_{\lambda=r_{1k}+r_{2k}+r_{3k}+1}^{r_{1k}+r_{2k}+r_{3k}+r_{4k}} - \alpha_4 (\tilde{x}_{k\lambda}^+ \tilde{x}_{k\lambda}^-)^{\beta_4} \right] \\ \times \prod_{i=1}^{A} \left( x_i^+ - \sum_{k=1 \atop \pi(k)=i}^{AB} \sum_{\lambda=1}^{l_k} \tilde{x}_{k\lambda}^+ \right)^{\alpha_{\text{remn}}} \\ \times \prod_{j=1}^{B} \left( x_j^- - \sum_{k=1 \atop \tau(k)=j}^{AB} \sum_{\lambda=1}^{l_k} \tilde{x}_{k\lambda}^- \right)^{\alpha_{\text{remn}}}. \quad (B11)$$

This will allow one to separate the $x^+$ and $x^-$ integrations. I will use

$$\int d\tilde{X}_{AB} = \int \prod_{k=1}^{AB} \left( \prod_{\lambda=1}^{l_k} d\tilde{x}_{k\lambda}^{+} d\tilde{x}_{k\lambda}^{-} \right) = \frac{\int \prod_{k=1}^{AB} \left( \prod_{\lambda=1}^{l_k} d\tilde{x}_{k\lambda}^{+} \right) \int \prod_{k=1}^{AB} \left( \prod_{\lambda=1}^{l_k} d\tilde{x}_{k\lambda}^{-} \right)}{d\tilde{X}_{AB}^{+} d\tilde{X}_{AB}^{-}}. \quad (B12)$$

where the upper limit $l_k$ is defined as

$$l_k = r_{1k} + r_{2k} + r_{3k} + r_{4k}. \quad (B13)$$

In addition, I define four intervals of integers, as

$$I_{1k} = [1..r_{1k}], \quad (B14)$$

$$I_{2k} = [r_{1k} + 1..r_{1k} + r_{2k}], \quad (B15)$$

$$I_{3k} = [r_{1k} + r_{2k} + 1..r_{1k} + r_{2k} + r_{3k}], \quad (B16)$$

$$I_{4k} = [r_{1k} + r_{2k} + r_{3k} + 1..r_{1k} + r_{2k} + r_{3k} + r_{4k}], \quad (B17)$$

such that

$$[1..l_k] = I_{1k} \cup I_{2k} \cup I_{3k} \cup I_{4k}. \quad (B18)$$