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

TL;DR

This paper analyzes inclusive hadron production at high energies within the CGC/saturation framework, highlighting the role of saturation momentum and thermal radiation contributions, and comparing different confinement models for data description.

Contribution

It introduces a model incorporating thermal radiation in the CGC approach and compares it with alternative confinement models to explain experimental data.

Findings

Thermal radiation improves data fit in the CGC model.

Gluon jet production mainly occurs near the saturation momentum.

Alternative confinement models can describe data without thermal emission.

Abstract

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

Tables

Table 1. Table 1: Fitted parameters of the model RESH , which we use in our estimates.

$\overline{\gamma }$	$N_0$	$\lambda$	$x_0$	$B_{CGC}$ ($Ge{V}^{-2}$)
0.6599 $\pm$ 0.0003	0.3358$\pm$ 0.0004	0.2063 $\pm$ 0.0004	$0.00105\pm {1.1310}^{-5}$	5.5

Table 2. Table 2: R = d 2 ​ σ / d 2 ​ p T ​ ( thermal radiation ) / d 2 ​ σ / d 2 ​ p T ​ ( sum ) 𝑅 superscript 𝑑 2 𝜎 superscript 𝑑 2 subscript 𝑝 𝑇 thermal radiation superscript 𝑑 2 𝜎 superscript 𝑑 2 subscript 𝑝 𝑇 sum R\,\,=\,\,d^{2}\sigma/d^{2}p_{T}(\mbox{thermal radiation})\Big{/}d^{2}\sigma/d^{2}p_{T}(\mbox{sum}) versus the values of energies and the value of m eff subscript 𝑚 eff m_{\rm eff} .

W (TeV)	$m_{\mathrm{eff}}$ = 0.5  GeV)	$m_{\mathrm{eff}}$ = 0.12  GeV
13	70%	46%
7	70%	43%
2.76	59%	13%
0.9	53%	7%

Equations

$$\frac{d\sigma }{dy{d}^2{p}_T}=\underset{\text{mbox}thermalradiation}{\underbrace{A_{\mathrm{therm}}{e}^{-\frac{m_T}{T_{\mathrm{th}}}}}}+\underset{\text{mbox}hardemission}{\underbrace{A_{\mathrm{hard}}\frac{1}{{\left(1+\frac{m_T^2}{T_{\mathrm{h}}^{2}n}\right)}^n}}}$$

$$T_{\mathrm{th}}=0.098{\left(\sqrt{\frac{s}{s_0}}\right)}^{0.06}\mathrm{GeV};T_{\mathrm{h}}=0.409{\left(\sqrt{\frac{s}{s_0}}\right)}^{0.06}\mathrm{GeV};$$

$$T_{\mathrm{th}}=\mathrm{c}\frac{Q_s}{2\pi }$$

$$Q_s\left(x,b\right)=Q_0\left(b\right){\left(\frac{1}{x}\right)}^{\lambda /2}=Q_s\left(x\right)S\left(b\right)$$

$$T_{\mathrm{th}}\propto T_{\mathrm{h}}\propto Q_s\propto {\left(\sqrt{\frac{s}{s_0}}\right)}^{\lambda /2}\sim {\left(\sqrt{\frac{s}{s_0}}\right)}^{0.1-0.112},$$

$$\frac{d{\sigma }_G}{dy{d}^2{p}_T}$$

$$N_G\left(y_i=\ln (1/x_i);r,b\right)=2N\left(y_i=\ln (1/x_i);r,b\right)-N^2\left(y_i=\ln (1/x_i);r,b\right)$$

$$x_1=\frac{p_T}{W}{e}^y;x_2=\frac{p_T}{W}{e}^{-y};$$

$$\frac{d{\sigma }_G}{dy{d}^2{p}_T}=$$

$$\frac{2C_F}{{\alpha }_s(\pi {)}^4}\frac{1}{p_T^2}\int d^{2}r{e}^{i{\mathbf{p}}_T\cdot \mathbf{r}}\int d^{2}b{Q}_s^2\left(x_1,b\right)e^{-z\left(r,b,x_1\right)}\frac{d^2{N}_G\left(z(r,b,x_1))\right)}{d{z}^2}\int d^2{b}^{\prime }{Q}_s^2\left(x_1,b\right)e^{-z\left(r,b^{\prime },x_2\right)}\frac{d^2{N}_G\left(z(r,b^{\prime },x_2)\right))}{d{z}^2}$$

$$\displaystyle=\frac{2C_{F}}{\alpha_{s}(\pi)^{3}}\,\frac{Q^{2}_{s}\left(x\right)}{p^{2}_{T}}\,\int dz\,e^{-z\left(r,b=0,x\right)}\,J_{0}\left(\frac{p}{Q(x)}e^{z}\right)\,\,\,\Bigg{(}\int d^{2}b\,S\left(b\right)\frac{d^{2}N_{G}\left(z(r,b,x_{1}))\right)}{dz^{2}}\Bigg{)}^{2}$$

