Thermal radiation and inclusive production in the KLN model for ion-ion collisions

We show that in order to obtain a successful description of the transverse momenta distribution for charged particles in ion-ion collisions, one must include a thermal emission term. The temperature of this emission $T_{\rm th}$ turns out to be proportional to the saturation scale, $T_{\rm th} = 1.8/2\pi \,Q_s$. The formalism for the calculation of the transverse momenta spectra in CGC/saturation approach is developed, in which two stages of the process are seen: creation of the colour glass condensate, and hadronization of the gluon jets. Our calculations are based on the observation that even for small values of $p_T$, the main contribution in the integration over the dipole sizes stems from the kinematic region in vicinity of the saturation momentum, where theoretically, we know the scattering amplitude. Non-perturbativ corrections need to be included in the model of hadronization. This model incorporates the decay of a gluon jet with effective mass $m^2_{\rm eff} = 2 Q_s \mu_{\rm soft}$ where $\mu_{\rm soft}$ denotes the soft scale, with the fragmentation functions at all values of the transverse momenta. We use the KLN model which, provides a simple way to estimate the cross sections for the different centrality classes. Comparing the results of this paper with the transverse distribution in the proton-proton scattering, we see two major differences. First, a larger contribution of the thermal radiation term is needed, in accord with higher parton densities of the produced colour glass condensate. Second, even changing the model for the hadronization, without a thermal radiation term, we fail to describe the $p_T$ spectrum. Consequently, we conjecture that theexistence of the thermal radiation term is independent of the model of confinement.


