The role of nuclear gluon distribution on particle production in heavy ion collisions

The transverse momentum spectra of hadrons is calculated from the unintegrated gluon distribution (UGD) within the $k_T$-factorization framework at small $x$. Starting from $pp$ collisions, the modification caused by the nuclear medium is incorporated in the UGD at high energies, which is related to the nuclear shadowing phenomenon. Moreover, we consider that particle production from minijet decaying is not enough to explain the $p_T$ spectra in $AA$ collisions due to collective phenomena that take place after the hard collision. The Boltzmann-Gibbs Blast Wave (BGBW) distribution is utilized in order to evaluate the distribution of particle production in equilibrium. Data from ALICE collaboration for $PbPb$ collisions at $\sqrt{s}=2.76$ TeV are analyzed and the nuclear modification factor for pion production is computed.


I. INTRODUCTION
The modification of nuclear structure functions at small Bjorken x, compared to those for free nucleon observed in DIS [1] (i.e, EMC effect, anti-shadowing and shadowing), can be attributed to distinct phenomena [2][3][4]. In particular, the shadowing seen in nuclear DIS [5] for x 0.01 is characterized by depletion of F A 2 with respect to F N 2 . This effect has been explored within the Color Glass Condensate (CGC) saturation approach [6], which in turn predicts a saturation scale, Q s (x). Such quantity establishes the region where the increasing of the unintegrated gluon distribution (UGD) on x is tamed. It is expected that for heavy ion collisions the saturation scale is enhanced compared to the nucleon case and it is deeply connected to shadowing corrections. Hence, the nuclear saturation scale characterises the dense system. In pA collisions, the depletion of cross section related to pp mode is also evident in the region of small p T and may be well described in the context of saturation formalism [7][8][9]. However, in AA collisions the particle production presents more complex behavior which can not be explained solely by nuclear effects in the gluon distribution at initial state. While the collinear factorization mechanism is well established at large Q 2 , the saturation formalism/CGC makes use of k T -factorization in the small x domain in order to describe minijet production of gluons where the UGD, φ(x, k T ), depends on the transverse momentum k T . This distribution is related to the dipole cross section, σ dp (x, r), being the latter directly extracted from DIS at small-x. This formalism has been employed to describe the p T spectra of produced hadrons at RHIC and LHC [10][11][12][13][14][15][16][17].
In our previous work [18], we have computed the p T spectra of produced hadrons in pp collisions over a large interval of the scaling variable, τ = p 2 T /Q 2 s ( √ s). In addition, a phenomenological parametrization for an UGD * lucas.moriggi@ufrgs.br † guilherme.peccini@ufrgs.br ‡ magnus@if.ufrgs.br was proposed, which has a power-like behavior at high k T . As a result, it was shown that scaling is a good approximation for describing the spectra at different collision energies, √ s. Regarding the nuclear case, the gluon distribution should be properly modified in order to include cold matter effects. In the context of k T -factorization approach or related formalism, distinct ways of obtaining the nuclear UGD were proposed [19][20][21][22][23][24][25][26][27]. In the present work, the nuclear effects are incorporated in the nuclear UGD by means of the Glauber multiple scattering theory as performed in [13,14]. In pp collisions, particle production can be obtained within the k T -factorization assuming the local hadron-partonduality (LHPD). The final hadron spectra are directly related to those ones from produced gluons at initial state in a good approximation [18,28]. As already pointed out, such scenario might not be appropriated for collisions of heavy ions where there are collective effects that modify the p T spectra of produced hadrons with respect to the initial state. However, it has been argued in [29][30][31] that p T spectra can be well described by making a time separation in the relaxation time approximation (RTA) formalism of Boltzmann transport equation [32] among produced hadrons at the initial hard collision (which is parametrized by Tsallis [33] or Hagedorn distributions [34]) and produced hadrons by the system in thermal equilibrium. The collective radial flux plays an important role on the distribution form. Moreover, models with two components are successful over a wide amount of data concerning particle production at high collision energies [35][36][37]. The hadronic spectra is decomposed into two parts, being the first one related to Boltzmann statistics and the second one based on power law behavior that captures the aspects of perturbative QCD predictions.
In this work, we propose a nuclear unintegrated gluon distribution that embeds the shadowing verified in DIS. The cross section for minijet production, which is driven by gluons within the k T -factorization framework, is obtained. Additionally, the effects caused by the medium at final state are incorporated by the formalism of Boltzmann equation in the relaxation time approximation. The produced particles at the initial hard collision are arXiv:2012.05388v1 [hep-ph] 10 Dec 2020 calculated by using the k T -factorization, whereas the p T spectra due to hydrodynamics expansion is given by BGBW distribution. The distribution parameters are determined from data of pion production at different centralities as measured by ALICE for √ s = 2.76 TeV [38] and the nuclear modification factor, R AA , is predicted. Similar analysis has been performed in [31] without taking into account nuclear shadowing. In that sense, the aim here is to understand the impact caused by the modification of gluon distribution at small x on the observed nuclear modification factor. This paper is organized as follows. In the next section, we present the details on the determination of the nuclear UGD, φ A (x, k T ), from the free nucleon distribution, φ p (x, k T ). This is achieved using the multiple scattering formalism as well as the predictions of the spectra of produced hadrons using the k T -factorization with these modifications. We also describe the hadron production at final state from the BGBW distribution and constrain the relevant parameters. In Section III, predictions are compared against data for pion production and nuclear modification factor. An analysis on the interpretation of the obtained parameters is performed. Finally, in Section IV we outline the main results and present conclusions.

