Hadron production in d-Au, p-Pb, Pb-Pb and Xe-Xe collisions at RHIC and LHC: onset of color glass condensate in light of an analytical solution of BK equation

TL;DR

This paper presents an analytical solution to the BK equation that effectively describes hadron production data from RHIC and LHC, supporting the onset of the color glass condensate in various collision systems.

Contribution

It introduces a simple, parameter-efficient analytical solution to the BK equation, enhancing understanding of gluon saturation effects in high-energy nuclear collisions.

Findings

Good description of RHIC and LHC data on differential yields

Quantitative predictions for nuclear modification factors in multiple collision systems

Supports the onset of color glass condensate at high energies

Abstract

An exact analytical solution to the Balitsky-Kovchegov equation is proposed incorporating Mellin transformation with saddle point approximation. The theoretical significance of this solution is its simplicity in form which requires very less number of parameters. Using this solution along with the idea of color glass condensate, a very good description of RHIC and LHC data on differential yield is obtained for d-Au, p-p and p-Pb collisions. Quantitative prediction for the nuclear modification factors for p-Pb, Pb-Pb and Xe-Xe collisions relevant for very recent LHC run 2 is yielded.

Equations

$${\partial }_Y\mathcal{N}={\bar{\alpha }}_s\chi (k^2)\mathcal{N}-{\bar{\alpha }}_s{\mathcal{N}}^2,$$

$${\partial }_Y\widetilde{\underline{\mathcal{N}}}=-{\bar{\alpha }}_s\tilde{\chi }(\gamma )\widetilde{\underline{\mathcal{N}}}-2\pi {\bar{\alpha }}_s\delta (i\gamma ),$$

$$\tilde{\chi }(\gamma )$$

$$=\pi \cot (\pi \gamma )-2{\mathcal{H}}_{-\gamma }^{(1)}$$

$${\tilde{\chi }}^{\prime }(\gamma )=-{\pi }^2\csc (\pi \gamma {)}^2+{\pi }^2/3-2{\mathcal{H}}_{-\gamma }^{(2)}.$$

$$\underline{\mathcal{N}}$$

$$\left[\frac{2\pi {\bar{\alpha }}_s\delta (i\gamma )}{\tilde{\chi }(\gamma )}+e^{\tilde{\chi }(\gamma )(y_0-Y)}\left({\widetilde{\underline{\mathcal{N}}}}_0-\frac{2\pi {\bar{\alpha }}_s\delta (i\gamma )}{\tilde{\chi }(\gamma )}\right)\right],$$

$$\widetilde{\underline{\mathcal{N}}}(1/2+i\nu ,y_0)\approx \widetilde{\underline{\mathcal{N}}}(1/2,y_0)e^{\frac{\nu \partial \ln \widetilde{\underline{\mathcal{N}}}(1/2+i\nu ,y_0)}{\partial \nu }{|}_{\nu =0}}.$$

$$\underline{\mathcal{N}}=$$

$$+\frac{{\bar{\alpha }}_s}{\lambda +1/8\Omega }\{1-e^{(\lambda +1/8\Omega )(y_0-Y)}\},$$

$$e^{(\lambda +1/8\Omega )(y_0-Y)}{|}_{Y\to y_0}\to 1.$$

$$\frac{d{N}_h^{pA\to hX}}{d\eta d^2{p}_t}=$$

$$\times {\tilde{N}}_F(x^{\prime },p_t^2)D_{h/q}(z,p_t^2)$$

$$\displaystyle+\frac{x_{F}}{z}f_{g/p}(x_{F}/z,p_{t}^{2})\tilde{N}_{A}(x^{\prime},p_{t}/z)D_{h/g}(z,p_{t}^{2})\bigg{]}$$

$${\tilde{N}}_{F(A)}(x,k)=\int d^2\text{r}{e}^{-i\text{ k.r}}(1-{\mathcal{N}}_{F(A)}(r,Y)),$$

$$\tilde{\mathcal{N}}_{A}(r,Y=0)=1-\text{exp}\bigg{[}-\frac{r^{2}Q_{s0}}{4}\ln\left(\frac{1}{\Lambda r}+e\right)\bigg{]}$$

