Thermal radiation and inclusive production in the CGC/saturation approach at high energies

In this paper, we discuss the inclusive production of hadrons in the framework of CGC/saturation approach. We argue, that the gluon jet inclusive production stems from the vicinity of the saturation momentum, even for small values of the transverse momenta $p_T$. Since in this region, theoretically, we know the scattering amplitude, we claim that we can provide reliable estimates for this process. We demonstrate, that in a widely accepted model for confinement, to describe the experimental data, we require a thermal radiation term. In this model the parton (quark or gluon) with the transverse momenta of the order of $Q_s$ decays into hadrons with the given fragmentation functions. However, we show that other approaches for the confinement could describe the data, without a need for the thermal emission


I. INTRODUCTION
In this paper we discuss the dynamics of generating multi-hadron processes at high energy in the framework of the Color Glass Condensate(CGC)/saturation approach (see Ref. [1] for the review). These processes occur at long distances and therefore to treat them theoretically, we need to develop a non-perturbative QCD approach. This is a very difficult and challenging problem, which is far from being solved. The CGC/saturation approach, being an effective QCD theory at high energies, states that the new phase of QCD: the dense system of partons (gluons and quarks) is produced in collisions with a new characteristic scale: saturation momentum Q s (W ), which increases as a function of energy W [2][3][4]. However, the transition from this system of partons to the measured state of hadrons is still an unsolved problem. At the moment, we need to use pure phenomenological input for the long distance non-perturbative physics, due to our lack of theoretical understanding of the confinement of quarks and gluons. 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 with the given fragmentation functions. The experimental data confirm this model of hadronization, which is 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. [5], which shows that the next-to-leading order QCD calculations with formation of the hadrons in accord with the fragmentation functions ( [6]), 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 multi -hadron production.
It turns out [7][8][9][10][11][12] that the experimental data [5,[13][14][15][16][17] at high energy, can be describes as the sum of two terms: GeV; We believe, on a qualitative level, that these two terms have a natural explanation in the CGC/saturation approach. The second term has a power-like decrease (∝ 1/p 2n T ) at large p T , as it should be in perturbative QCD. The appearance of a thermal term 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 [18][19][20][21] and, recently, the new idea that the quantum entanglement is at the origin of the parton densities has been added to these arguments [22]. The resulting picture is presented nicely in Ref. [11], to which we refer our readers. The brief sketch below, is intended to indicate the main ideas that originate in the CGC approach.
For proton-proton scattering in the lab. frame, the parton configuration in QCD is formed long before the interaction at distances 1/(mx), where m -denotes the proton mass, and x the fraction of longitudinal momentum carried by parton which interacts with the target. However, before the collision, the wave function of this partonic fluctuaction is the eigenfunction of the Hamiltonian and, therefore, the system has zero entropy. The interaction with the target of size R 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, 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 [22] with substantial entanglement entropy. After this rapid quench, the interaction of the gluons change the Hamiltonian. In the CGC approach, all partons with rapidity larger than that of a particular gluon y i , live longer than this parton. They can be considered as the source of the classical field that emits this gluon. It was shown that after the quench, the fast gluons create the longitudinal chromo-electrical background field, which leads to the thermal distribution of the produced gluons.
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 [18,21] that with the semi-classical estimates [21] for the constant c = 1.2. The saturation scale Q s depends on x and the impact parameters (b) of the reaction, and has the following form [2][3][4]: 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 the reasonable estimate forᾱ S (Q s ) ≈ 0.2, one can see that λ turns out to be large, about 0.8-1. The phenomenological description of the hard processes both for nucleus interactions [23] and DIS( see Ref. [24] and references therein), give the value of λ = 0.2 − 0.24. We therefore see, that in the CGC approach we expect that Eq. (1) and Eq. (2) show that both temperatures have dependence on energy in accord with the CGC result but, on the other hand, this dependence contradicts the CGC prediction. Especially, the second term in Eq. (1), which corresponds to the contribution of the hard processes looks strange. Indeed, the above interpretation of the thermal radiation cannot be considered as the conventional one, 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 from a fitting procedure. The goal of the paper is to re-visit inclusive production in the CGC/saturation approach for a more thorough consideration, and to show that the thermal term with the temperature given by Eq. (3), is needed for describing the experimental data at high energies. It should be noted that in the first attempt [25] to compare the CGC prediction with the experiment at W = 7 T eV , the thermal term was not required.
The paper is organized as follows. In the next section we discuss the general procedure for the calculation of the gluon inclusive production in CGC/saturation approach. In section III, we consider the evolution equation for the theory with a simplified BFKL kernel. We show that the solution to this equation confirms our key idea, that the main contribution to the inclusive production stems from the kinematic region in the vicinity of the saturation scale. Since theoretically we know the scattering amplitude in this region, we demonstrate that we are able to provide reliable estimates for this process. In section IV we develop the saturation model which we need to use due to the long standing unsolved problem i.e. the behaviour of the scattering amplitudes at large impact parameter. In section V we compare our estimates with the experimental data, and demonstrate that within our model for confinement: the fragmentation function for the gluon jets, needs to have a thermal radiation term with temperature, which is proportional to the saturation scale Q s . In the Conclusions we summarize our results.

II. INCLUSIVE PRODUCTION IN CGC/SATURATION APPROACH: GENERALITIES
The formula for the gluon jet production in proton-proton collisions has the following general form (see Ref. [26] for the proof): where N G (y 1 = ln(1/x 1 ); r, b) can be found from the amplitude of the dipole-proton scattering N (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 . 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). Introducing a new variable z (r, b, x) = ln τ 2 , we can re-write Eq. (5) at In Eq. (8) we have taken x = x 1 = x 2 since y = 0, and introduce a new variablep T = p T /Q s (x).
To calculate the hadron distributions, we need to take into account the decay of the jet into hadrons. The formula has the form We take the fragmentation function D π G from Ref. [30] which has the form with α = 0.899, β = 1.57 and γ = 4.91.
We note that Eq. (8) as well as Eq. (9) leads 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 [31] 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 jet = 2 (Q s Θ(Q s − p T ) + p T Θ(p T − Q s )) m eff which has these two limits. Θ(x) denote the step function. Using the same idea we replace Eq. (7) by

A. Equations
The dipole scattering amplitude is the solution of the Balitsky-Kovchegov [32] non-linear equation, which has the following form: 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 .
The linear part of this equation reduces to the BFKL equation [33] with the eigenfunctions r 2 γ which corresponds to the value of the kernel: where ψ(z) = d ln Γ(z)/dz is the Euler ψ-function (see Ref. [34] formula 8.36).
At first sight, we need to solve this equation to obtain the inclusive cross section. However, we need to know the dependence of the dipole amplitude on the impact parameter, which cannot be found from Eq. (12). The failure to reproduce the correct large b behaviour of the scattering amplitude, is a long standing problem of non-perturbative QCD contributions [35] to the BK equation. As a consequence we are doomed to use a phenomenological input in addition to Eq. (12). Hence, we suggest the following strategy: to use a simplified form of Eq. (12) and to study in this approach the main features of the inclusive production. After such an investigation, we will select the model which satisfies both the BK equation, and reproduces the correct behaviour at large b.
The BFKL kernel of Eq. (13) includes the summation over all twist contributions. In the simplified approach we restrict ourselves to the leading twist term only, which has the form instead of the full expression of Eq. (13).
As indicated in Eq. (14) we have two types of logs: To sum these logs it is necessary to modify the BFKL kernel in different ways in the two kinematic regions, as shown in Eq. (14).
• τ = r Q s 1 For the perturbative QCD region of τ 1, the logs originate from x 2 02 ∼ x 2 12 x 2 01 ≡ r resulting in the following form of the kernel K ( The non-linear BK equation in this region can be written as 1 Inside the saturation region where τ > 1 the logs originate from the decay of a large size dipole into one small size dipole and one large size dipole. However, the size of the small dipole is still larger than 1/Q s . This observation can be translated in the following form of the kernel Inside the saturation region the BK equation takes the form . The advantage of the simplified kernel of Eq. (14) is that, in the Double Log Approximation (DLA) for τ < 1, it provides a matching with the DGLAP evolution equation [37].
B. Solutions

Perturbative QCD: linear equation
First we discuss the solution to Eq. (16) in the perturbative QCD region, where one can neglect the contributions of the shadowing corrections. The equation has the form: One can recognize that Eq. (19) is the DGLAP equation [37] in the double log approximation (DLA). It has the following solution (see for example Ref. [1]) where we use the following notation: The solution of Eq. (20) provides the boundary condition for the solution inside the saturation region: As was expected [29], in the vicinity of the saturation scale ( ζ 8ξ s ), the amplitude exhibits geometric scaling behavior, being a function of only one variable ζ and it has the following form [38]: where γ cr denotes the critical anomalous dimension which in the DLA is equal to 1 2 .

Solution in the region τ < 1
The solution of Eq. (23) assumes that the value of the constant N 0 in this equation is small, and we can neglect the non-linear term. However, we need to estimate the non-linear corrections in this region, due to the value of N G , which we need to evaluate for the calculation of the inclusive production (see Eq. (6)), explicitly accounts for the N 2 term. Eq. (16) can be simplified assuming the geometric scaling behaviour of the amplitude in the vicinity of the saturation scale [29]. It has the form: We can obtain the solution to Eq. (24) introducing dn/dz = p (n) . For this function we obtain the equation: which has the solution Since p ∝ n for small n, as we have seen above, we conclude that C 1 = 0. The amplitude n can be found by solving the algebraic equation: The first correction to the scattering amplitude of order n 2 has the following form: 3. Solution in the region τ > 1 In the saturation region it has been shown [27] that the scattering amplitude manifests geometric scaling behavior, which is supported by the experimental data [28]. Therefore, Eq. (18) can be re-written in the form: Introducing we re-write Eq. (29) in the form: Integrating we obtain which can be resolved as Taking the derivatives on both sides of Eq. (33) we reduce this equation to the form: Therefore, we can find φ as the solution to the following equation: The value of C has to be found from matching with the region τ < 1. For small φ 0 C = 0. Indeed, in this case the solution at small φ has the following form: which coincides with Eq. (23) for the region τ < 1 at small φ 0 = N 0 1. For the values of N 0 , which are not very small, we need to find the value of C from the matching conditions, that can be taken from Eq. (28):

Inclusive production
In our approach we can estimate the value of function I (p) in Eq. (8). Taking N 0 = 0.05, and considering this as a small value, we can use Eq. (35) with C = 0. In Fig. 1-a we present the results of the numerical estimates for the amplitude and its derivatives. One can see that N (z) tends to 1 at large z, but both N z and N zz become very small at large z. The integrand in Eq. (8) for the inclusive production is equal to e −z N zz (z) and it is shown in Fig. 1-b. From this figure, we note that the main contribution stems from z < 0, where our amplitude has the simple form of Eq. (23). It turns out that the behaviour for z > 0 can be approximated by the simple formula for the scattering amplitude (see Fig. 1-b): Hence, one can see that the situation for inclusive production looks quite different to that for the description of the DIS structure function. ∇ 2 T in Eq. (8) generates an extra factor e −z which is large at z < 0, and enhances the contribution of the perturbative QCD region. In principle, we expected that for p T < Q s the region of large distances will contribute, and the physical observable will be sensitive to the theoretical expectations in this region. Fig. 1-b shows that it is not the case for inclusive production. To illustrate the influence of the saturation region (z > 0), we calculate the function I (p) at p T = 0. In I (p T = 0) we expect the largest contribution to come from the saturation region. We found, using the numerical solution for the scattering amplitude, that the contribution of the perturbative QCD region at z < 0 to I (p T = 0) gives 85%, while only 15% comes from the saturation domain with z > 0, for the range (−8 < z < 8). One can see that the integral diverges at z → −∞, so, we need to generalize the behaviour of φ in the perturbative QCD region (see Eq. (20)) by replacing Takingᾱ S = 0.2 we find that integration over negative z gives 85% of the total contribution. Hence, we conclude that the following expression provides a good approximation of the function I(z) in Eq. (8): The solution φ = φ 0 e 1 2 z is correct only at small values of N 0 . If N 0 is not very small (say is about 1/3) we need to take the integral over φ in Eq. (35) keeping C = −N 3 0 /2 in the dominator(see Eq. (37)). The solution for φ(z) at small z has the form One can see that φ 2 (z) is close to φ(z) = N 0 e 1 2 z , even for N 0 = 0.3 − 0.4.

IV. IMPACT-PARAMETER DEPENDENT CGC DIPOLE MODEL
The result of the assay in the previous section can be formulated as follows: in inclusive production the main contribution comes from the vicinity of the saturation scale, or in the region of perturbative QCD, while contributions from long distances can be neglected. The behaviour of the amplitude in the vicinity of the saturation momentum is predicted theoretically, and has the form [38]: whereγ = 1 − γ cr and γ cr = 0.37 in the leading order is the solution to the equation [1]: χ (γ) is given by Eq. (13). The advantage of Eq. (43) 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 [24,39] on the market that describe the final set of the HERA experimental data on deep inelastic structure functions [40] . They have different forms for S (b): Ref [39] : The ansatz of Eq. (47) is preferable, since it leads to S , which is in accord with the Froissart theorem [41] . However, we choose Eq. (46) which allow us to do several integrations analytically. Using Eq. (48) ∇ 2 T N takes the following form after integrating over b: Before discussing the details about our model, we would like to outline, which features of Eq. (48) stems from the theorem, and which from phenomenological assumptions. The expression for τ ≤ 1, as we have mentioned (see Eq. (43)) follows from the theory. However, the calculation of N G (see Eq. (6)) takes into account the term of the order N 2 . As we have demonstrated in section III-B-2 the corrections of this order appears in this region in the scattering amplitude and has to be included in calculation of N G . We did not take them into account, because we view N = N 0 τ 2γ as a phenomenological expression that describes the DIS data [24]. However, we check how large the contribution to the inclusive production from the corrections of Eq. (27) and Eq. (28) is. These corrections are shown by the dashed curve in Fig. 2. One can see that the value of such a correction is not more than 2%.
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 the integration over b, stems from the model for S (b).
For τ ≥ 1 we have discussed the form of Eq. (48) in the previous section, and have given strong arguments for such an expression. The b integration is performed with the phenomenological S (b).
In our estimates we use the values of the parameters from Ref. [24](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. (11) even at W = 13 TeV, is about 10 −5 , which is in the region that has been measured at HERA.   [24], which we use in our estimates.
In Fig. 3 we plot functions d 2 b ∇ 2 T N (z(r, y, b)) and I (τ ) of Eq. (8). One can see that the both functions provide the main contribution from the region τ < 1. The largest contribution to the function I (p = 0) in the region of integration τ > 1, is 13.5% . For large p T we expect that the contribution of the region τ > 1 to the function I (p T ) is small. In this region we have an analytical expression for this function: In Fig. 4 we plot the numerical result for I (p T ) with functions given by Eq. (48) and the analytical expression of Eq. (50). Note that forp T ≥ 2 (or p T ≥ 2 Q s (x) ) both functions coincide.
The dashed line in Fig. 4 illustrates the result of the analytical integration over τ , assuming that the region for τ > 1 does not contribute. The result of this integration has the form: * Note that we introduce the extra factor 1 2 in the definition of the saturation scale since we use τ = r Qs, while in Ref. [24] τ is defined as τ = rQs/2.  N (z(r, y, b)) ( Fig. 3-a) and I (τ ) of Eq. (8) (Fig. 3-b) versus τ in the model, given by Eq. (48). the values of the parameters are taken from   Fig. 4 we can conclude that the saturation region with τ ≥ 1, only provides around 10-15% of the contribution forp T ≤ 1, but the sharp cutoff changes the behaviour for large values ofp T . We know that in the vicinity of the saturation scale the scattering amplitude in the momentum representation has the following behaviour: Therefore, from Eq. (50) we can determine the value of constant in Eq. (52) from the value of N 0 . For large p T Q s we cannot use neither Eq. (22) nor Eq. (52), since we have to take into account the violation of the geometric scaling behaviour in perturbative QCD (see Eq. (20)). Eq. (20) indicates that part of this violation can be taken into account by replacing in Eq. (22): In the DLA Eq. (53) gives Eq. (40). Eq. (53) is used in Ref. [24] for fitting the HERA experimental data.
We 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. 5. 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 distances turns out to be r = 1/(2Q s (x)), and for largep T they are of the order of 1/(4p T ). Hence, we suggest to use in Eq. (53) the calculated τ (p T ) forp ≤ 4 and 1/(4p T ) for τ ≥ 4.
In Fig. 6 we show the behaviour of the γ eff at different energies. One can see that this dependence is essential for describing the experimental data. Indeed, the value of n in the hard term in Eq. (1) is n = 3.1. As we have seen above (see Eq. (50) 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. The p T and W dependence of γ eff of Eq. (53) is shown in Fig. 6. We see that even this behaviour of γ eff leads only to 1/p 4−5 T at p T ≈ 7 GeV. Hence, we expect that the CGC/saturation approach in the form of the model, will not be able to describe the hadron spectra at large p T . However, at large p T we are outside of the vicinity of the saturation scale, and have to perform the perturbative QCD calculation which, as we have mentioned, describes the experimental data [6]. It should be emphasized, that the discussed problems of the CGC approach, has no influence on the behaviour of the transverse momentum distributions at p T ≤ 7 GeV, and on the value of the contribution of the thermal radiation, which is negligibly small at p T ∼ 7 GeV.

V. COMPARISON WITH THE EXPERIMENTAL DATA
As we have discussed in the introduction, our goal is to answer the question. Do we need the thermal radiation term to describe the experimental data in the framework of the CGC/saturation approach? Our answer is yes. We calculate the cross section for gluon production using Eq. (8). As we have discussed in function I p T Qs the main contribution comes from p T ∼ Q s (see Fig. 4). However, the factor 1/p 2 T in front in Eq. (8) stems from the gluon propagator [33], 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 [31]: Therefore, we calculate gluon production using Eq. (8) in which we use Eq. (55) to replace the factor 1/p 2 T .  [24] has a different assumption regarding the behaviour of the amplitude in the saturation region. However, ∇ 2 N in this function is not continuous at τ = 1 and leads to a very small contribution for τ ≥ 1 (see Fig. 7). Hence, in the original model we neglect the contribution from the saturation region. As we have seen in Fig. 4 such a procedure results in 10-15% accuracy of the calculations in the entire region of p T . It should be mentioned that our assumptions about the behaviour of ∇ 2 N in the saturation region, are based on the solution to the BK equation in the vicinity of the saturation scale. While the model of Ref. [24] reproduces only the the theoretical behaviour of Ref. [36] at large τ or, in other words, in the region which does not contribute to ∇ 2 N . We note that the inclusive production calculated from the CGC/saturation approach, has a different form than the hard term of Eq. (1) (see Fig. 8).
We have partly explained that γ eff from Fig. 6 is not able to describe the shape of p T distribution at high p T . On the other hand, we predict the value of the cross section, while in Eq. (1) this value was a fitted parameter. We will discuss below the dependence of the cross section on the value of m eff . Fig. 9 illustrates how well we describe the data for the transverse momentum distribution at the LHC. T th ∝ Q S was taken from Eq. (3) , however, it turns out that the experimental data can be described with c = 2.3 which is almost twice larger than estimated in Ref. [21].
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 as the most reliable estimate. In Fig. 10 we present the calculation at W = 7 TeV with m eff = 0 and m eff = 0.06 GeV without the thermal emission term.
We see that at small values of the effective mass, we can describe the experimental data without the thermal radiation   term. It should be stressed that we do not need the so called K -factor, to include the next-to-leading order corrections. Even for the multiplicity distribution at W = 13 TeV we are able to describe the data using σ in = σ tot − σ el − σ dif f from Ref. [42]. 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. [43] for the measurement and Ref. [44]) and reference therein for theoretical discussions). For the minimal mass of µ = m π = 0.14 GeV we obtain m eff = 0.2 GeV .
Discussing hadron production we have to construct a model for the process of hadronization. Our model is the production of the gluon jets with the hadronization, which is given by the fragmentation functions. We showed that in this model for confinement, we obtained a reasonable description of the experimental data, with the thermal radiation and with the temperature of Eq. (3), which is predicted in the CGC approach. It is possible that our hadronization model is too primitive, and for the gluon with the transverse momenta of the order of Λ QCD , we should not apply the CGC formulae which are based on the perturbative QCD approach. If we cut our gluon spectra at p T = Λ QCD , we obtain a good description of the experimental data, without the thermal radiation. An alternate picture could be the following: 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 [45] and the propagator acquires the form 1/(p 2 T + m 2 g ). This mass provides the infrared cutoff in the gluon spectrum. Therefore, the same rescatterings in the produced CGC medium, which generates the thermal spectrum can be a reason for blocking the small gluon p T ≈ T th ≈ 0.12 − 0.14 GeV . In this picture we will not see any thermal emission in the spectrum of hadrons. Note, that p T ∼ 0.12 − 0.14 GeV corresponds to the m eff ≈ 0.06 GeV and to the distribution of Fig. 10-b.
In Fig. 11 we present the estimates with the gluon propagator 1/(p 2 T + m 2 ) with m = T th . One can see that we are able to successfully describe the data without the thermal radiation term.

VI. CONCLUSIONS
The main result of the paper is, that we show the need for thermal emission within a particular model for confinement: the parton (quark or gluon) with the transverse momenta of the order of Q s decays into the hadrons with the given fragmentation functions. The temperature of this emission turns out to be equal to 2.3/(2 π) Q s , as was expected in the CGC/saturation approach. Note, that the coefficient c in Eq. (3) turns out to be in almost two times larger than predicted in Ref. [21].
We develop the formalism for the calculation of the transverse momenta spectra in CGC/saturation approach, which is based on the observation that even for small values of p T the main contribution stems from the kinematic region in vicinity of the saturation momentum, where theoretically, we know the scattering amplitude. We suggest to  [5,13]. For the description of the data at W=13 TeV we use the model of Ref. [42], for σin = σtot − σ el − σ dif f . take into account the behaviour of ∇ 2 N in the saturation region in the form: and demonstrate that this suggestion follows from the solution to the non-linear Baltsky-Kovchegov equation, for the simplified BFKL kernel. 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 also to mention, that we have calculated the inclusive production takinḡ α S = 0.25. This value is less thatᾱ S (Q s ) = 0.3 which appears more natural in Eq. (8). Forᾱ S =ᾱ S (Q s ), we need to introduce a K-factor of about 1.3 -1.5. The value of the thermal radiation term contribution depends on the value of the effective mass m eff , however, in the region of possible values for this mass 0.12 − 0.5 GeV we need to account for the thermal emission to describe the spectrum at low p T .
We show that a different mechanism of confinement that blocks the emission of gluons with p T ≤ Λ QCD or/and that generates the gluon mass m = T th , is able to describe the experimental data without the thermal radiation.
In our approach we are able to evaluate the kinematic region that we can use the formalism of the deep inelastic structure functions. The structure function is related to the scattering amplitude being d 2 bN (Y, r, b). However, as we have discussed the inclusive production is determined by function N G (see Eq. (6)). In the region where we can neglect the N 2 -term in N G , we can safely perform the integration over b, and obtain the expression for the inclusive production through the structure functions. In Fig. 12 we show Eq. (50) and the first term of this equation which is the gluon structure function in the vicinity of the saturation scale. One can conclude that for p T > 2 Q s we can safely use the gluon structure function which has been measured for DIS at HERA. On the other hand, the thermal emission comes from the region p T < 2 Q s and, therefore, the existence of this phenomenon depends completely on the inclusive prodiction in CGC/saturation approach in this region.