$$\equiv \frac{2C_F}{{\alpha }_s(\pi {)}^3}\frac{1}{{\tilde{p}}_T^2}\int dz{J}_0\left({\tilde{p}}_T{e}^z\right)I\left(z\right)=\frac{2C_F}{{\alpha }_s(\pi {)}^3}\frac{1}{{\tilde{p}}_T^2}\mathcal{I}\left(\tilde{p}\right)$$

$$\frac{d{\sigma }^{\pi }}{dy{d}^2{p}_T}={\int }_0^{1}d{x}_G\frac{d{\sigma }_G}{dy{d}^2{p}_T}\left(\frac{p_T}{x_{\pi }}\right)D_G^{\pi }\left(x_{\pi }\right)$$

$$D_G^{\pi }\left(x^{\pi }\right)=2.17z^{\alpha }(1-z{)}^{\beta }\left(20(1-z{)}^{\gamma 1}+1\right);$$

$$x=\frac{p_T}{W}=\frac{Q_s\Theta (Q_s-p_T)+p_T\Theta (p_T-Q_s)}{W}.$$

$$\frac{\partial N\left(Y;{\mathbf{x}}_{01},\mathbf{b}\right)}{\partial Y}$$

$$K\left({\mathbf{x}}_{01};{\mathbf{x}}_{02},{\mathbf{x}}_{12}\right)⟶\chi \left(\gamma \right)=2\psi (1)-\psi (\gamma ))-\psi (1-\gamma )$$