II. THEORETICAL FRAMEWORK AND MAIN PREDICTIONS
The nuclear UGD may be obtained from the nucleon distribution by using the Glauber-Mueller [39,40] approach for multiple scattering. It has been carried out, for instance, in Ref. [20]. In this case, the dipole scattering matrix in configuration space, r, can be determined from the cross section for dipole scattering off a proton, where T A (b) is the thickness function which depends on the impact parameter b. In the present work, a Woods-Saxon-like parametrization for the nuclear density [41] with normalization d 2 bT A (b) = A has been applied for a lead nucleus. The related nuclear UGD is given by being H 0 {f (r)} = rdrJ 0 (k T r)f (r) the Hankel transform of order zero. For the proton case, a homogeneous target with radius R p is considered so that the dependence on impact parameter is factorized as S dp (x, r, b) = S dp (x, r)Θ(R p − b). In the limit of large dipoles, S dp (x, r) → 0, and the cross section reaches a maximum, σ 0 = 2πR 2 p . Within the parton saturation formalism, the gluon distribution should have a maximum around k T = Q s (x). One of the features of this formalism is the presence of geometric scaling in the observables, which become dependent on the ratio Q 2 /Q s (x) rather than of Q 2 and x separately.
It has been proposed in [18] a gluon distribution based on geometric scaling of high p T spectra of produced hadrons in pp collisions, where the scaling variable is defined as τ = k 2 T /Q 2 s (x) and the parameter δn establishes the power-like behavior of the spectra of produced gluons at high momentum τ 1. The cross section for dipole scattering in coordinate space, r, may be written as in which τ r = rQ s (x) is the scaling variable in the position space and ξ = 1+δn. Therefore, by placing σ dp (x, r) in Eq. (1) the nuclear gluon distribution is obtained directly from Eq. (2). The spectra of produced gluons in the initial hard collision, given an impact parameter b, can be calculated in the k T -factorization formalism [42], Above, p T is the transverse momentum of the produced gluon and x A and x B are the gluon momentum functions in the nucleus A and B, respectively. They are expressed in terms of the rapidity y in the following way: In the LHPD approximation, the spectra of produced hadrons is directly related to the minijet originated in gluons. In that case, we consider the hadron being produced with momentum p T h = z p T , The parameters K and z are the same as those obtained for the spectra pp → π 0 + X. Also, it is important to notice that Eq. (5) diverges for p T → 0, so that one needs to apply a cut associated with the jet mass p 2 The suppression of hadron production in nuclear collisions due to the nuclear modifications on gluon distribution may be quantified through the following ratio: It can be seen that for small dipoles, i.e., r → 0 (or, equivalently, for high k T ), it is possible to expand Eq. (1), which leads to S dA ∼ T A (b)σ 0 S dp (x, r) and (b) → 1 for any value of b. This scenario is strictly valid for the case of particle production from minijet yield that is originated in the initial hard interaction. On the other hand, in heavy ion collisions the initial hard scattering is followed by the formation of Quark-Gluon Plasma (QGP) and the evolution of the system until the freeze-out in the hadronic phase has an important effect on final spectra. In [29] it was proposed that the evolution of particle distribution due to their interaction with the medium is set up by the Boltzmann transport equation within the relaxation time approximation (RTA), being t r the relaxation time and t f the time of freeze-out. Thereby, the hard initial distribution (t = 0) evolves until the final distribution f f in at t = t f . The equilibrium distribution, f eq , is characterized by the equilibrium temperature, T , and the relaxation time, t r , being the latter responsible for determining the amount of time until the system reaches equilibrium. Following [31], we have considered that particle distribution in equilibrium can be evaluated by the Boltzmann-Gibbs-Blast-Wave (BGBW) model [43], being I 0 and K 1 the Bessel functions of first and second kinds of order zero and one, respectively. The quantity m T is the transverse mass, m T = p 2 T h + m 2 h , and the velocity profile ρ is given by where β s is the maximum velocity expansion of the surface with average transverse velocity β = 2 2+m β s . The Tsallis distribution has been used in Refs. [29-31, 44, 45] in order to constrain the initial distribution without taking into account nuclear shadowing. That distribution can be deeply understood in the context of the fractal structures present in QCD or in general Yang-Mills theories [46,47]. It implies the need of Tsallis statistics (TS) to describe the thermodynamics aspects of the fields and the entropic index of the TS can be obtained in terms of the field fundamental parameters (see Ref. [48] for a recent review). In our approach, we have shown that one can describe pion and charged hadrons distributions through the k T -factorization framework. By utilizing the UGD of Eq. (3), such formalism produces a distribution with similar features of Tsallis distribution, which for τ ≥ 1 can be approximated (neglecting the nuclear effects) by The nuclear modification is introduced into the nuclear UGD and the hard initial distribution, f in , is taken from Eqs. (5) and (7). The saturation scale and the power index, δn, were parametrized in the following way [18]: in which the parameters a, b and x 0 were obtained by fitting the HERA data (see Ref. [18] for details). Given these considerations, the hadron production in nuclear collisions is expressed as the following sum: The distribution of produced particles in thermal equilibrium is given by Eq. (10) and the nuclear modification factor for each centrality class reads as where T AB is the mean value of nuclear overlap for a given centrality. It is relevant to stress out that this definition is the same of that one utilized experimentally for the determination of R AA in each centrality class. Furthermore, both the pp and AA spectra are calculated from the same model described before [18]. Albeit the parametrization for pp cross section has been made for the sum of charged hadrons and neutral pions, in nuclear collisions we have restricted the analysis for pion spectra. The reason is that proton production has strong influence on the region of middle p T spectra in AA collisions [49,50]. Such phenomenon is known as baryon anomaly and other mechanisms are needed in order to explain it.