$$R_{pA}=\frac{d^2{N}^{pA}/d\eta d{p}_t}{⟨T_{pA}⟩d^2{\sigma }^{pp}/d\eta d{p}_t},$$

$$\frac{d{N}^{ch}}{dyd{p}_t}=\frac{2\pi K}{C_F}\frac{{\alpha }_s}{p_t^2}\int d{k}_T^2{\varphi }_{A_1}(x_1,k_t^2){\varphi }_{A_2}(x_2,(p-k{)}_t^2),$$

$$\varphi (y,k)=\int \frac{d^{2}r}{2\pi r^2}{e}^{ir.k}\mathcal{N}(y,r).$$

Full text

Heavy ion collisions performed previously at RHIC BRAHMS collaboration et al. (2004); STAR collaboration et al. (2006) and recently at LHC CMS collaboration et al. (2018); ALICE collaboration et al. (2018a, b) probe the color glass condensate (CGC) Weigert (2005) regime of QCD governed by the nonlinear effects in the colliding nuclei wave functions as well as gluon saturation. Pb-Pb collision at LHC opens up an excellent opportunity to study the hot-dense hypothetical QCD medium arguably called quark gluon plasma (QGP). Various theoretical frameworks Kharzeev and Levin (2001); Kharzeev et al. (2005); Armesto et al. (2005a) that are modeled to study these QCD processes, primarily based on the concept of saturation scale ($Q_{sA}$) that characterizes the onset of nonlinear coherent phenomena. The growth of the saturation scale with increasing energy is yielded by the Balitsky-Kovchegov (BK) nonlinear evolution equation Balitsky (1996); Kovchegov (1999). This portrays a direct link between the gluon saturation dynamics and the experimentally measured hadron yields.

The BK equation in position space is characterized by the dipole scattering amplitude, $\mathcal{N}(r,Y)$ on a target, where the target is specified by the initial condition. For large and homogeneous targets, Fourier transformation of the scattering amplitude to the momentum space yields

[TABLE]

where $\bar{\alpha}_{s}=N_{c}\alpha_{s}/\pi$ and the quantity $\mathcal{N}(k,Y)=\int_{0}^{\infty}\frac{d^{2}r}{2\pi r^{2}}e^{ir.k}\mathcal{N}(r,Y)$ is related to unintegrated gluon distribution (UGD) of the target. The operator $\chi$ is the well known LL BFKL kernel in momentum space. The nonlinear term in (1) is not acted on by the kernel operator and thereby posses a particularly simple representation.

To solve the BK equation (1) in Mellin space we recast (1) in terms of its inverse function $\underline{\mathcal{N}}=\mathcal{N}^{-1}$ such that the semilinear PDE transforms into a linear one and then diagonalize the same by Mellin transformation

[TABLE]

where $\delta(i\gamma)$ is a Dirac-delta function. The Mellin transformation of the scale invariant kernel $\chi$ in (2) can be formulated in terms of the first order harmonic number $\mathcal{H}_{-\gamma}^{(1)}$ as follows,

[TABLE]

implies,

[TABLE]

Inverse Mellin transformation on the solution of (2) yields

[TABLE]

with $\tilde{\underline{\mathcal{N}}}_{0}=\tilde{\underline{\mathcal{N}}}(k,y_{0})$ where the initial condition $y_{0}$ is chosen large enough to ensure the validity of BK equation. An arbitrary scale $k_{0}$ is introduced for dimensional purpose. Note that $\tilde{\chi}(\gamma)$ is symmetric around $\gamma=1/2$ and along the contour $\gamma=1/2+i\nu$, with $-\infty<\nu<\infty$ has its maximum at the ”saddle” point $\gamma=c=1/2$. This indicates that the $\gamma\sim 1/2$ region dominates the small-$x$ behavior of the $N(k,Y)$ which seem to recover the familiar LL BFKL characteristic function behavior. To produce the iterative solution of (2) we therefore set $\gamma=1/2+i\nu$ in (5) and expand $\tilde{\chi}(\gamma)$ about $\nu=0$ i.e. $\tilde{\chi}(1/2+i\nu)=\lambda-1/2\Omega\nu^{2}+\mathcal{O}(\nu^{2})$ with $\lambda=\bar{\alpha}_{s}4\ln 2$ and $\Omega=\bar{\alpha}_{s}28\zeta(3)$, $\zeta$ being Reimann zeta function. We also need to expand $\tilde{\underline{\mathcal{N}}}(1/2+i\nu,y_{0})$ around $\nu=0$,