$$\displaystyle\chi\left(\gamma\right)\,\,=\,\,\left\{\begin{array}[]{l}\,\,\,\frac{1}{\gamma}\,\,\,\,\,\,\,\,\,\,\mbox{for}\,\,\,\tau\,=\,rQ_{s}\,<\,1\,\,\,\,\,\,\mbox{summing}\left(\ln\left(1/(r\,\Lambda_{\rm QCD})\right)\right)^{n};\\ \\ \,\,\,\frac{1}{1\,-\,\gamma}\,\,\,\,\,\mbox{for}\,\,\,\tau\,=\,rQ_{s}\,>\,1\,\,\,\,\,\mbox{summing}\left(\ln\left(rQ_{s}\right)\right)^{n};\end{array}\right.$$

$$\int d^2{x}_{02}K\left({\mathbf{x}}_{01};{\mathbf{x}}_{02},{\mathbf{x}}_{12}\right)\to \pi x_{01}^2{\int }_{r^2}^{\frac{1}{{\Lambda }_{QCD}^2}}\frac{d{x}_{02}^2}{x_{02}^4}$$

$$\frac{\partial^{2}n\left(Y;\boldsymbol{x}_{01},\boldsymbol{b}\right)}{\partial Y\,\partial\ln\left(1/(x^{2}_{01}\,\Lambda^{2}_{QCD})\right)}\,\,\,=\,\,\frac{\bar{\alpha}_{S}}{2}\,\Big{(}2n\left(Y;\boldsymbol{x}_{01},\boldsymbol{b}\right)\,\,-\,n^{2}\left(Y;\boldsymbol{x}_{01},\boldsymbol{b}\right)\Big{)}$$

$$\int K\left({\mathbf{x}}_{01};{\mathbf{x}}_{02},{\mathbf{x}}_{12}\right)d^2{x}_{02}\to \pi {\int }_{1/Q_s^2(Y,b)}^{x_{01}^2}\frac{d{x}_{02}^2}{x_{02}^2}+\pi {\int }_{1/Q_s^2(Y,b)}^{x_{01}^2}\frac{d|{\mathbf{x}}_{01}-{\mathbf{x}}_{02}{|}^2}{|{\mathbf{x}}_{01}-{\mathbf{x}}_{02}{|}^2}$$

$$\frac{{\partial }^2\tilde{N}\left(Y;{\mathbf{x}}_{01},\mathbf{b}\right)}{\partial Y\partial \ln r^2}={\bar{\alpha }}_S\left\{\left(1-\frac{\partial \tilde{N}\left(Y;{\mathbf{x}}_{01},\mathbf{b}\right)}{\partial \ln x_{01}^2}\right)\tilde{N}\left(Y;{\mathbf{x}}_{01},\mathbf{b}\right)\right\}$$

$$\frac{{\partial }^{2}n\left(Y;{\mathbf{x}}_{01},\mathbf{b}\right)}{\partial Y\partial \ln \left(1/(x_{01}^2{\Lambda }_{QCD}^2)\right)}={\bar{\alpha }}_{S}n\left(Y;{\mathbf{x}}_{01},\mathbf{b}\right)$$

$$N\left(Y;x_{01},b\right)=N_0\exp \left(\sqrt{-{\xi }_s\xi }+\xi \right)\stackrel{\tau \to 1;\zeta \to 0_{-}}{\to }{N}_0{e}^{\frac{1}{2}z}\exp \left(-\frac{z^2}{8{\xi }_s}\right)$$

$${\xi }_s=4{\bar{\alpha }}_S\left(Y-Y_{\min }\right);\xi =\ln \left(x_{01}^2{Q}_s^2\left(Y=Y_{\min };b\right)\right);z={\xi }_s+\xi ;$$

$$N\left(Y;z=0_{-}(\xi =-{\xi }_s),b\right)=N_0\left(b\right);\frac{\partial \ln N\left(Y;z=0_{-}(\xi =-{\xi }_s),b\right)}{\partial \zeta }=\frac{1}{2};$$

$$N\left(Y;r,b\right)=N_0{\left(r^2{Q}_s^2\left(Y,b\right)\right)}^{1-{\gamma }_{cr}}$$

$$\frac{d^{2}n\left(z\right)}{dz^{2}}\,\,\,=\,\,\frac{1}{8}\,\Big{(}2n\left(z\right)\,\,-\,n^{2}\left(z\right)\Big{)}$$

$$\frac{d\,p^{2}}{d\,n}\,=\,\frac{1}{4}\Big{(}2n\,\,-\,n^{2}\Big{)}$$

$$p^{2}\,\,=\,\,\frac{1}{4}\Big{(}n^{2}\,\,-\,\frac{n^{3}}{3}\Big{)}+{\rm C}_{1}$$

$$2\int^{n}_{n_{0}}\frac{dn^{\prime}}{n^{\prime}}\frac{1}{\sqrt{1-n/3}}\,\,=\,\,z\,\,\,=\,\,\,2\Bigg{(}\ln\frac{\left(1-\sqrt{1-\frac{n}{3}}\right)}{\left(1+\sqrt{1-\frac{n}{3}}\right)}\,\,-\,\,\ln\frac{\left(1-\sqrt{1-\frac{n_{0}}{3}}\right)}{\left(1+\sqrt{1-\frac{n_{0}}{3}}\right)}\Bigg{)}$$

$$n\left(z\right)\,\,=\,\,n_{0}\,e^{\frac{1}{2}\,z\,\,-\,\,\frac{n_{0}}{6}\Big{(}e^{\frac{1}{2}\,z}\,-\,1\Big{)}}$$

Full text

DGLAP and BFKL evolution, double parton distributions, Bose-Einstein correlations, shadowing corrections, non-linear evolution equation, CGC approach.

††preprint: TAUP - 3030/18

1. I Introduction

2. II Inclusive production in CGC/saturation approach: generalities

3. III Non-linear evolution for the leading twist BFKL kernel

4. III.1 Equations

5. III.2 Solutions

6. III.2.1 Perturbative QCD: linear equation

7. III.2.2 Solution in the region $\tau\,<\,1$

8. III.2.3 Solution in the region $\tau\,>\,1$

9. III.2.4 Inclusive production

10. IV Impact-parameter dependent CGC dipole model

11. V Comparison with the experimental data

12. VI Conclusions

13. VII Acknowledgements

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.KOLEB 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$GLR ; MUQI ; MV . 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.ALICE13 , which shows that the next-to-leading order QCD calculations with formation of the hadrons in accord with the fragmentation functions (NLOFIT ; NLOFIT1 ), is able to describe the transverse momentum spectra for the LHC range of energies. However, such a description is only successful for large values of $p_{T}\,>\,3\,GeV$NLOFIT or $p_{T}\,>\,5GeV$NLOFIT1 , while we assume that one can use these fragmentation functions in the region of small $p_{T}$ as well. In a sense, at present, this model is the best that we can propose to describe multi - hadron production.

It turns outFIT0 ; FIT1 ; FIT2 ; FIT3 ; BAKH ; FPV that the experimental data ALICE13 ; ALICE2 ; ATLAS1 ; ATLAS2 ; CMS1 ; CMS2 at high energy, can be describes as the sum of two terms:

[TABLE]

with

[TABLE]

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 ($\propto 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. Therefore, the appearance of the first term in Eq. (1) is not related to the equilibrium of the produced system of partons, and we used the word ‘thermal’ just to characterize the form of $p_{T}$ dependence of this term. The origin of thermal radiation in the framework of the CGC approach was clarified a decade ago KHTU ; DUHA ; CKS ; KLT and, recently, the new idea that the quantum entanglement is at the origin of the parton densities has been added to these argumentsKHLE . The resulting picture is presented nicely in Ref.BAKH , 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 $\Delta t\propto 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 stateKHLE 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 KHTU ; KLT that

[TABLE]

with the semi-classical estimates KLT 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 formGLR ; MUQI ; MV :

[TABLE]

The value of $\lambda$ can be calculated theoretically and measured experimentally. The leading order QCD evaluation leads to $\lambda=4.9\bar{\alpha}_{S}$,where $\bar{\alpha}_{S}$ denotes the running QCD coupling. Plugging in the reasonable estimate for $\bar{\alpha}_{S}\left(Q_{s}\right)\approx 0.2$, one can see that $\lambda$ turns out to be large, about 0.8-1. The phenomenological description of the hard processes both for nucleus interactionsDKLN and DIS( see Ref.RESH and references therein), give the value of $\lambda=0.2-0.24$. We therefore see, that in the CGC approach we expect that

[TABLE]

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 that this dependence should be proportional to the saturation scale. 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 attemptLERE to compare the CGC prediction with the experiment at $W\,=\,7\,TeV$, 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.KTINC for the proof):

[TABLE]

where $N_{G}\left(y_{1}=\ln(1/x_{1});r,b\right)$ can be found from the amplitude of the dipole-proton scattering $N\left(y_{i}=\ln(1/x_{i});r_{;}b\right)$ :

[TABLE]

where $r$ denotes the size of the dipole, $b$ it’s impact parameter and

[TABLE]

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 $\bar{\alpha}_{S}\,=\,\alpha_{S}N_{c}/\pi$ with the number of colours equals $N_{c}$. $\alpha_{S}$ denotes the running QCD coupling, $\nabla^{2}_{T}$ the Laplace operator with respect to $r$, it is equal to $\nabla^{2}_{T}\,=\,\frac{1}{r}\frac{d}{dr}\left(r\frac{d}{dr}\right)$.

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: $\tau\,=\,r\,Q_{s}\left(W,b\right)$. Introducing a new variable $z\left(r,b,x\right)\,\,=\,\,\ln\left(\tau^{2}\right)$, we can re-write Eq. (5) at $y=0$ in the form

[TABLE]

In Eq. (II) we have taken $x=x_{1}=x_{2}$ since $y=0$, and introduce a new variable $\tilde{p}_{T}\,=\,p_{T}/Q_{s}\left(x\right)$.

To calculate the hadron distributions, we need to take into account the decay of the jet into hadrons. The formula has the form

[TABLE]

We take the fragmentation function $D^{\pi}_{G}$ from Ref.FRFU which has the form

[TABLE]

with $\alpha=0.899$, $\beta=1.57$ and $\gamma=4.91$.

We note that Eq. (II) 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.This divergency also indicates that we need to reformulate our assumptions about confinement, since using the fragmentation functions does not suppress the divergency at low $p_{T}$. We believe that the reason for this divergency, is the fact that we neglected the mass of the jet of hadrons that stem from the decay of the gluon. The simple estimates KHLEKLN give for a gluon with the value of the transverse momentum $p_{T}$, the mass of the jet $m^{2}_{\rm jet}\,=\,2p_{T}\,m_{\rm eff}$, where $m_{\rm eff}=\sqrt{m^{2}+k_{T}^{2}+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}\approx\,k_{T}$ is the longitudinal momentum of this hadron. Since most pions stem from the decay of $\rho$-resonances we expect that $m_{\rm eff}\approx m_{\rho}$.

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}_{\rm jet}\,=\,2Q_{s}\,m_{\rm eff}$. However, for rare gluons with $p_{T}\gg Q_{s}$ we still have $m^{2}_{\rm jet}\,=\,2p_{T}\,m_{\rm eff}$ . For numerical estimates we use $m^{2}_{\rm jet}\,=\,2\left(Q_{s}\Theta(Q_{s}-p_{T})+p_{T}\Theta(p_{T}-Q_{s})\right)m_{\rm eff}$ which has these two limits. $\Theta(x)$ denote the step function. Using the same idea we replace Eq. (7) by