III. RESULTS AND DISCUSSION
Our analysis on p T spectra is limited by the kinematic window determined by geometric scaling compatible with that observed at HERA [51][52][53][54]. For √ s = 2.76 TeV, one should have p T h 10 GeV. At first, the results of the nuclear modification factor are presented including only shadowing effect for different values of b in Eq. (8). In Fig. 1 it can be seen that R shadow P bP b considerably enhances until a maximum point around p T h ∼ 2 GeV, which is the well known Cronin peak [55]. This is a result of multiple scattering and the ratio further decreases until the limit R AA = 1 at high p T h . The position of this peak is set by the saturation scale, Q s (x), and by the mean value of the gluon momentum fraction carried by the hadron, since the pp cross section only depends on the scaling variable τ h = τ z 2. 33 . In such case we use z = 0.345 which is the fitted value for π 0 spectra in pp collisions for distinct values of √ s. Higher values of z should lead to the shift of this peak towards higher p T . In more central collisions, the shape of R AA has little dependence on the impact parameter, whereas for more peripheral collisions the nuclear effects are weaker. Other models based on geometric scaling were proposed in order to determine the nuclear shadowing within the dipole approach (as done in Ref. [10], for instance). However, we did not get good results when utilizing this picture because the resulting modification factor grows rapidly with p T . This issue has already been discussed in Ref. [56]. Furthermore, it is shown in [7] that multiple scattering formalism produces good results for the nuclear modification factor in pA collisions.
The determination of hadron final distribution, Eq. (15), is carried out by fitting the parameters of BGBW distribution in Eq. (10). It has been considered a linear expansion profile, that is, m = 1 in Eq. (11), and the parameters t f /t r , T and β s are taken from data of π +− spectra at √ s = 2.76 TeV for each centrality class. The fit results are presented in Tab. I. The values obtained for t f /t r are very close to those ones in [31]. Such quantity lowers in terms of the centrality, indicating a higher relaxation time in peripheral collisions. This fact was understood as a result of initial distribution closer to equilibrium in more central collisions. Fig. 2 displays the resulting curve compared with data of p T spectra from π +− in [38]. The dotted and dashed lines represent the two contributions in Eq. (15), namely the hard initial distribution and the distribution of produced particles in equilibrium, respectively. It can be realized that for more central collisions (up to 10 − 20%) the region of small p T is dominated by BGBW-like thermal spectra, while for p T 4 the leading mechanism is that one of produced particles in the hard initial collision through minijet yield initiated by gluons. For larger centralities, the contribution f eq becomes smaller even in the region of small p T and then the nuclear effects are basically determined by the shadowing of gluon distribution in the initial state. Fig. 3 shows the nuclear modification factor. The pp cross section has been calculated in [18]. For the sake of consistency, the values of T AB that were employed are the same as in [38] for the determination of the experimental nuclear modification. In more central collisions, the modification on spectra is more intense than what is expected, which is explained by the shadowing seen in Fig. 1. It indicates that in this case particle production is caused by a distinct mechanism. In more peripheral collisions, the decreasing of cross section in P bP b collisions with respect to the pp one at small p T h is smaller and more compatible with what is expected if one considers only the shadowing of gluon distribution. The usual treatment for R AA is based on pQCD, in which effects of energy loss are absorbed in medium-modified parton fragmentation function in a dynamically expanding medium. There are several prescriptions on how to include radiative energy loss and some of them also introduce a collisional contribution [57][58][59][60]. Our predictions agree with these theoretical approaches at large p T and describe correctly the peak at low p T .
The inclusive multiplicity of charged hadrons dN/dy was also calculated by integrating the spectra over p 2 T h . In order to compare the pion distribution with data of dN/dy at y = 0 for charged hadrons, we have taken into account a correction of 5% relative to the contributions of charged kaons and protons. Fig. 4 presents our results compared against data from ALICE [61] as a function of the participants number, N part . The line is an interpolation among the results from each centrality class. Interestingly, one may note that f in has a slower increase in terms of N part , whereas f eq grows rapidly in the central region. Fig. 5 displays the relation between the temperature and the nuclear overlap for each centrality. While the values obtained for T and β are close to those ones taken from fits using only the BGBW model as in [38] for central collisions, (i.e., T ∼ 0.1 GeV and β ∼ 0.6 are practically constants) in peripheral collisions the values of T tend to be lower. Besides, since T and β are anti correlated, we have an increasing of β in more peripheral collisions. These differences occur due to the fact that in more central collisions the second term in Eq.
(15) is the leading one at small p T h , whereas for more peripheral collisions its contribution is much smaller. This effect can be interpreted in the following way: the first one is that the nuclear modification can be explained with good approximation only by the effect of nuclear shadowing. In this case, the fit of BGBW distribution parameters have large uncertainties. Instead, if one assigns a physical meaning for the temperature decreasing that occurs in more peripheral collisions, the possible interpretation is that collective expansion that takes place in there is similar to what happens in more central collisions for smaller collisions energies in which T is signifi- cantly lower. This fact may occur since the mean number of collisions, N coll = σ in ( √ s) T AB , grows for more central collisions and higher energies. In [35] it is pointed out that the temperature obtained from fits utilizing a Boltzmann-like distribution enhances in terms of the energy initial density (it depends on N part and √ s) until it reaches a bound, and this was understood as a QGP phase transition temperature for hadrons. In [62] it is shown that the dependence of T in terms of √ s can be parametrized as T = T lim 1 − 1 A+Be x . This form can be applied for our case. Fig. 5 presents the fitted line for this function with respect to T AB .

IV. SUMMARY AND CONCLUSIONS
In this work, we investigated the role of nuclear shadowing incorporated in the gluon distribution through the saturation/CGC formalism applied to the spectra of produced gluons in heavy ion collisions at high energies. Through the Boltzmann equation formalism within the relaxation time approximation, the contributions from the hard initial distribution and from particle production after QGP formation are separated. The former is obtained from nuclear modifications at initial state taking into consideration the shadowing of gluon distribution proposed previously to describe the p T spectra in pp collisions. These modifications were incorporated by the multiple scattering formalism in the color dipole picture. The second part of the spectra considers effects of plasma formation until the freeze-out at the final state and has been parametrized by BGBW distribution having parameters fitted from ALICE data. We verified that the inclusion of shadowing introduces modifications in the fitted parameters, especially in more peripheral collisions where the nuclear effects can be mostly explained by the modifications of gluon distribution at the initial state.