[TABLE]

Substituting these expansions back to (5) and performing $\nu$ integration along with various Dirac-delta operations one may arrive at

[TABLE]

where we define an arbitrary scale $\ln\tilde{k_{s}}^{2}=i\frac{d}{d\nu}[\ln\underline{\tilde{\mathcal{N}}}(y_{0},1/2+i\nu)]|_{\nu=0}$. Above equation serves as the analytical solution of (2), however, the original solution $\mathcal{N}$ of the BK equation (1) would be the immediate reciprocal of (7) as we recall $\underline{\mathcal{N}}=\mathcal{N}^{-1}$ from our assumption.

A few comments on the solution (7) are in order. The first term of (7) forecasts the characteristic $\sim e^{\lambda Y}$ behavior of $\mathcal{N}$, so does UGD, similar to LL BFKL evolution. On the other hand, the second term modulates the power growth arising from the linear contribution thereby suppressing the gluon evolution. Interestingly, at pre asymptotic rapidities ($Y\sim y_{0}$), the nonlinear contribution vanishes since the driving term

[TABLE]

Thus at not so large $Y$, the BK solution (7) recovers usual linear BFKL behavior. Surprisingly the nonlinear contribution in the solution is irrelevant of the momentum scale, $k^{2}$ and it depends particularly on the rapidity only.

Equation (7) suggests the exact asymptotic solution of the BK equation and thus geometrical scaling appears as a universal property of these kinds of nonlinear equations. However, the detailed phenomenological studies on the gluon evolution in context of geometrical scaling is beyond the scope of this letter. Instead, in this work, we try to address a more realistic aspect of the solution relating our theory to the ongoing collider physics experiments.

In high energy asymmetric collisions, there are at least two particular approaches leading to hadron production. The first one is the $k_{t}$ factorization formalism, appropriate for central rapidity region, where both the target and the projectile are characterized in terms of UGD. But at more forward rapidities, $k_{t}$ factorization fails to grasp the dominant contribution to the cross section coming from valence quarks of the projectile. In this situation, so called hybrid formalism is most pronounced where the large-$x$ degrees of freedom of the projectile is characterised by usual collinear gluon distribution. In hybrid formalism, the differential cross section for forward hadron production with transverse momentum, $p_{t}$ and pseudo rapidity $\eta$, in proton-nucleus collisions Dumitru et al. (2006); Albacete and Marquet (2010) is given by

[TABLE]

where the kinematics are $x_{F}=\sqrt{m_{h}^{2}+p_{t}^{2}}/\sqrt{S_{NN}}$ and $x^{\prime}=x_{F}e^{-2y_{h}}/z$ with $\sqrt{S_{NN}}$ being the collision energy per nucleon. For highly relativistic particles the difference between rapidity $y_{h}$ and pseudorapidity $\eta$ can be neglected, $\eta\approx y_{h}$. The projectile is characterized by the collinear PDFs $f_{i/p}$, however its partons eventually acquire a large transverse momentum $p_{t}$ after interaction with the small-$x$ fields of the saturated nucleus causing multiple scattering which is described by the UGDs in the adjoint or fundamental representation $\tilde{N}_{A}(F)$. The term $D_{h/i}$ refer to the final-state hadron fragmentation function which characterizes the hadronization of the scattered parton. The lack of impact parameter dependence in the calculation is accounted for by the normalization factor C.