[TABLE]

Hence, we finally deal with the following model for the confinement ( hadronization): the decay of the gluon jet with the effective mass that we have discussed above, and with fragmentation functions of Eq. (10) at all values of the transverse momenta.

III Non-linear evolution for the leading twist BFKL kernel

III.1 Equations

The dipole scattering amplitude is the solution of the Balitsky-Kovchegov BK non-linear equation, which has the following form:

[TABLE]

$N\left(Y;\boldsymbol{x}_{01},\boldsymbol{b}\right)$ denotes the dipole scattering amplitude. $\boldsymbol{x}_{01}=\boldsymbol{x}_{1}-\boldsymbol{x}_{0}\equiv\boldsymbol{r}$ the size of the dipole. The kernel $K\left(\boldsymbol{x}_{01};\boldsymbol{x}_{02},\boldsymbol{x}_{12}\right)$ describes the decay of the dipole with size $x_{01}$ into two dipoles of size: $\boldsymbol{x}_{02}$ and $\boldsymbol{x}_{12}=\boldsymbol{x}_{01}-\boldsymbol{x}_{02}$. The linear part of this equation reduces to the BFKL equationBFKL with the eigenfunctions $\left(r^{2}\right)^{\gamma}$ which corresponds to the value of the kernel:

[TABLE]