I. INTRODUCTION
In this paper we continue to discuss the processes of multi-particle generation at high energy, in the framework of the Color Glass Condensate(CGC)/saturation approach (see Ref. [1] for the review). The main ideas of our approach have been discussed in our paper (see Ref. [2]) for hadron-hadron scattering at high energies. Here, we consider heavy ion collisions in which, we believe, the distinctive features of the CGC/saturation approach manifest themselves in the clearest way. For theoretical descriptions of these processes we have to develop a non-perturbative approach, since these processes occur at long distances. In particular, we have to deal with the unsolved problem of the confinement of quarks and gluons. Fortunately, in the framework of the CGC/saturation approach, the most basic features of the processes of the multi-particle generation stem from the production of the new phase of QCD: the dense system of partons (gluons and quarks) with a new characteristic scale: saturation momentum Q s (W ), which increases as a function of energy W [3][4][5]. In ion-ion collisions the new phase of QCD is produced with a larger density than in hadron-hadron scattering, and we believe, that the essential properties of this phase will manifest themselves in a clear way. However, the transition from this system of partons to the measured state of hadrons, is still an unsolved problem.
Due to our lack of theoretical understanding of the confinement of quarks and gluons, at the moment, we need to use a pure phenomenological input for the long distance non-perturbative physics. In particular, we wish to use phenomenological fragmentation functions. Hence, our model for confinement is that the parton (quark or gluon) with the transverse momenta of the order of Q s decays into hadrons according to the given fragmentation functions. Experimental data supports this model of hadronization, which provides the foundation of all Monte Carlo simulation programs, and leads to descriptions of the transverse momenta distribution of the hadrons, at the LHC energies.
As an example, we refer to Ref. [6], which shows that next-to-leading order QCD calculations with formation of the hadrons in accord with the fragmentation functions ( [7]), is able to describe the transverse momentum spectra, for the LHC range of energies. In a sense, at present, this model is the best that we can propose to describe multihadron production * .
The space-time picture of the high energy interaction in the CGC/saturation approach, is as follows: The parton configuration in QCD is formed long before the interaction at distances R A /x, where R A -denotes the nucleus radius, and x the fraction of longitudinal momentum carried by the parton which interacts with the target. However, before the collision, the wave function of this partonic fluctuation is the eigenfunction of the Hamiltonian and, therefore, the system has zero entropy. The interaction with the target of size R A destroys the coherence of the parton wave function of the projectile. The typical time, which is needed for this, is of the order of ∆t ∝ R A , and is much smaller than the lifetime of all faster partons in the fluctuation. Hence, this interaction can be viewed as a rapid quench of the entangled partonic state [8] with substantial entanglement entropy. After this rapid quench, the interaction of the gluons change the Hamiltonian. Since in the CGC/saturation approach all partons with rapidity larger than that of a particular gluon y i , live longer than this parton, so they can be considered to be the source of the classical field that emits this gluon. It has been shown that after the quench, the fast gluons create the longitudinal chromo-electrical background field. Moving in this field the gluon accelerates and emits gluons which have the thermal distribution ( the first term of Eq. (2) below) [9][10][11]. The temperature of this distribution is intimately related to the saturation momentum, which provides the only dimensional scale in the colour glass condensate. It determines both the strength of the longitudinal fields and the ultraviolet cutoff on the quantum modes, resolved by the collision. It turns out [9][10][11][12][13] that with the semi-classical estimates [10] for the constant c = 1.2.
The appearance of a thermal emission in a high energy proton-proton collision is a remarkable feature of the interaction, since the number of the secondary interactions in proton-proton collisions is rather low, and cannot provide the thermalization due to the interaction in the final state. The origin of thermal radiation in the framework of the CGC approach was clarified a decade ago [9,10,12,13] and, recently, the new idea that the quantum entanglement is at the origin of the parton densities, has been added to these arguments [8].
The goal of this paper is to re-visit the process of inclusive production in the CGC/saturation approach, for more thorough consideration, and to show that the thermal term with the temperature given by Eq. (1), is needed for describing the experimental data for ion-ion collisions at high energies. At first sight, it looks as we are pushing at an open door, since it has been shown [11,[16][17][18][19][20] that the experimental data [6,[21][22][23][24][25] at high energy both for hadron-hadron and ion-ion scattering , can be described as the sum of two terms: The first term in Eq. (2) is the desired thermal radiation, while the second term describes hadron production in the hard processes, showing the power-like decrease at high p T . The same dependence on energy of both T th and T h supports the main idea of CGC, i.e. both parameters are related to the saturation momentum (Q s ). However, it turns out that this dependence differs from the saturation scale Q s ∝ s s0 λ .
The value of λ can be calculated theoretically and measured experimentally. The leading order QCD evaluation leads to λ = 4.9ᾱ S , whereᾱ S denotes the running QCD coupling. Plugging in a reasonable estimate forᾱ S (Q s ) ≈ 0.2, λ turns out to be large, about 0.8-1. The phenomenological description of the hard processes both for nucleus interactions [14] and DIS( see Ref. [15] and references therein), give the value of λ = 0.2 − 0.24. Thus, in the CGC approach we expect Hence, in spite of the fact that Eq. (2) and Eq. (3) show that both temperatures have dependence on energy in accord with the CGC result, this dependence contradicts the CGC prediction. Especially, the different Q s energy dependence in the second term of Eq. (2), which corresponds to the contribution of the hard processes, looks strange if not wrong, since the CGC approach to the hard processes has been confirmed both theoretically and experimentally, and we know that the typical scale in these processes, is the saturation momentum. The second remark is related to the value of the hard contribution. In the CGC approach it should be calculated theoretically, and not be determined by a fitting procedure.

A. General formulae
The general formula for the gluon jet production in ion-ion collisions in the CGC/saturation approach, has the following form (see Ref. [27] for the proof): where N Ai G (y 1 = ln(1/x 1 ); r, b) can be found from the amplitude of the dipole-nucleus scattering N Ai (y i = ln(1/x i ); r ; b) : where r denotes the size of the dipole, b it's impact parameter and where y denotes the rapidity of the produced gluon in c.m.f. and W the c.m.s. energy of the collision. In this paper we consider the gluon production at y = 0. C F = (N 2 c − 1)/2N c andᾱ S = α S N c /π with the number of colours equals N c . α S denotes the running QCD coupling, ∇ 2 T the Laplace operator with respect to r, it is equal to ∇ 2 T = 1 r d dr r d dr . In our paper [2] we found that the main contribution at high energies stem from the specific kinematic region in the vicinity of the saturation scale. Indeed, at high energies and sufficiently small values of p T the dipole amplitudes are in the saturation region, where the parton densities are large and the dipole scattering amplitude displays geometric scaling behaviour, being a function of only one variable: τ = r Q s (W, b) [28][29][30]. Deep in the saturation region the dipole amplitude tends to approach 1, but ∇ 2 Consequently the main contribution in Eq. (4) stems from the kinematic region where τ ∼ 1 or, in other words, from the vicinity of the saturation scale. The most attractive features of this observation is the fact that we know theoretically the behaviour of the scattering amplitude in this region. It has the form: whereγ = 1 − γ cr and γ cr = 0.37 in the leading order is the solution to the equation [1]: χ (γ) is given by Note that ω =ᾱ S χ (γ) is the eigenvalue of the BFKL equation [32]. The saturation momentum Q s (A, Y, b) has the following dependence on energy [3][4][5] where λ =ᾱ S χ (γ cr ) /(1 − γ cr ) in the LO of perturbative QCD. We will discuss the dependence of the saturation momentum on A below. It should be stressed that Eq. (7) together with Eq. (10) give the correct behaviour N ∝ exp (−µb) of the scattering amplitude at large impact parameters, which is determined by the non-perturbative behaviour of the saturation momentum at the initial energy (Y 0 ). It should be compared with the general case, for which at the moment we are not able to modify the main equations in a such way to obtain this behaviour [33]. Note, that we can find N G of Eq. (5) only if we know the impact parameter behaviour of the scattering amplitude. Hence, we conclude that ∇ 2 T N A1 G (y 1 = ln(1/x 1 ); r, b) and ∇ 2 T N A1 G (y 2 = ln(1/x 2 ); r, b) in Eq. (4), as well as the integral over r, can be calculated in the framework of CGC/saturation approach, and there is no need to introduce the non-perturbative corrections due to the unknown physics at long distances (see Refs.[67? ] for example) in the dipole scattering amplitude. However, plugging the scattering amplitude of Eq. (7) in Eq. (4), one can see that This divergence can only be tamed in the process of hadronization, which has to be treated using a non-perturbative approach to QCD.
At the moment, as has been alluded to in the introduction, we can treat the hadronization only using phenomenological models. Our model consists of two elements. First, every gluon decays into the jet of hadrons with a known fragmentation function: We take the fragmentation function D π G from Ref. [36] which has the form with α = 0.899, β = 1.57 and γ = 4.91. It should be stressed, that we assume that Eq. (11) holds for any value of the transverse momenta of gluons, while it is shown that this equation gives a good description of the experimental data for rather large transverse momenta: p T > 3 GeV [7] and p T > 5 GeV [38].
We note that Eq. (4) as well as Eq. (11) lead to a cross section which is proportional to 1/p 2 T . This behaviour results in a logarithmic divergency of the integral over p T , or in other words gives an infinite number of produced pions at fixed rapidity. The reason for this problem, is that we neglected the mass of the jet of hadrons that stems from the decay of the gluon. The simple estimates [39] give for a gluon with the value of the transverse momentum p T , the mass of the jet m 2 jet = 2p T m eff , where m eff = m 2 + k 2 T + k 2 L − k L , m is the mass of the lightest hadron in the jet, k T is it's transverse momentum and k L ≈ k T is the longitudinal momentum of this hadron. Since most pions stem from the decay of ρ-resonances we expect that m eff ≈ m ρ .
As we have mentioned in the introduction, our model for confinement is that of the CGC approach, the typical momentum for the produced gluon is the saturation momentum. Hence, most hadrons are created in the jets with the mass m 2 jet = 2Q s m eff . However, for rare gluons with p T Q s we still have m 2 jet = 2p T m eff . For numerical estimates we use m 2 ) m eff which has these two limits. Θ(x) denotes the step function. Using the same idea we replace Eq. (6) by Concluding, we see that our model of the hadronization processes includes two ingredients: the fragmentation function of Eq. (12) which gives us the number and p T distribution of the produced hadron from one gluon; and the replacement 1/p 2 T in Eq.

B. Dipole-nucleon scattering amplitude
For completeness of presentation, in this subsection we give a brief review of our theoretical and phenomenological approach to the contribution of the dipole-nucleon scattering amplitude to the inclusive gluon production, that has been discussed in Ref. [2].
In the vicinity of the saturation scale the dipole-nucleon scattering amplitude has a general form of Eq. (7) which leads to Inside of the saturation domain we suggest using N = 1 − exp −φ 0 τ 2γ which gives We need to match this formula with Eq. (14) at τ = 1. For N 0 1 we obtain φ 0 = N 0 . Ref. [2], checked that Eq. (15) describes to a good accuracy the exact solution for the non-linear Balitsky-Kovchegov [40] equation, for the leading twist BFKL kernel. We wish to remind the reader that the scattering amplitude in this kinematic region only gives a small contribution (not more than 15%) event at p T = 0.
In this region we can safely use the perturbative QCD approach for the scattering amplitude. Therefore, we need to solve the BFKL evolution equation in this region which has the following form: where N (Y ; x 01 , b) denotes the dipole scattering amplitude. x 01 = x 1 − x 0 ≡ r the size of the dipole. The kernel K (x 01 ; x 02 , x 12 ) describes the decay of the dipole with size x 01 into two dipoles of size: x 02 and x 12 = x 01 − x 02 . After integrating Eq. (16) over b the equation reduces to the BFKL equation [32] with the eigenfunctions r 2 γ . Therefore, the general solution takes the form: where χ (γ) is given by Eq. (9) and ξ = ln In the definition of ξ we introduce a new momentum scale which characterizes the value of the initial condition. For simplicity, we have that We can view this momentum as the saturation momentum at In the vicinity of the saturation scale, the integral over γ can be evaluated using the method of steepest descent, with the equations for the saddle point γ SP ≡ γ cr : Dividing the first equation by the second one, we obtain the value forγ = 1 − γ cr which is the solution of Eq. (8).
The first equation gives the value of the saturation momentum Expanding Ψ (Y, ξ, γ) at γ →γ we obtain Plugging Eq. (23) into Eq. (20) and integrating over γ − γ cr , the resulting solution can be written in the form In the method of steepest descend, the saddle point for the integration turns out to be 4. Impact-parameter dependent CGC dipole model As we have mentioned, the advantage of Eq. (24) is that we can introduce the correct behaviour of the amplitude at large impact parameter, by imposing the phenomenological decrease in saturation momentum for large b, by writing it in the form: In the LO BFKL [1] λ =ᾱ S χ(γcr) γ . Parameters N 0 and Q 0 , as well as function S (b) should in future be taken from non-perturbative QCD calculations but, at the moment, has to be determined from a fit to experimental DIS data. We have two models [15,42] † on the market that describe the final set of the HERA experimental data on deep inelastic structure functions [45] . They have different forms for S (b): Ref [42] : The ansatz of Eq. (29) is preferable, since it leads to S , which is in accord with the Froissart theorem [46] . However, we choose Eq. (28) which allows us to do several integrations analytically. Using Eq. (30), ∇ 2 T N takes the following form, after integrating over b: Before discussing the details of our model, we would like to outline, which features of Eq. (30) stem from the theorem, and which from phenomenological assumptions. The expression for τ ≤ 1, as we have mentioned (see Eq. (7)) follows from the theory. However, the calculation of N G (see Eq. (5)) takes into account the term of the order N 2 . The corrections of the order of N 2 has to be estimated for τ < 1 and, in principle, they change Eq. (??). In Ref. [2] we discussed these corrections, but we do not take them into account, in Eq. (30) as we view N = N 0 τ 2γ as a phenomenological expression that describes the DIS data [15].
The fact that the impact parameter behaviour of the saturation momentum determines the b-dependence of the scattering amplitude, comes from theory, while the particular form and result of integration over b, stems from the model for S (b).
For τ ≥ 1 we have discussed the form of Eq. (30) in the previous section (see also Ref. [2]). The b integration is performed with the phenomenological S (b).
In our estimates we use the values of the parameters from Ref. [15](see Table 1). In this paper the HERA data were fitted in the wide range of Q 2 from 0.75 GeV 2 to 650 GeV 2 . The expression for Q s (x) in this model is taken in the form ‡ Q s (x) = 1 2 It should be noted, that the value of x from Eq. (13) even at W = 13 TeV, is about 10 −5 , which is in the region that has been measured at HERA.   [15], which we use for our estimates.

Inclusive production for pT Qs
We can significantly simplify our calculation for p T Q s . Indeed, for such large p T the integral over r in Eq. (4) is concentrated in the region τ 1, where we can safely use the simple expression of Eq. (14). The integral can be † Actually, the same set of the data was described in the model of Ref. [43], but this model does not include the correct behaviour deep in the saturation region [44], and we do not discuss it here. ‡ Note that we introduce the extra factor 1 2 in the definition of the saturation scale, since we use τ = r Qs, while in Ref. [15] τ is defined as τ = rQs/2. calculated explicitly and leads to the following We know that in the vicinity of the saturation scale the scattering amplitude in the momentum representation has the following behaviour: Therefore, from Eq. (32) we can determine the value of the constant in Eq. (33) from the value of N 0 . Eq. (32) allows us to take into account the violation of the geometric scaling behaviour , given by Eq. (24) and Eq. (25). We can make such a replacement directly in the momentum representation, since r ∝ 1/p T . However, we need to find the coefficient in front of p T and, perhaps, an additional constant. We calculate the average τ using the expression: The results of these estimates are shown in Fig. 1. One can see at atp T → 0 τ = 0.478 ≈ 1 2 while at largep T it is proportional to 1/(4p T ). Forp T = 0, or more generally for p T Q s , the typical distance turns out to be r = 1/(2Q s (x)), and for largep T it is of the order of 1/(4p T ). Hence, we suggest to use in Eq. (25) the calculated τ (p T ) forp ≤ 4 and 1/(4p T ) for τ ≥ 4 and to substitute in Eq. (32) : wherep T = p T /Q s . In Re. [2] it was shown that this replacement results in the γ eff > 1 at largep T . This dependence is essential for describing the experimental data. Indeed, the value of n in the hard term in Eq. (2) is n = 3.1. As we have seen above (see Eq. (32) for example) at large p T the inclusive cross section is proportional to 1/p 4γ eff T . For γ eff =γ it is impossible to obtain a decrease of about 1/p 6 T , as indicated by the data, while Eq. (35) makes such a description possible (see Ref. [2] for detail discussion).

C. Dipole-nucleus scattering amplitude
For the dipole-nucleus scattering amplitude in the vicinity of the saturation scale τ ∼ 1 we have the same behaviour as for the dipole-nucleon amplitude: The main question that we will try to answer in this section is how the value of Const in Eq. (36) is related to N 0 for the dipole-nucleon scattering amplitude.
The principle difference between scattering with a nucleus and a nucleon is shown in Fig. 2. For dipole-nucleon scattering the initial condition at Y = Y 0 ≡ Y min is such that N (Y = Y 0 , r, b) 1(see Fig. 2-a). On the other hand, for dipole-nucleus scattering even at Y = Y 0 the shadowing corrections are large, and we impose the initial conditions inside the saturation region (see Fig. 2-b). These initial conditions lead to different solutions for dipolehadron and dipole-nucleus amplitude in the saturation region (see Refs. [47][48][49][50]). In particular, for ξ < 0, where we expect geometric scaling behaviour of the scattering amplitude, and for ξ > 0 and no such behaviour is seen. However, a glance at Fig. 2, shows that for Y − Y min 1 for ∇ 2 T N A (Y, r, b) we do not expect violation of the geometric scaling behaviour of the scattering amplitude, since, as has been discussed in the previous section, this observable gives the main contribution in the vicinity of the saturation scale shown by red line in Fig. 2. As we have discussed for the case of the dipole-nucleon scattering, the behaviour of the scattering amplitude in the vicinity of the saturation scale is determined by the solution to the linear BFKL equation (see Eq. (16)).   Fig. 2-a shows the kinematic regions for the dipole-nucleon amplitude, while in Fig. 2-b, we show the kinematic regions for dipole-nucleus scattering. Note that the saturation domain in this case can be divided in two subregions: (i) for ξ < 0, where we expect geometric scaling behaviour of the scattering amplitude, and ξ > 0 where there is no such behaviour. ξs = λ (Y − Ymin) and ξ = − ln τ 2 = ln r 2 Q 2 s . z = ξs − ξ and x = ξ + ξ.
Therefore, to find the value of Const in Eq. (36), we need to solve the BFKL equation with the correct initial condition. This condition is given by the McLerran-Venugopalan formula [5], which we use in the simplified form; The general solution to the BFKL equation of Eq. (16) has the same form as in Eq. (18): where Q s (Y 0 ) has been introduced in Eq. (19). We recall that ξ = − ln r 2 Q 2 s (Y 0 ) . Using Eq. (37) we can find n A in γ, we takes the following form: Plugging Eq. (39) into Eq. (38) we obtain the following solution: where ξ A = − ln r 2 Q 2 s (A, Y 0 , b) . Solving Eq. (21), we obtain the solution in the vicinity of saturation, in the form of Eq. (36) scale with In Fig. 3 we plot the value of R as function of q = Q 2 s (A, Y 0 , b) Q 2 s (Y 0 ). One can see that R is very close to 1. The chosen range of q we will be discussed in the next section. In further estimates, we take Const = N 0 . As we discussed in the previous section, Eq. (35) which takes into account the violation of the geometric scaling behaviour, is essential for the description of the experimental data, since the region of perturbative QCD gives a large contribution. Bearing this in mind, we take into account the contribution of n A in γ, into effective γ eff as well as Eq. (35). Considering the resulting Ψ in Eq. (20) in the form we obtain the following expression for the effective γ: where κ is defined in Eq. (35). Recall that τ = r Q s (A, Y ). The value of D A is plotted in Fig. 3-b, as a function of In the region of q = 1.5 − 4 D A ≈ 4, and we will argue below that this region contributes to the ion-ion collisions at the LHC energies. Fig. 4 shows the interaction of two nucleus A 1 and A 2 . From this picture we see that we need to know how many nucleons interact, and how this interaction occurs. In particular, one nucleon can interact with several nucleons from another nucleus, and so on. We use the KLN model to answer all similar questions (see Refs. [14,39,[52][53][54][55]). We choose this model since, on one hand, this model not only successfully described the RHIC data and it's predictions for the experimental data at the LHC, are quite good. On the other hand, it is based on the Glauber approach [56,57], which successfully describes the nucleus-nucleus interaction over a wide energy range. In Refs. [52,58] the Glauber approach is developed for the hadron production for nucleus-nucleus collisions, and it demonstrates how we can calculate (i) the number of nucleons that participate in the interactions (number of participants N part ), their density ρ part at fixed values of b; and (2) how to correlate the typical values of b with the centrality c, which is measured experimentally, and which characterizes the percentile of events with the largest number of produced particles (as registered in detectors). Experimentally collisions are grouped into event(centrality) classes, with the most central class defined by events with the highest multiplicity(smallest forward energy), which corresponds to small values of the impact parameter. Basically centrality c (N ) is equal to In our approach, we use the estimates of Ref. [52] for the density of the participants at for different centrality classes § .

III. KLN MODEL FOR ION-ION SCATTERING AT HIGH ENERGY
The key idea of the KLN model [39,52] is, that the saturation momentum is proportional to the density of the participants, Q s ∝ ρ part . The arguments for this stem from the simple equation for the saturation momentum[1, 3-5] where Y = ln (1/x), N c denotes the number of colours, R A the radius of the nucleus and xG A is the gluon structure function of the nucleus. For the kinematic region where the gluon densities are small we can safely consider G A Q 2 , x = A xG N Q 2 , x , where G N is the gluon structure function for a nucleon. Plugging this relation in Eq. (45) as well as R A = R N A 1/3 , § For comparison with the ALICE experimental data we use the results of the estimates, given by Ref. [59], in which the procedure of Ref. [52] is used in the Glauber Monte Carlo.
where R N is the radius of the nucleon, we obtain that where ρ A is the density of the nucleons in the nucleus in the transverse plane. The KLN model generalizes Eq. (46) proposing We suggest to re-write Eq. (47) in the form Factor 1 2 reflects the fact that we are dealing with the density of those nucleons in a single nucleus, which will participate in the collision at a given impact parameter b, or in a definite centrality class. In Eq. (48) we write the expression for the nucleon saturation momentum in the frame.
One can see that for the DIS with a nucleus with N part = A we obtain from Eq. (48) in accord with Refs. [60][61][62][63]. Using the estimstes of ρ part from Ref. [52] we obtain that the range of q in Fig. 3 for the lead-lead scattering is q = 1 − 3.06.
For the gluon transverse momenta p T distributions, we suggest the following formula for the ion-ion collisions in the definite centrality class with the number of participants N part : where N G is related to the scattering dipole-nucleon amplitude (see Eq. (5)), and n G denotes the multiplicity of the emitted gluons. σ in denotes the inelastic cross section for the nucleon-nucleon scattering at corresponding energy. One can see, that in Eq. (50) we consider each participant as a nucleon which interacts with another nucleon from the different nucleus, at impact parameter b (see Fig. 4). The saturation momentum for this scattering is given by Eq. (48), for the effective γ eff we take Eq. (43), which includes corrections from the interaction with the other nucleons. We need also to fix the dependence of the saturation scale on the impact parameter of the nucleon-nucleon interaction (b in Fig. 4). Finally, the equation for the saturation scale, which we use in our estimates of gluon production in nucleon-nucleon scattering, takes the form: where Q s and S (b ) are given by Eq. (27) and Eq. (28), respectively. At first sight, we do not need Eq. (50), since we can use the master equation ( see Eq. (4)) and using the experimental information of gluon structure functions of nucleus (see for example Ref. [64]), we can calculate the inclusive cross section as we did for proton-proton scattering in Ref. [2]. Indeed, for the inclusive production of gluons for an ion-ion collision, we can proceed in this manner, but the factorization of Eq. (4) is proven only in the case, when we did not make any additional selections of the events summing over all accompanying hadrons. Considering the different centrality classes we make the additional selection on the multiplicity of produced hadrons. It is sufficient to refer to the AGK cutting rules [65], to see that these selections violate the factorization of Eq. (4). The violation of the AGK cutting rules in QCD (see Ref. [1]) makes the situation even worse. Hence, the KLN model gives us the simple way to estimate the cross sections for the different centrality classes.

IV. COMPARISON WITH THE EXPERIMENT
We calculate the cross section for gluon production using Eq. (50). As we have discussed in section II-B, the typical values of r that contribute to the integral in Eq. (50) are rather small, of the order of Q s . This means that we can safely use the CGC/saturation approach or/and the perturbative QCD estimates for this integral. We do not need to incorporate modifications of the gluon propagators, for example due to the confinement [66,67], in our calculations.
The factor 1/p 2 T in front in Eq. (50) stems from the gluon propagator [32], and it is affected both by the hadronization, and by interactions with co-movers in the parton cascade. In our approach to the confinement problem, we first need to take into account, the effect of the mass of produced gluon jet due to hadronization, which changes the gluon propagator [39]: Therefore, we calculate gluon production using Eq. (50), in which we use Eq. (52) to replace the factor 1/p 2 T . Therefore, our model of the hadronization consists of two ingredients: the decay of the gluon into hadron jet with the fragmentation function (see Eq. (11)), and to account for the mass of this jet, using Eq. (52). It is instructive to note, that in our approach we see two stages of the multi-particle production in an explicit way: the creation of the Colour Glass Condensate (calculation of the integral over r in Eq. (50)) and the stage of hadronization ( Eq. (52), and fragmentation functions).
For calculating the integral over r in Eq. (50) we use Eq. (30) for p T ≤ 2 Q s (x) and Eq. (32) for larger values of p T (see Ref. [2] for more details). The values of σ in for proton-proton scattering at high energies we take from Ref. [68], which describes all available data on measured soft cross sections.
We compare with the experimental data of ALICE collaboration [69][70][71][72] on the tansverse momentum distribution of the charged hadrons in different centrality classes at two energies W = 5 TeV and W = 2.78 TeV. We add to Eq. (50) the thermal term (see first term of Eq. (2)), and find that in our model of hadronization we need this term to describe the experimental data. The value of this contribution essentially depends on the mass of the gluon jet in the framework of our model. It turns out that the temperature T th is proportional to the value of the saturation scale in the given centrality class in accord with Eq. (1), but the coefficient c turns out to be 1.5 larger than it is predicted in Ref. [10] (c = 1.8).
In Fig. 5, Fig. 6 and Fig. 7 we show that we describe the data fairly well in the region p T ≤ 10 GeV . It should be noted that we do not need any K-factor which would account for the higher order corrections in the framework of CGC/saturation approach. However, this conclusion we need to take with a degree of scepticism, since the value of σ in was taken from the model [68], and the value ofᾱ S in Eq. (50) was takenᾱ S = 0.25. The main uncertainty in σ in = σ tot − σ el − σ dif f is the value of σ dif f , which at high energies is about 15 -20% of the total cross sections. The value ofᾱ S Q 2 s ≈ 0.3 − 0.4. Therefore, these two uncertainties can lead to the K ≈ 2. Fig. 5, Fig. 6 and Fig. 7 show that we describe the behaviour of the transverse momentum distribution at p T ≥ 2 by Eq. (50), rather well. We wish to note that this description stems from the expression for γ eff (see Eq. (43) ), in which the last term comes from the corrections related to the nucleus target. In the description of p T distribution, these corrections are considerable.
The rate of thermal radiation is shown in Table  2, Note that the contribution of the thermal radiation increases with the growth of energy. The value of the CGC term depends on the value of the m eff . We believe that most of the pions are produced from ρ resonances and we consider m eff = 0.5 GeV . We recall that the simple formula for m eff = µ 2 + k 2 T + k 2 L − k L leads to m eff = 0.5 GeV if µ is equal to the mass of ρ-resonance since the value of k T = k L = 0.45 GeV (see Ref. [73] for the measurement and Ref. [74]) and reference therein for theoretical discussions). For the minimal mass of µ = m π = 0.14 GeV we obtain m eff = 0.2 GeV . To illustrate the influence of the mass of gluon jet we estimate the contribution of Eq. (50) with m eff = 0.14 M eV . One can see that even such small mass cannot describe the spectrum without the thermal radiation term.
In Fig. 8 we plot the estimates for the p T distribution with m eff = 0. It turns out that we need to add the thermal emission with R = 42%. This should be compared with the proton-proton scattering [2], where for m eff = 0 we do not need to add the thermal emission.
In general, comparing table II with the estimates for the ratio R for proton-proton scattering (see Table II of Ref. [2]) we see that the contributions of the thermal radiation are larger for ion-ion collisions for the classes with small centrality. It supports the CGC picture in which the longitudinal fields, that are the sources of the thermal radiation, are stronger in the denser gluon states which are produced in ion-ion collisions. Discussing hadron production, we have to construct a model for the process of hadronization. Our model where the production of the gluon jets with the hadronization, which is given by the fragmentation functions and by the mass of the gluon jet, requires a thermal radiation term. Our estimates for m eff = 0 show that, possibly, our claim does not depend on the model of the hadronization in the case of the ion-ion collisions, since the description of the experimental data in Fig. 8-a requires a thermal radiation term of the order of 46%. In Fig. 8-b we plot the p T spectrum for a different model, with the gluon propagator 1/(p 2 T + m 2 ) with m = T th , instead of Eq. (52). It turns out that R = 46% is needed to describe the data. The last case corresponds to a different hadronization model : the propagator of the gluon with transverse momentum p T in the CGC medium with the temperature T th , acquires a mass m g ∝ T th [75] and the propagator has the form 1/(p 2 T + m 2 g ). This mass provides the infrared cutoff in the gluon spectrum. However, we need to take our estimates in Fig. 8-b with a grain of salt, since Ref. [75] predicts that the gluon mass will be m = g T but with small value of g. However, Fig. 8-a supports our claim, that the requirement of the thermal emission does not depend on the details of the assumption on the confinement of quarks and gluon.

V. CONCLUSIONS
The main result of the paper is, that we show that a thermal emission term is required to describe the transverse momenta distribution for charged particles in ion-ion collision. The temperature of this emission T th turns out to proportional to the saturation scale (see Eq. (1)) with coefficient c = 1.8, which is 1.5 times larger than predicted in Ref. [10].
We develop the formalism for the calculation of the transverse momenta spectra in CGC/saturation approach, in which we clearly see two stages of the process: the creation of colour glass condensate; and the hadronization stage. Our calculations are based on the observation that even for small values of p T the main contribution in the integration over r in Eq. (4) and in Eq. (50) stems from the kinematic region in the vicinity of the saturation momentum, where theoretically, we know the scattering amplitude. In other words, it means that we do not need to introduce the non-perturbative corrections due to the unknown physics at long distances (see Refs. [66,67] for example) in the dipole scattering amplitude. The non-perturbative corrections have to be included to describe the process of hadronization, which we discuss in the model. This model incorporates the decay of the gluon jet with the effective mass m 2 eff = 2Q s µ soft where µ soft is the soft scale, and with the fragmentation functions of Eq. (12), at all values of the transverse momenta.
It should be emphasized that we reproduce the experimental data without any K-factor, which is used for accounting of the higher order corrections. We wish to mention, that we have calculated the inclusive production takingᾱ S = 0.25. This value is less thatᾱ S (Q s ) = 0.3 − 0.4 which appears more natural in Eq. (50). Forᾱ S =ᾱ S (Q s ), we need to introduce a K-factor of about 1.3 -1.6.
We use the KLN model [14,39,[52][53][54][55] which gives us the simple way to estimate the cross sections for the different centrality classes. At first sight, we do not need Eq. (50), since we can use the master equation ( see Eq. (4)), and using the experimental information of gluon structure functions of nucleus (see for example Ref. [64]) we can calculate the inclusive cross section, as we did for proton-proton scattering in Ref. [2]. However, considering the different centrality classes, we make the additional selection on the multiplicity of produced hadrons, which violates the factorization of Eq. (4). The KLN model suggests a way to evaluate the contribution of the different centrality classes.
scattering [2]. Therefore, we suspect that that the existence of the thermal radiation term does not depend on the model of confinement.