The nuclear UGDs, $\tilde{N}_{F(A)}$ entering (8) are given by the two dimensional Fourier transform of the imaginary part of the forward dipole scattering amplitude in the fundamental (F) and adjoint (A) representation

[TABLE]

where $r$ is the dipole size and rapidity $Y=\ln(x_{0}/x)$. The small-$x$ dynamics of the dipole amplitudes $\mathcal{N}_{F(A)}(r,Y)$ in (9) is given by the solution of BK equation (7) (inverse) but in momentum space, however, it is straightforward to switch between position and momentum domains via Fourier transform. The initial conditions for the evolution of dipole-nucleus amplitude are taken from the semiclassical Mclerran-Venugopalan (MV) model McLerran and Venugopalan (1998),

[TABLE]

where $Q_{s0}$ is the initial saturation scale probed by quarks. The constant $e$ serves as an infrared regulator and $\Lambda\approx 0.24$ GeV. The scattering amplitude in fundamental representation, $\mathcal{N}_{F}$ can also be parametrized as in (10) but with the replacement $Q_{s}^{2}\rightarrow Q_{s}^{2}C_{F}/C_{A}=4/9Q_{s}^{2}$ Kharzeev et al. (2004).

With this setup, we obtain a very good description of RHIC and LHC run 2 data for hadron yields of different hadron species in p-p, d-Au and p-Pb collisions. The invariant differential yields in d-Au collision at RHIC energy $\sqrt{S_{NN}}=200$ GeV and pseudo-rapidity $\eta$ = 2.2 and 3.2 for negatively charged hadrons ($h^{-}$) (BRAHMS collaboration BRAHMS collaboration et al. (2004)) and $\eta$ = 4 for neutral pions ($\pi^{0}$) (STAR collaboration STAR collaboration et al. (2006)) is shown in Fig. 1. In analogy to Armesto et al. (2005b), we obtain a quark saturation scale $Q_{s0}\approx 0.43~{}\text{GeV}^{2}$ for gold nucleus at $x_{0}=0.01$. Corresponding gluon saturation scale $Q_{s0}^{2}$ is found to be $Q_{s0}^{2}\approx 0.97~{}\text{GeV}^{2}$. For entire calculations performed with hybrid formalism, we adopt the DSS NLO fragmentation functions de Florian et al. (2007) and CT14 NLO collinear PDFs Dulat et al. (2016) as default. Our results show a good agreement with RHIC data. However, for the better description of neutral pion data, a normalisation factor $C=0.6$ is multiplied, while charged hadron yields are sketched without any normalisation (i.e. C=1). Considering the sizable uncertainty in the data, the precise value of the normalisation constant does not have much meaning. In Fig. 1, the $\pi^{0}$ and $h^{\pm}$ yields for p-Pb collisions and $h^{\pm}$ yield for p-p collisions measured in ALICE at $\sqrt{S_{NN}}=5.02~{}\text{TeV}$ ALICE collaboration et al. (2018a) are also shown. Corresponding pseudo rapidity windows for $h^{\pm}$ and $\pi^{0}$ are $|\eta|<0.8$ and $|\eta|<0.4$ respectively. Here for all calculation the normalisation is set to 1. The theoretical measurement are clearly consistent within the uncertainties over the entire $p_{t}$ range.

The modification of forward particle production is quantified with the nuclear modification factor,

[TABLE]

where $N^{pA}$ is the charged particle yield in pA collisions and $\sigma^{pp}$ is the particle production cross-section in p-p collision. The average nuclear overlap function, $\braket{T_{pA}}=\braket{N_{coll}}/\sigma_{inel}^{NN}$ depends on the collision centrality determined from the Glauber model Loizides et al. (2018). For instance, at $\sqrt{S_{NN}}=5.02$ TeV, $\braket{T_{pPb}}\approx 0.098\quad\text{mb}^{-1}$ for minimum bias. Note that computing the ratio (11) removes sensitivity to the normalization factor C. Our $R_{pPb}$ calculations scaled by $\braket{T_{pPb}}$ are displayed in Fig. 2 at $\sqrt{S_{NN}}=5.02~{}\text{TeV}$ subject to the pseudo rapidity windows $-0.8<\eta<0.8$ for neutral pions ($\pi^{0}$) and $-1.265<\eta<0.335$ for charged pions ($\pi^{+}+\pi^{-}$) agree well with the ALICE data within the uncertainties over the complete $p_{t}$ domain. Both theory and data indicate a sizable suppression around 30% at $p_{t}=1~{}\text{GeV}$, whereas the ratio $R_{pPb}$ is seen tending to unity for $p_{t}$ above 2 GeV for both neutral pions and charged pions.