where $\psi(z)=d\ln\Gamma(z)/dz$ is the Euler $\psi$-function (see Ref.RY 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 KOWI 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

[TABLE]

instead of the full expression of Eq. (13).

As indicated in Eq. (17) we have two types of logs: \Big{(}\bar{\alpha}_{S}\ln\left(r\,\Lambda_{QCD}\right)\Big{)}^{n} in the perturbative QCD kinematic region where $r\,Q_{s}\left(Y,b\right)\,\,\equiv\,\,\tau\,\,\ll\,\,1$; and \Big{(}\bar{\alpha}_{S}\ln\left(r\,Q_{s}\left(Y,b\right)\right)\Big{)}^{n} inside the saturation domain ($\tau\,\gg\,1$), where $Q_{s}\left(Y,b\right)$ denotes the saturation scale. To sum these logs it is necessary to modify the BFKL kernel in different ways in the two kinematic regions, as shown in Eq. (17).

$\bullet\,\,{\it\tau\,\,=\,\,r\,Q_{s}\ll\,1}$

For the perturbative QCD region of $\tau\,\ll\,1$, the logs originate from $x^{2}_{02}\sim x^{2}_{12}\,\ll\,x^{2}_{01}\equiv r$ resulting in the following form of the kernel $\displaystyle{K\left(\boldsymbol{x}_{01};\boldsymbol{x}_{02},\boldsymbol{x}_{12}\right)}$ LETU

[TABLE]

The non-linear BK equation in this region can be written as

[TABLE]

for $n\left(Y;\boldsymbol{x}_{01},\boldsymbol{b}\right)\,=\,N\left(Y;\boldsymbol{x}_{01},\boldsymbol{b}\right)/x^{2}_{01}$ .

$\bullet\,\,{\it\tau\,\,=\,\,r\,Q_{s}\gg\,1}$

Inside the saturation region where $\tau\,\,>\,\,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

[TABLE]

Inside the saturation region the BK equation takes the form

[TABLE]

where $\widetilde{N}\left(Y;\boldsymbol{x}_{01},\boldsymbol{b}\right)\,\,=\,\,\int^{x^{2}_{01}}dx^{2}_{02}\,N\left(Y;\boldsymbol{x}_{02},\boldsymbol{b}\right)/x^{2}_{01}$ .

The advantage of the simplified kernel of Eq. (17) is that, in the Double Log Approximation (DLA) for $\tau<1$, it provides a matching with the DGLAP evolution equationDGLAP .

III.2 Solutions

III.2.1 Perturbative QCD: linear equation

First we discuss the solution to Eq. (19) in the perturbative QCD region, where one can neglect the contributions of the shadowing corrections. The equation has the form:

[TABLE]

One can recognize that Eq. (22) is the DGLAP equationDGLAP in the double log approximation (DLA). It has the following solution (see for example Ref.KOLEB )

[TABLE]

where we use the following notation:

[TABLE]

The solution of Eq. (23) provides the boundary condition for the solution inside the saturation region:

[TABLE]

As was expectedGS3 , in the vicinity of the saturation scale ( $\zeta\ll 8\xi_{s}$), the amplitude exhibits geometric scaling behavior, being a function of only one variable $\zeta$ and it has the following formMUT :

[TABLE]

where $\gamma_{cr}$ denotes the critical anomalous dimension which in the DLA is equal to $\frac{1}{2}$.

III.2.2 Solution in the region $\tau\,<\,1$

The solution of Eq. (26) 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. (19) can be simplified assuming the geometric scaling behaviour of the amplitude in the vicinity of the saturation scaleGS3 . It has the form:

[TABLE]

We can obtain the solution to Eq. (27) introducing $dn/dz=p\left(n\right)$ . For this function we obtain the equation:

[TABLE]

which has the solution

[TABLE]

Since $p\propto 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:

[TABLE]

The first correction to the scattering amplitude of order $n^{2}$ has the following form:

[TABLE]

III.2.3 Solution in the region $\tau\,>\,1$

In the saturation region it has been shown GS1 that the scattering amplitude manifests geometric scaling behavior, which is supported by the experimental dataGS2 . Therefore, Eq. (21) can be re-written in the form:

[TABLE]

Introducing

[TABLE]

we re-write Eq. (32) in the form:

[TABLE]

Integrating we obtain

[TABLE]

which can be resolved as

[TABLE]

Taking the derivatives on both sides of Eq. (36) we reduce this equation to the form:

[TABLE]

Therefore, we can find $\phi$ as the solution to the following equation:

[TABLE]

The value of $C$ has to be found from matching with the region $\tau\,<\,1$. For small $\phi_{0}$   $C=0$. Indeed, in this case the solution at small $\phi$ has the following form:

[TABLE]

which coincides with Eq. (26) for the region $\tau\,<\,1$ at small $\phi_{0}\,=\,N_{0}\,\ll\,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. (31):

[TABLE]

III.2.4 Inclusive production

In our approach we can estimate the value of function ${\cal I}\left(\tilde{p}\right)$ in Eq. (II). Taking $N_{0}=0.05$, and considering this as a small value, we can use Eq. (38) 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^{\prime}_{z}\,=\,dN/dz$ and $N^{\prime\prime}_{\!\!zz}\,=\,d^{2}N/dz^{2}$ become very small at large $z$. The integrand in Eq. (II) for the inclusive production is equal to $e^{-z}\,N^{\prime\prime}_{\!\!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. (26). It turns out that the behaviour for $z>0$ can be approximated by the simple formula for the scattering amplitude (see Fig. 1-b):

[TABLE]

Hence, one can see that the situation for inclusive production looks quite different to that for the description of the DIS structure function. $\nabla^{2}_{T}$ in Eq. (II) 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 ${\cal I}\left(\tilde{p}\right)$ at $p_{T}=0$. In ${\cal I}\left(\tilde{p}_{T}=0\right)$ 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 ${\cal I}\left(\tilde{p}_{T}=0\right)$ 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\to-\infty$, so, we need to generalize the behaviour of $\phi$ in the perturbative QCD region (see Eq. (23)) by replacing

[TABLE]

with

[TABLE]

Eq. (43) is Eq. (23) re-written in the explicit form.

Taking $\bar{\alpha}_{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. (II):

[TABLE]

The solution $\phi=\phi_{0}e^{\frac{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 $\phi^{\prime}$ in Eq. (38) keeping $C=-N^{3}_{0}/2$ in the dominator(see Eq. (III.2.3)). The solution for $\phi(z)$ at small $z$ has the form

[TABLE]

One can see that $\phi_{2}(z)$ is close to $\phi(z)=N_{0}\,e^{\frac{1}{2}\,z}$, even for $N_{0}=0.3-0.4$.

We wish to emphasize that the contribution to ${\cal I}\left(\tilde{p}_{T}\right)$ at sufficiently short distances does not lead to the suppression of this function at $\tilde{p}_{T}\,=\,0$. Therefore, the cross section is still divergent at $p_{T}\to 0$, which results in the large production of soft gluons in the framework of CGC/saturation approach, and is suppressed by the process of hadronization.

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 formMUT :

[TABLE]

where $\bar{\gamma}\,=\,1-\gamma_{cr}$ and $\gamma_{cr}\,=\,0.37$ in the leading order is the solution to the equationKOLEB :

[TABLE]

$\chi\left(\gamma\right)$ is given by Eq. (13).

The advantage of Eq. (49) 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:

[TABLE]

In the LO BFKL KOLEB $\lambda\,\,=\,\,\bar{\alpha}_{S}\,\frac{\chi\left(\gamma_{cr}\right)}{\bar{\gamma}}$. Parameters $N_{0}$ and $Q_{0}$, as well as function $S\left(b\right)$ 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 modelsRESH ; CLP *** The DIS data has also been described by the saturation model of Ref.RSKV , as this model does not produce the theoretically correct behaviour deep in the saturation region, we do not consider it in this paper. on the market that describe the final set of the HERA experimental data on deep inelastic structure functionsHERADATA . They have different forms for $S\left(b\right)$:

[TABLE]

The ansatz of Eq. (53) is preferable, since it leads to $S\left(b\right)\xrightarrow{b\gg 1/m}\,\exp\left(-mb\right)$, which is in accord with the Froissart theoremFROI . However, we choose Eq. (52) which allow us to do several integrations analytically. Using Eq. (57) $\nabla^{2}_{T}N$ takes the following form after integrating over $b$:

[TABLE]

with $\phi_{0}\,e^{-2\phi_{0}}\,=\,N_{0}\,\left(1-N_{0}\right)$.

Before discussing the details about our model, we would like to outline, which features of Eq. (57) stems from the theory, and which from phenomenological assumptions. The expression for $\tau\,\leq\,1$, as we have mentioned (see Eq. (49)) 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}\tau^{2\bar{\gamma}}$ as a phenomenological expression that describes the DIS dataRESH . However, we check how large the contribution to the inclusive production from the corrections of Eq. (30) and Eq. (31) 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\left(b\right)$.

For $\tau\,\geq\,1$ we have discussed the form of Eq. (57) in the previous section, and have given strong arguments for such an expression. The $b$ integration is performed with the phenomenological $S\left(b\right)$.

In our estimates we use the values of the parameters from Ref.RESH (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†††Note that we introduce the extra factor $\frac{1}{2}$ in the definition of the saturation scale since we use $\tau=r\,Q_{s}$, while in Ref.RESH $\tau$ is defined as $\tau=rQ_{s}/2$.

[TABLE]

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.

In Fig. 3 we plot functions $\int d^{2}b\,\nabla^{2}_{T}N\left(z(r,y,b)\right)$ and $I\left(\tau\right)$ of Eq. (II). One can see that the both functions provide the main contribution from the region $\tau\,<\,1$. The largest contribution to the function ${\cal I}\left(\tilde{p}=0\right)$ in the region of integration $\tau\,>\,1$, is 13.5% . For large $p_{T}$ we expect that the contribution of the region $\tau\,>\,1$ to the function ${\cal I}\left(\tilde{p}_{T}\right)$ is small. In this region we have an analytical expression for this function:

[TABLE]

In Fig. 4 we plot the numerical result for ${\cal I}\left(\tilde{p}_{T}\right)$ with functions given by Eq. (57) and the analytical expression of Eq. (59). Note that for $\tilde{p}_{T}\,\geq\,2$ (or $p_{T}\,\geq\,2\,Q_{s}(x)$ ) both functions coincide.

The dashed line in Fig. 4 illustrates the result of the analytical integration over $\tau$, assuming that the region for $\tau>1$ does not contribute. The result of this integration has the form:

[TABLE]

From Fig. 4 we can conclude that the saturation region with $\tau\,\geq\,1$, only provides around 10-15% of the contribution for $\tilde{p}_{T}\,\leq\,1$, but the sharp cutoff changes the behaviour for large values of $\tilde{p}_{T}$.

We know that in the vicinity of the saturation scale the scattering amplitude in the momentum representation has the following behaviour:

[TABLE]

Therefore, from Eq. (59) we can determine the value of constant in Eq. (61) from the value of $N_{0}$.

For large $p_{T}\,\gg\,Q_{s}$ we cannot use neither Eq. (25) nor Eq. (61), since we have to take into account the violation of the geometric scaling behaviour in perturbative QCD (see Eq. (23)). Eq. (23) indicates that part of this violation can be taken into account by replacing in Eq. (25):

[TABLE]

In the DLA Eq. (62) gives Eq. (43). Eq. (62) is used in Ref.RESH for fitting the HERA experimental data.

We make such a replacement directly in the momentum representation, since $r\propto 1/p_{T}$. However, we need to find the coefficient in front of $p_{T}$ and, perhaps, an additional constant. We calculate the average $\tau$ using the expression:

[TABLE]

The results of these estimates are shown in Fig. 5. One can see at at $\tilde{p}_{T}\to 0$ $\langle\tau\rangle\,=\,0.478\approx\frac{1}{2}$ while at large $\tilde{p}_{T}$ it is proportional to $1/(4\tilde{p}_{T})$.

For $\tilde{p}_{T}=0$, or more generally for $p_{T}\ll Q_{s}$, the typical distances turns out to be $r=1/(2Q_{s}(x))$, and for large $\tilde{p}_{T}$ they are of the order of $1/(4p_{T})$. Hence, we suggest to use in Eq. (62) the calculated $\langle\tau\rangle\left(\tilde{p}_{T}\right)$ for $\tilde{p}\leq 4$ and $1/(4\tilde{p}_{T}$) for $\tau\geq 4$.

In Fig. 6 we show the behaviour of the $\gamma_{\rm 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. (59) for example) at large $p_{T}$ the inclusive cross section is proportional to $1/p_{T}^{4\gamma_{\rm eff}}$. For $\gamma_{\rm eff}=\bar{\gamma}$ 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 $\gamma_{\rm eff}$ of Eq. (62) is shown in Fig. 6. We see that even this behaviour of $\gamma_{\rm eff}$ leads only to $1/p^{4-5}_{T}$ at $p_{T}\approx$ 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 NLOFIT . 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}\leq$ 7 GeV, and on the value of the contribution of the thermal radiation, which is negligibly small at $p_{T}\sim 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. (II). As we have discussed in function ${\cal I}\left(\frac{p_{T}}{Q_{s}}\right)$ the main contribution comes from $p_{T}\sim Q_{s}$ (see Fig. 4). However, the factor $1/p^{2}_{T}$ in front in Eq. (II) stems from the gluon propagator BFKL , 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 propagatorKHLEKLN :

[TABLE]

Therefore, we calculate gluon production using Eq. (II) in which we use Eq. (64) to replace the factor $1/p^{2}_{T}$.

For calculating $I\left(\tau,W\right)$ in Eq. (II) we use Eq. (57) for $p_{T}\,\leq\,2\,Q_{s}(x)$ and Eq. (59) for larger values of $p_{T}$. It should be noted that the original model of Ref.RESH has a different assumption regarding the behaviour of the amplitude in the saturation region. However, $\nabla^{2}N$ in this function is not continuous at $\tau=1$ and leads to a very small contribution for $\tau\geq 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 $\nabla^{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.RESH reproduces only the the theoretical behaviour of Ref.LETU at large $\tau$ or, in other words, in the region which does not contribute to $\nabla^{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 $\gamma_{\rm 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_{\rm eff}$.

Fig. 9 illustrates how well we describe the data for the transverse momentum distribution at the LHC.

$T_{\rm th}\,\propto\,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.KLT .

The rate of thermal radiation is shown in Table 2, in which $R=\int d^{2}p_{T}\,d^{2}\sigma^{\rm charged}_{\mbox{therm. rad.}}/d^{2}p_{T}/\int d^{2}p_{T}\,d^{2}\sigma^{\rm charged}_{\mbox{sum.}}/d^{2}p_{T}$. 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_{\rm eff}$. We believe that most of the pions are produced from $\rho$ resonances and we consider $m_{\rm eff}=$ 0.5  GeV as the most reliable estimate. In Fig. 10 we present the calculation at W = 7 TeV with $m_{\rm eff}=0$ and $m_{\rm 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 $\sigma_{in}\,\,=\,\,\sigma_{tot}-\sigma_{el}-\sigma_{diff}$ from Ref.GLP . We recall that the simple formula for $m_{\rm eff}=\sqrt{\mu^{2}+k^{2}_{T}+k^{2}_{L}}-k_{L}$ leads to $m_{\rm eff}=0.5\,GeV$ if $\mu$ is equal to the mass of $\rho$-resonance since the value of $k_{T}=k_{L}=0.45\,GeV$ (see Ref.ATLASJET for the measurement and Ref.KLRX ) and reference therein for theoretical discussions). For the minimal mass of $\mu=m_{\pi}=0.14\,GeV$ we obtain $m_{\rm 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 $\Lambda_{\rm 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}=\Lambda_{\rm 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_{\rm th}$, acquires a mass $m_{g}\,\,\propto\,\,T_{\rm th}$BRPI 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}\approx T_{\rm th}\approx 0.12-0.14\,GeV$. In this picture we will not see any thermal emission in the spectrum of hadrons. Note, that $p_{T}\sim 0.12-0.14\,GeV$ corresponds to the $m_{\rm eff}\approx 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_{\rm th}$. One can see that we are able to successfully describe the data without the thermal radiation term. Such a description should only be considered with a grain of salt, since $m=gT$ with small $g$ in Ref.BRPI , and realistic estimates will overshoot the data.

In all the examples above, we change our model for the hadronization, adding the mass of the gluon jet. To illustrate our claim, that the existence of the thermal radiation crucially depends on the model of confinement, we plot in Fig. 12 the inclusive spectra for two different models: in Fig. 12-a we assume that the gluons with $p_{T}<Q_{s}$ do not take part in the hadronization, and in Fig. 12-b only gluons with $p_{T}>\Lambda_{QCD}$ produce the jet of hadrons. Fig. 12-a shows that we need the thermal radiation term with the contribution of 56% to describe the data, while in Fig. 12-b the data do not require the thermal emission.

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 hadrons with the given fragmentation functions. The temperature of this emission turns out to be equal to $2.3/(2\,\pi)\,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.KLT .

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. 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.VAZW ; KHALE 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}_{\rm eff}=2Q_{s}\mu_{\rm soft}$ where $\mu_{\rm soft}$ is the soft scale, and with the fragmentation functions of Eq. (10) at all values of the transverse momenta.

We suggest to take into account the behaviour of $\nabla^{2}N$ in the saturation region in the form:

[TABLE]

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 for higher order corrections. We wish also to mention, that we have calculated the inclusive production taking $\bar{\alpha}_{S}=0.25$. This value is less that $\bar{\alpha}_{S}\left(Q_{s}\right)\,\,=\,\,0.3$ which appears more natural in Eq. (II). For $\bar{\alpha}_{S}=\bar{\alpha}_{S}\left(Q_{s}\right)$, 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_{\rm 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}\,\leq\Lambda_{\rm QCD}$ or/and that generates the gluon mass $m=T_{\rm th}$, is able to describe the experimental data without the thermal radiation.

Hence, we state that the existence of the thermal term in the $p_{T}$ spectrum of produced hadrons, depends crucially on the model for hadronization.

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 $\int d^{2}bN\left(Y,r,b\right)$. 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. 13 we show Eq. (59) 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.

VII Acknowledgements

We thank our colleagues at Tel Aviv university and UTFSM for encouraging discussions. Our special thanks go to Keith Baker and Dmitry Kharzeev for fruitful discussions on the subject which prompted the appearance of this paper.

This research was supported by Proyecto Basal FB 0821(Chile), Fondecyt (Chile) grant 1180118 and by CONICYT grant PIA ACT1406.