We have also checked our BK analytical solution + CGC framework to testify charged hadron production particularly at mid rapidities using $k_{t}$ factorization formalism. At mid rapidities, $x_{F}/z$ becomes too small to valid hybrid formalism in (8). The $k_{t}$ factorization formula is valid for small values of $x$, for both the projectile and targets. But since it includes only gluonic degrees of freedom, this formalism is insufficient at very forward rapidities or large $p_{t}$. Within $k_{t}$ factorization approach Albacete (2007) the pseudo rapidity density of charged particle can be calculated via.

[TABLE]

where the variables $x_{1,2}=p_{t}/\sqrt{S_{NN}}e^{\pm y}$ and $K$ is accounted for the normalisation purpose. The UGD $\phi$, in (12) is related to collinear gluon distribution, $\phi(x,k)\propto\frac{d(xG(x,k^{2}))}{d^{2}k}$ and is given in terms of the dipole scattering amplitude

[TABLE]

The CMS and ALICE data CMS collaboration et al. (2018); ALICE collaboration et al. (2018b) for the nuclear modification factor $\text{R}_{\text{AA}}$ for primary charged particles in Pb-Pb ($\sqrt{S_{NN}}=5.02~{}\text{TeV}$) and Xe-Xe ($\sqrt{S_{NN}}=5.44~{}\text{TeV}$) collisions are shown in Fig. 3 for the centrality bins 0-100%, 50-70%, 10-30% and 0-5%. For minimum bias (0-100%), the ratio $R_{pA}$ is also included in Fig. 3 to sketch a quantitative comparison with $R_{AA}$. The pronounced suppression in our theoretical prediction is seen at $p_{t}<2~{}\text{GeV}$ for p-Pb system and at $8<p_{t}<100~{}\text{GeV}$ for Pb-Pb system which is consistent with CMS data. In the peripheral (50-70%) bin, $R_{AA}$ measurements at $\sqrt{S_{NN}}=2.76~{}\text{TeV}$ is also shown to examine the sensitivity of the ratio towards probe energy. The measured distribution $R_{AA}$ at 2.76 TeV and 5.02 TeV are found to be quantitatively similar at intermediate $p_{t}$, but above about 10 GeV, $R_{AA}$ at 5.02 TeV tend to be slightly smaller than the same at 2.76 TeV. In Fig. 3, $R_{AA}$ for Xe-Xe collision at $\sqrt{S_{NN}}=5.44~{}\text{TeV}$ is also sketched for two centrality bins 0-10% and 10-30%. The theoretical prediction seem to be in rough agreement with CMS data at moderate $p_{t}$, however, towards high $p_{t}$, our theory tend to predict much stronger suppression compared to the data. From this analysis it is clear that our theory holds good for minimum bias measurements, while it might be inappropriate for most central events. This is expected as our theory is impact parameter $b$ independent and thus one have to improve the theory for the $b$ dependence of the UGD to study different centrality bins which is beyond the scope of this letter.

In summary, we have proposed an analytical solution to BK equation incorporating Mellin transformation with saddle point approximation. The main key feature of this solution is its simplicity in form with relatively very less number of parameters, which might be handy for investigation of various characteristics of such parton evolution processes. At pre asymptotic $y$, the solution is found to reconstruct the characteristics linear BFKL solution. The feasibility of the solution is tested at both RHIC and LHC energies. As per the good consistency between the theory and data is concerned, the analytical solution of BK equation could be a very reliable framework for exploring quark gluon plasma at LHC as well as future colliders.