Suppression vs enhancement of heavy quarkonia in pA collisions

TL;DR

This paper models heavy quarkonia production in pA collisions using the dipole approach, revealing that nonzero heavy quark separation dominates nuclear effects and explaining experimental suppression patterns across energies.

Contribution

It introduces a dominant higher twist correction based on heavy quark separation and proposes a new mechanism for charmonium enhancement, aligning with experimental data.

Findings

Nonzero heavy quark separation is the main source of nuclear effects.

Nuclear suppression remains stable from RHIC to LHC energies.

The dipole model explains the suppression of $ ext{psi}(2S)$ relative to J/psi.

Abstract

We describe production of heavy quarkonia in pA collisions within the dipole approach, assuming dominance of the perturbative color-singlet mechanism (CSM) in the $p_T$-integrated cross section. Although accounting for a nonzero heavy $Q$-$\bar Q$ separation is a higher twist correction, usually neglected, we found it to be the dominant source of nuclear effects, significantly exceeding the effects of leading twist gluon shadowing and energy loss. Moreover, this contribution turns out to be the most reliably predicted, relying on the precise measurements of the dipole cross section at HERA. The nuclear suppression of quarkonia has been anticipated to become stronger with energy, because the dipole cross section steeply rises. However, the measured nuclear effects remain essentially unchanged within the energy range from RHIC to the LHC. A novel production mechanism is proposed, which enhances the charmonium yield. Nuclear effects for the production of $J/\psi$, $\psi(2S)$, $\Upsilon(1S)$ and $\Upsilon(2S)$ are calculated, in agreement with data from RHIC and LHC. The dipole description offers a unique explanation for the observed significant nuclear suppression of $\psi(2S)$ to $J/\psi$ ratio, related to the nontrivial features of the $\psi(2S)$ wave function.

Figures

Equations

$$l_c=\frac{1}{q_L}=\frac{2E_{\bar{c}c}}{M_{\bar{c}c}^2}.$$

$$l_f\sim \frac{2E_{\bar{c}c}}{M_{\psi (2S)}^2-M_{J/\psi }^2}\gg l_c$$

$$A(z_1,z_2)=\int d^2{r}_1{d}^2{r}_2{\Psi }_f^{†}({\vec{r}}_2)G({\vec{r}}_2,z_2;{\vec{r}}_1,z_1){\Psi }_{in}({\vec{r}}_1),$$

$$i\frac{\partial }{\partial z_2}G\left(z_2,{\vec{r}}_2;z_1,{\vec{r}}_1\right)=$$

$$\left[\frac{m_c^2-{\Delta }_r}{2E_{\bar{c}c}{\alpha }_c{\bar{\alpha }}_c}+V\left(r,z_2\right)\right]G\left(z_2,{\vec{r}}_2;z_1,{\vec{r}}_1\right),$$

$$A(z_1,z_2)\Rightarrow A_0(z_1,z_2)=\int d^{2}r{\Psi }_f^{†}(\vec{r}){\Psi }_{in}(\vec{r}),$$

$$ϵ(x_2,\Delta z)=\frac{{\left|A(z_1,z_2)\right|}^2}{{\left|A_0(z_1,z_2)\right|}^2},$$

$$x_{1,2}=\frac{M_T}{\sqrt{s}}{e}^{\pm y},$$

$$E=\frac{M_T^2}{2m_N{x}_2}.$$

$$G({\vec{r}}_2,z_2;{\vec{r}}_1,z_1)=\frac{N}{2\pi \sinh (\Omega \Delta z)}$$

$$N^2$$

$$\Omega$$

$$\displaystyle G(\vec{r}_{2},z_{2};\vec{r}_{1},z_{1})\Bigr{|}_{free}$$

$${\bar{c}}^i{c}_j+N\to {\bar{c}}^k{c}_l+X,$$

$$S(r,z)$$

$$O(r,z)$$

$$S(r,z)$$

$$O(r,z)$$

$$n_{coll}^{\bar{c}c}(r,B)={\sigma }_{\bar{c}c}(r,E_{\bar{c}c})T_A(B),$$

$$T_A(B)=\underset{-\infty }{\overset{\infty }{\int }}dz{n}_A(B,z),$$

$$R_{pA}^{J/\psi }(s,y)$$

$${\sigma }_4(\rho ,{\alpha }_g)=\frac{1}{2}\left[{\sigma }_{gg}(\rho )+{\sigma }_{gg}({\alpha }_g\rho )+{\sigma }_{gg}({\bar{\alpha }}_g\rho )\right],$$

$${\sigma }_{gg}(\rho )=\frac{9}{4}{\sigma }_{\bar{q}q}(\rho ).$$

$$x_{1,2}=\frac{\sqrt{M_{\bar{c}c}^2+p_T^2}}{\sqrt{s}}{e}^{\pm y_{\bar{c}c}},$$

$${\sigma }_{\bar{q}q}(r,x_2)={\sigma }_0\left\{1-\exp \left[\frac{{\pi }^2{r}^2{\alpha }_s({\mu }^2)x_{2}g(x_2)}{3{\sigma }_0}\right]\right\},$$

$$W(\vec{\rho },\vec{r},{\vec{r}}^{\prime })=\frac{d{\sigma }_{pp}^{J/\psi }}{dy{d}^2\rho d^{2}r{d}^2{r}^{\prime }}{\left[\frac{d{\sigma }_{pp}^{J/\psi }}{dy}\right]}^{-1},$$

$$R_{pA}^{(1N)}(s,y)$$

$$S_A^{(1N)}(B,z_1,\vec{\rho },\vec{r},{\vec{r}}^{\prime })$$

$$T_{-}(B,z)$$

$$T_{+}(B,z)$$

$$R_{pA}^{(2N)}(s,y)=\frac{{\sigma }^{(2N)}(pA\to J/\psi X)}{A\sigma (pp\to J/\psi X)},$$

$$\sigma (gp\to (\bar{c}c{)}_{\{8^{-}\}}X)=\sum _{\mu ,\bar{\mu }}\underset{0}{\overset{1}{\int }}d\alpha d{\alpha }^{\prime }\int d^{2}r{d}^2{r}^{\prime }$$

$${\Sigma }_{g\to \{8^{-}\}}\approx \frac{5}{8}\left[{\sigma }_{\bar{q}q}\left(\frac{\vec{r}+{\vec{r}}^{\prime }}{2}\right)-{\sigma }_{\bar{q}q}\left(\frac{\vec{r}-{\vec{r}}^{\prime }}{2}\right)\right].$$

Full text

I Introduction

Inelastic interactions of a heavy quark pair propagating through a nucleus is a higher twist effect, $\sim 1/m_{c}^{2}$, which is therefore usually neglected, while leading twist gluon shadowing is believed to be the main source of nuclear suppression at high energies. However, a considerable nuclear suppression of $J/\psi$ production in $pA$ collisions has been observed in the pioneering measurements serpukhov ; na3 , even though the energy range of these experiments was too low to explain the observed nuclear effects by gluon shadowing. These data provided the first evidence for the importance of higher twist effects, which certainly remain essential at higher energies and should contribute to the strong nuclear suppression observed in $pA$ collisions at Fermilab e866 , RHIC phenix-psi and LHC alice1 ; alice2 . Although higher-twist effects as a possible explanation of the observed nuclear suppression of $J/\psi$ was proposed in pire , no numerical evaluation was done.

In what follows for concreteness we consider charmonium production, and mainly $J/\psi$, unless otherwise stated. However, the developed techniques will be also applied to the calculation of nuclear effects in the production of radial excitations and bottomium states.

Charmonium suppression related to the non-zero size, $r\sim 1/m_{c}$, of a perturbatively produced $\bar{c}c$ dipole, is a higher twist effect, which vanishes in the limit of high quark masses. Quantitatively, however, it turns out to be the main contributor to the nuclear effects in charmonium production observed so far nontrivial ; rhic-lhc . At this point we should emphasise that this higher twist effect is the best known part of nuclear effects. The dipole cross section has been thoroughly measured in deep-inelastic scattering (DIS) at HERA, as function of the dipole energy and size. Therefore, the higher-twist part of dipole attenuation in nuclear matter, which is responsible for charm nuclear shadowing, is pretty well known, and leaves not much room for other mechanisms, when is compared with data nontrivial ; rhic-lhc . On the other hand, leading twist gluon shadowing, which makes nuclear medium more transparent for dipoles, has been poorly fixed by data so far, ranging from a very weak florian1 ; florian2 up to dramatically strong effect eps08 , even breaking the unitarity bound klps .

These effects lead to reduction of the $J/\psi$ production rate, while the magnitude of shadowing (both, leading and higher twist), as well as the break-up dipole cross section, steadily rise with energy. Therefore, it looks natural to anticipate a stronger suppression of $J/\psi$ produced in $pA$ collisions at the LHC compared with RHIC, as was predicted in rhic-lhc ; andronic . However, the measurements alice1 ; alice2 , unexpectedly revealed energy independent magnitude of $J/\psi$ suppression, which remains unchanged through the huge energy range between RHIC and LHC. This contradiction creates a serious challenge, because as was mentioned above, the dipole phenomenology is well fixed by HERA data, leaving little freedom in its predictions. In spite of the large uncertainly in the gluon shadowing case, it cannot reduce the problem, because its magnitude also rises with energy.

Here we identify a novel mechanism, which enhances charmonium production and explains the observed anomalous energy dependence. This mechanism was proposed and employed in hkz for the explanation of the EMC experiment puzzling data emc on nuclear photoproduction of $J/\psi$. The observed nuclear enhancement was related to non-Glauber double color-exchange interactions with different bound nucleons. In fact, multiple color-exchange interactions of a dipole propagating through a nuclear medium lead to a non-vanishing survival probability of the dipole, and even to an enhancement in specific channels. This is demonstrated in Sect. II, based on the results presented in Appendix A.

Nevertheless, the opacity expansion shows that the mean number of color-exchange interactions is rather small even in heavy nuclei, although it rises with energy. Still, the main contribution to $J/\psi$ production is expected to be provided by the single color-exchange interaction considered in Sect. III. The cross section on a nucleon tends to cancel in the nucleus-to-proton ratio, but the nuclear attenuation factor depends on the features of the parton ensemble propagating through the nucleus. Therefore the description of $J/\psi$ production in $pp$ collisions is essential within the dipole approach, since it allows to calculate the distribution function of the produced partons in impact parameter space. The details of the calculations are presented in Appendix B. The next term of opacity expansion, the double color-exchange interaction, is described in Sect. IV. The specific challenge here is the calculation of nuclear attenuation factors for $\bar{c}c$ pairs in certain color states. We found the correction $R^{(2N)}$ to the nuclear ratio $R_{pA}$ to be significant.

Other nuclear effects are also included in the calculations. Gluon shadowing corrections, evaluated in Sect. V, are found to be negligibly small at the RHIC energy, but rather significant at the energy of LHC, especially at forward rapidities. Energy loss corrections are considered and introduced in the calculation in Sect. VI. The nonperturbative source of energy loss, related to the energy sharing problem at forward rapidities, occurs on a soft scale and brings major corrections to the nuclear effects, as is described in Sect. VI.1. The perturbative mechanism of energy loss, described in Sect. VI.1, is related to the phenomenon of saturation, which generates a new scale, the saturation momentum, or nuclear broadening. We found the related energy loss to be a quite weak effect, being strongly suppressed by the ratio of the saturation scale to the quarkonium mass squared. This suppression has been missed in previous evaluations, which grossly overestimated this effect of energy loss.

Special interest has always been paid to nuclear effects in the production of radial excitations, considered in Sect. VII. The quarkonium wave function participating in the convolution with the produced $\bar{c}c$ wave packet, has a node, which leads to a partial compensation between small and large dipole separations. Nuclear color filtering modifies the convolution and can lead to illuminating effects, as was found in the photoproduction of $\psi(2S)$. The dynamics of hadroproduction is more involved and we arrived at a stronger suppression of $\psi(2S)$ compared to $J/\psi$. Nuclear effects in $\psi(2S)$ production is a sensitive test of the dipole description of the production mechanism. It provides a unique explanation of the strong suppression of the $\psi(2S)$ to $J/\psi$ ratio in $pA$ collisions.

The developed dipole description of nuclear effects in charmonium production can be easily extended to heavier quarkonia. In sect. VIII we perform calculations for the production of $\Upsilon(1S)$ and $\Upsilon(2S)$, in good accord with available data.

II Propagation of $\bar{c}c$ dipoles in nuclear medium

II.1 Characteristic length scales

Two general amplitudes of $\bar{c}c$ production at different points separated by longitudinal distance $\Delta z$, have a relative phase shift $\Delta\phi=q_{L}\Delta z$ in the nuclear rest frame, where the longitudinal momentum transfer is $q_{L}=M_{\bar{c}c}^{2}/2E_{\bar{c}c}$. Correspondingly, the longitudinal length scale $l_{c}=1/q_{L}$, usually called coherence length brod-mueller ; kz91 , reads,

[TABLE]

If the coherence length exceeds the nuclear dimension, one cannot localize the coordinate of the $\bar{c}c$ pair production, in which case the pair propagates through the whole nucleus. This regime occurs at the energies of RHIC and LHC (except for large negative rapidities).

The $\bar{c}c$ dipole produced with small transverse separation $r_{T}\sim 1/m_{c}$, expands and eventually forms the charmonium wave function on a much longer length scale, called formation length brod-mueller ; kz91 ,

[TABLE]

where the masses in the denominator correspond to the first radial excitation $\psi(2S)$ and the $J/\psi$. This can be interpreted in terms of the uncertainty principle as the time interval required to disentangle between the two hadronic masses, while the originally created $\bar{c}c$ pair has no certain invariant mass and no wave function.

II.2 Fluctuating dipoles

It is clear that at sufficiently high energies the dipole separation does not fluctuate during propagation through the nucleus due to Lorentz time dilation. In this regime the calculations are significantly simplified, so we intend to figure out the kinematic constraints for employing such a ”frozen” regime.

The evolution of a $\bar{c}c$ dipole propagating through a medium can be described summing up all possible trajectories of the quarks between the initial and final states. The amplitude of dipole propagation between longitudinal coordinates $z_{1}$ and $z_{2}$, with initial and final transverse separations $\vec{r}_{1}$ and $\vec{r}_{2}$ respectively, is given by the matrix element of the Green function

[TABLE]

where $\Psi_{in}(\vec{r}_{1})$ and $\Psi_{f}(\vec{r}_{2})$ are the initial and final $\bar{c}c$ distribution amplitudes respectively.

The Green function satisfies the two-dimensional light-cone equation kz91 ; kst1 ; kst2 ; klss ; kpss-psi ,

[TABLE]

with the boundary condition, $G\left(z_{2},\vec{r}_{2};z_{1},\vec{r}_{1}\right)_{\Delta z=0}=\delta(\vec{r}_{2}-\vec{r}_{1})$. Here $\alpha_{c}$ and $\bar{\alpha}_{c}=1-\alpha_{c}$ are the fractional light-cone momenta of $c$ and $\bar{c}$ respectively. In what follows we fix $\alpha=1/2$ because the charmonium wave function strongly peaks at this value kth-psi ; klss ; kpss-psi . The real part of the light-cone potential $\mbox{Re}\,V(r)$ describes the binding effects, while $\mbox{Im}\,V(r,z)$ is related to the absorption effects, i.e. multiple inelastic interactions of the dipole with the medium.

The goal of this section is to figure out the kinematic range of validity of the ”frozen” approximation, which corresponds to the high-energy limit, where the formation length Eq. (2) is much longer than the path length of the dipole in the medium, $l_{f}\gg\Delta z$. In this ”frozen” dipole regime the Green function approaches the limit $G(\vec{r}_{2},z_{2};\vec{r}_{1},z_{1})\Rightarrow\delta(\vec{r}_{1}-\vec{r}_{2})$, and correspondingly the amplitude Eq. (3) takes the form,

[TABLE]

In order to quantify the deviation from the ”frozen” approximation we evaluate the ratio,

[TABLE]

where $x_{1,2}$ are the fractional light-cone momenta of the colliding gluons, $gg\to\bar{c}c$,

[TABLE]

Here $M_{T}=\sqrt{M^{2}_{\bar{c}c}+p_{T}^{2}}$, $p_{T}$ and $y$ are the transverse invariant mass, transverse momentum and rapidity (in the $NN$ collision c.m.) of the produced $\bar{c}c$ pair, respectively. Notice that the dipole energy in the nuclear rest frame is directly related to the value of $x_{2}$ ,

[TABLE]

Anticipating that the validity of the ”frozen” approximation means that the result is not sensitive to the details of the binding potential, we evaluate $\epsilon(x_{2},\Delta z)$ in a harmonic oscillator potential modelkz91 ; kst2 , $\mbox{Re}\,V(r)=(2\omega^{2}m_{c}^{2}/E)r^{2}$, where $\omega=(M_{\psi^{\prime}}-M_{J/\psi})/2\approx 0.3\,\mbox{GeV}$. The imaginary part is related to the absorption rate, $\mbox{Im}\,V(r,z)=C(x_{2})r^{2}\,n_{A}(z)/2$, where the nuclear density is assumed to be constant, $n_{A}=0.15\,\mbox{fm}^{-3}$, and the coefficient $C(x_{2})$, which controls the dipole cross section at small dipole separations, was calculated in broadening with the parametrization gbw of the dipole cross section. In this case Eq. (4) has analytic solution kz91 ; kst2

[TABLE]

Here

[TABLE]

With this solution we evaluated $\epsilon(x_{2},\Delta z)$, Eq. (6), fixing $\Delta z=5\,\mbox{fm}$ and using the oscillatory $J/\psi$ wave function, as well as the initial distribution function with the mean separation $\langle r^{2}\rangle\sim 1/m_{c}^{2}$. The results are depicted as function of $x_{2}$ by a solid curve in Fig. 1.

We see that the ”frozen” approximation is valid with a high precision up to rather large values of $x_{2}\sim 0.1$, and works reasonably well even at larger $x_{2}$, matching the Glauber regime. These results confirm the observation made earlier kz91 , that the nuclear effects in $J/\psi$ photoproduction remain constant down to quite low energies and are close to the results of the Glauber approximation.

It is instructive to compare this with free $\bar{c}c$ pair propagation with no binding potential and no absorption. In this case the free Green function is simplified,

[TABLE]

The corresponding ratio $\epsilon(x_{2},\Delta z)$ is depicted by a dotted curve in Fig. 1. We see that even in this extreme case of free expansion the ”frozen” approximation is still accurate up to $x_{2}\sim 0.1$, far more than is needed for the description of available data for $J/\psi$ production at RHIC and LHC. Of course at larger $x_{2}$ the result significantly deviates from the ”frozen” limit, because the quarks freely fly away from each other.

II.3 Breakup and restoration of colorless dipoles

According to the conventional wisdom, supported by eikonal-type models, the survival probability of a colorless $\bar{c}c$ dipole propagating through a nuclear medium is exponentially falling with respect to the propagation path length. This is expected to be a result of color-exchange interactions with the surrounding bound nucleons, which break-up the dipole. However, as is demonstrated below, this is not correct, a high-energy dipole has a finite survival probability even in the limit of full absorption, the so called ”black disk” regime alesha ; k-zam .

If the dipole energy is sufficiently high, the regime of ”frozen”dipoles, described above, remains valid in the medium. Indeed, multiple color-exchange interactions of the dipole with the bound nucleons do not affect the dipole transverse separation, and the interactions only change the color indices of the quark pair, leading to breakup of the dipole, which becomes colored,

[TABLE]

as is illustrated in Fig. 2.

Such interactions also destroy the target nucleons, so they occur incoherently and should be described in terms of the density matrix ${}^{k}_{l}U^{i}_{j}(\vec{r};\vec{r}^{\,\prime};z)$. The evolution of the density matrix of a high-energy dipole propagating through the nuclear matter is described in Appendix A. Here we present the results for the probabilities of production of the final dipole in either color singlet, $S(r)$, or color octet, O(r), states ($\vec{r}=\vec{r}^{\,\prime}$).

After propagation through nuclear matter of thickness, $\Delta T_{A}=\int_{z_{1}}^{z}dz^{\prime}\,n_{A}(z^{\prime})$, where $n_{A}(z^{\prime})$ is the nuclear density along the propagation trajectory, the probability of finding the dipole in a color singlet (S) or octet (O) states, reads (see derivation in Appendix A),

[TABLE]

Here $S_{in}(r)$ is the size-distribution function of the initial color-singlet dipole; $\sigma_{\bar{q}q}(r)$ is the universal dipole-nucleon cross section zkl , which depends on transverse dipole separation and implicitly on the dipole energy or Bjorken $x_{2}$ (unless specified otherwise). This cross section is difficult to predict theoretically, but it is well known from phenomenology, fitted to DIS and photoproduction data. A concrete parametrization will be specified later.

Even if the initial state is a color-octet dipole with the size distribution function $O_{in}(r)$, evolution in the medium may end up with production of either a color singlet, or octet,

[TABLE]

We see from (13) and (14) that for a large number of inelastic collisions of the $\bar{c}c$ dipole, $\sigma_{\bar{c}c}(r)\Delta T_{A}\gg 1$, the probabilities of production of color-singlet and octet states approach the universal values, $1/9$ and $8/9$ respectively, independently of the color structure of the incoming $\bar{c}c$ pair. This could be easily anticipated, since both quarks become completely unpolarized in color after multiple rotations in the color space. All possible $9$ color states ($N_{c}^{2}$) of the $\bar{c}c$ are produced with equal probabilities, and only one of them is a singlet, while the 8 others ($N_{c}^{2}-1$) are octets.

II.4 Opacity expansion

The mean number of inelastic (color-exchange) collisions of a $\bar{c}c$ dipole of transverse quark separation $r$, propagating through the nucleus, is,

[TABLE]

where the nuclear thickness function at impact parameter $B$ reads,

[TABLE]

and $n_{A}(B,z)$ is the nuclear density.

For the energy dependence of $\sigma_{\bar{c}c}(r,E_{\bar{c}c})$ we rely on parametrizations in the saturated form gbw ; bgbk ; amir for $\sigma_{\bar{c}c}(r,x)$, fitted to HERA data on the proton structure function $F_{2}(x,Q^{2})$. We are interested in rather low values of $Q^{2}\sim M_{\bar{c}c}^{2}$, for which even the simple parametrization gbw works well hikt . The value of target fractional momentum of a target gluon $x_{2}=e^{-y}\sqrt{(M_{\bar{c}c}^{2}+p_{T}^{2})/s}$, controls the magnitude of the dipole cross section. Here $y$ is the rapidity of the produced $\bar{c}c$ pair; $p_{T}$ is its transverse momentum, which is of the order of the mean value, because we are interested in the $p_{T}$-integrated cross sections.

The dipole cross section steeply rises with energy at small separations, $\sigma_{\bar{q}q}(r,x)\sim(1/x)^{0.3}$. At energy $\sqrt{s}=200\,\mbox{GeV}$ and at the measured so far rapidity range 0<y\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}2, the mean number of collisions preceding the production of the final colorless $\bar{c}c$, is $n^{\bar{c}c}_{coll}\sim 0.05-0.1$. Correspondingly, at energy $\sqrt{s}=5\,\mbox{TeV}$ and 0<y\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}3, $n^{\bar{c}c}_{coll}\sim 0.1-0.2$.

In view of such a small probability of interaction, we keep only the two lowest order terms in the opacity expansion: (i) single-step direct production kth-psi ; rhic-lhc of charmonium by the projectile gluon interacting with a bound nucleon, $gN\to\{\bar{c}c\}_{\psi}X$, with coordinates $(z,\vec{B})$, with no preceding or final state interactions; (ii) a double-step process hkz , with the production of a color-octet dipole, $gN\to\{\bar{c}c\}_{8}X$, in the first collision, and the final creation of $J/\psi$ in the second collision, $\{\bar{c}c\}_{8}+N\to J/\psi+X$.

Correspondingly, the ultimate observable to be calculated, the nucleus-to-proton ratio, gets contributions from two terms,

[TABLE]

We assume here that all cross sections are $p_{T}$-integrated.

The first term $R^{(1N)}_{pA}$, single-step production, was evaluated for the production of $\chi_{2}$ in kth-psi and for $J/\psi$ in nontrivial ; rhic-lhc . While this term alone reproduces RHIC data reasonably well, the nuclear suppression predicted in rhic-lhc for the LHC, turned out to be too strong compared with the latest measurements alice1 ; alice2 . Data show that nuclear suppression of $J/\psi$ remains nearly unchanged within the wide energy range from RHIC to LHC. This is impossible for $R^{(1N)}_{pA}$, because the dipole cross section, well constrained by HERA data, rises steeply with energy, leading to a stronger nuclear attenuation of dipoles and smaller values of $R^{(1N)}_{pA}(s,y)$ at higher energies. Therefore, the observed similarity of nuclear effects at RHIC and LHC indicates the onset of a new mechanism, which enhances $J/\psi$ production at LHC energies. A natural candidate for such a mechanism is the double-step term in (17), which indeed gives a positive contribution, which rises with energy faster than the single-step term.

At this point a word of caution is in order. The above estimates for the opacity expansion assumed the same interaction cross section for each of multiple collisions. If, however, the double-step production is dominated by the color-singlet mechanism, the term $R_{2}$ turns out to be a ratio of different mechanisms. Moreover, gluon radiation in the color-singlet model (CSM) brings an additional factor $r$ (dipole size) into the amplitude, and then the $r$-dependences of $R_{1}$ and $R_{2}$ become similar. Therefore, one can make reliable conclusions about the relative values of the two terms in Eq. (17) only after performing detailed calculations, presented below.

III Single-step $J/\psi$ production

At first glance, if we assume that charmonium is produced on a bound nucleon in the same way as on a free one (see however Sect. VI), the production cross section on a nucleon should cancel in the first term $R^{(1N)}$ of the nuclear ratio Eq. (17), as happens in Glauber-type models. However, attenuation of the projectile and produced partonic ensembles propagating through the nucleus depends on the mutual transverse separations between the partons, which are controlled by the production mechanism. Therefore, the nucleon cross section of $J/\psi$ production does not cancel out and affects the nuclear ratio $R^{(1N)}$, which becomes model dependent.

First of all, one should specify the model for $J/\psi$ production, $pp\to J/\psi X$. Currently the most successful parameter-free description of data has been achieved within the color-singlet model (CSM) proposed in csm1 ; csm2 , with further developments and applications in csm3 ; csm4 ; csm5 . Production of $J/\psi$ is treated in CSM perturbatively, as a result of glue-glue fusion resulting in production of a colorless $S$-wave $\bar{c}c$ pair and a gluon as is illustrated in Fig. 3 (left). Gluon radiation allows the $\bar{c}c$ dipole to have $S$-wave symmetric wave function (see below).

Another popular approach, called color octet model, is based on the non-relativistic QCD effective field theory nrqcd ; raju ; Fujii:2013gxa ; Ducloue:2015gfa ; Ducloue:2016pqr . The main assumption of the model is that color neutralization occurs via evaporation of soft gluons on a long time scale, of the order of the formation time Eq. (2). Such an unjustified assumption has obvious problems. The initial color-octet $\bar{c}c$ pair is produced perturbatively at a hard scale $Q^{2}\sim 4m_{c}^{2}$, with no soft gluonic field with frequencies $k_{T}<m_{c}$. The laking field is regenerated via perturbative radiation of gluons, making possible $J/\psi$ production in color octet to singlet transition $\{\bar{c}c\}_{8^{-}}\to g\{\bar{c}c\}_{1^{+}}$, which is a part of the CSM (see details and notations below). In this way the $\bar{c}c$ pair can survive as a color octet and evolve its virtuality down to low scale of the order of the inverse mean $J/\psi$ radius, and then radiate gluons non-perturbatively (color evaporation). However, the probability of scale evolution without gluon radiation, neutralizing the dipole color, is suppressed by a Sudakov-like factor, which is ignored in the color-octet model.

Moreover, the idea of preferable color neutralization at a soft scale, enhanced by a large value of the QCD coupling, does not seem to be correct either. Indeed, according to the Low theorem low the matrix element of a process with soft radiation is proportional to the process amplitude without radiation, which is impossible for $J/\psi$ production 111We thank Yuri Dokshitzer for this remark.. Besides, this model has low predictive power, because it fits the unknown parameters to the data to be explained. In view of all that, we will consider the color-octet model as a dominant mechanism of $J/\psi$ production.

Another alternative to the CSM is the possibility of producing $J/\psi$ without gluon radiation, but via $1+2$ gluon fusion, where the two gluons originate either from the beam, or target. However, evaluation of the cross section motyka results in an order of magnitude smaller production rate in comparison to the CSM. We disregard this contribution in what follows. Nevertheless, a precaution is required for $J/\psi$ production at very forward rapidities, where CSM is suppressed by the shrinking phase space for gluon radiation.

III.1 Initial state shadowing vs final state attenuation

As we already discussed in Sect. II.1, at sufficiently high energies any short time interval is subject to Lorentz time dilation and becomes long. Even a hard collision, which is characterized by a very short time scale $\tau\sim 1/Q$ in its c.m. frame, may last long time (see (1)) in the target rest frame, longer than the nucleus dimension. In this limit $J/\psi$ production can be treated as a result of interaction of the $|\bar{c}cg\rangle$ Fock component of the incoming gluon with the whole nucleus.

Formally one can derive this adding up the two amplitudes depicted in Fig. 3. The first one corresponds to the direct production of the final $S$-wave colorless $\bar{c}c$ pair symmetric in spacial and spin variable, denoted by $\{1^{+}\}$. Another contribution, depicted by the right picture in Fig. 3, contains diffractive on-mass-shell production of the projectile gluon fluctuation $g\to\bar{c}cg$, preceding the color-exchange interaction. In order to end up with the production of a $J/\psi$, the $\bar{c}c$ pair in this fluctuation should be a $P$-wave color octet state, asymmetric in spacial-spin variables, which we denote by $\{8^{-}\}$ kt-hf . This color octet pair undergoes color-exchange interactions with the same bound nucleon, as in the first term of the amplitude, and switches to the final colorless $\{1^{+}\}$ state.

While the color-exchange interaction occurs on different nucleons incoherently, the diffractive production on different nucleons is a coherent process. If the coherence length (inverse longitudinal momentum transfer) is much longer than the nucleus radius, the result is equivalent to interaction of a $|\bar{c}cg\rangle$ fluctuation with the whole nucleus hkz ; ivan-stan .

Nuclear effect calculations in the CSM have been performed so far in momentum representation csm1 ; csm2 ; csm3 ; csm4 ; csm5 , which makes them hardly possible, and in fact this is the reason why the dipole representation for high-energy interactions was first proposed in zkl , and extensively used, in particular for charmonium production off nuclei hikt2 ; kth-psi ; kt-hf ; nontrivial ; rhic-lhc ; kpss-psi . On the other hand, multiple interactions in a nucleus factorize in impact parameter representation, which is then the most appropriate for calculation of the nuclear effects.

III.1.1 CSM via dipoles: the size distribution

First of all, one should formulate the CSM in terms of dipole interactions. As is explained in detail in Ref. kt-hf , the cross section of the process $g+p\to\bar{c}cg+X$ is given by the cross section of of the 4-body dipole $|gg\bar{c}c\rangle$, $\sigma_{4}(\vec{r},\vec{\rho},\alpha,\alpha_{g})$, where $\vec{r}$ is the $c$-$\bar{c}$ transverse separation; $\vec{\rho}$ is the transverse distance between the center of gravity of the $\bar{c}c$ and the radiated gluon. The second gluon in the 4-body dipole is the time inverted initial gluon, whose transverse position coincides with the center of gravity of the whole system. The fraction of the light-cone momentum of the initial gluon, carried by the final gluon is $\alpha_{g}$; and the fractional momenta of $c$ and $\bar{c}$ inside the produced colorless dipole, projected to the $J/\psi$ wave function, are $\alpha$ and $\bar{\alpha}=1-\alpha$ respectively.

Notice that that the mean values of $\langle r^{2}\rangle$ and $\langle\rho^{2}\rangle$ are controlled by different mass scales. While the former is related to the heavy quark mass, $r\sim 1/m_{c}$, the latter is controlled by a semi-hard scale, related to the nonperturbative dynamics. It has been determined by phenomenological analysis of data kst2 ; spots , with fixed $m_{g}^{2}\approx 0.5\,\mbox{GeV}^{2}$, which can be treated as an effective gluon mass squared. The calculations are significantly simplified, if the small $\langle r^{2}\rangle$ is neglected compared with $\langle\rho^{2}\rangle$. Then the color-octet $\bar{c}c$ pair can be treated as point-like, i.e. is equivalent to a gluon, so the $\sigma_{4}$ takes the form of a 3-gluon dipole cross section kt-hf ,

[TABLE]

where

[TABLE]

We remind that all dipole cross sections depend also implicitly on $x_{2}$, related to the rapidity $y$ of the produced $\bar{c}c$,

[TABLE]

where $M_{\bar{c}c}$, $p_{T}$ and $y_{\bar{c}c}$ are the mass, transverse momentum and rapidity of the $\bar{c}c$ pair produced in glue-glue fusion, with subsequent decays to $J/\psi$ and gluon. For the sake of simplicity we will associate them with the mass and rapidity of the detected $J/\psi$, unless specified otherwise.

For further calculations we need to make a choice of parametrization of the dipole cross section $\sigma_{\bar{q}q}(r,x_{2})$. Hereafter we rely on the parametrization bgbk fitted to HERA data 222More recent analyses, which also include impact parameter dependence of the elastic dipole amplitude are now available amir . For our puposes a $b$-integrated cross section is sufficient,

[TABLE]

where the parameter $\sigma_{0}$ and the scale $\mu^{2}$ are defined in bgbk .

The presence of the gluon density in (21) shows that this parametrization corresponds to the Pomeron contribution to the dipole cross section. This is the reason why it describes well the DIS data only at sufficiently small $x_{2}<0.01$ bgbk . At larger $x_{2}$ the Reggeon contribution, which corresponds to valence quarks in $F_{2}(x,Q^{2})$, increases, and the Pomeron alone fails to describe data. This problem, however, is relevant only for light quarks, which dominate in the $F_{2}(x_{2},Q^{2})$ measured at HERA. For $\bar{c}c$ dipoles the Reggeon term, corresponding to valence $\bar{q}q$ exchanges, is suppressed by the OZI rule ozi1 ; ozi2 ; ozi3 , which suppresses valence charm component in the proton. Smallness of such a component (intrinsic charm stan1 ; stan2 ) is confirmed by data pumplin , so it can be neglected.

For $\bar{c}c$ dipoles the Reggeon term is suppressed by the OZI rule and can be neglected. ozi1 ; ozi2 ; ozi3 Thus, the parametrization Eq. (21) for $\bar{c}c$ dipoles can be safely extended up to $x_{2}\sim 0.1$, where the coherence length Eq. (1) shrinks down to the nucleon size.

The cross section of $J/\psi$ production is derived in Appendix B and is given by Eq. (B.16). Since the amplitude contains the projection to the $J/\psi$ wave function, the cross section contains integrations over $\vec{r}$ and $\vec{r}^{\,\prime}$. On the other hand, the radiated gluon is not registered, and integration over its transverse momentum produces a delta function $\delta(\vec{\rho}-\vec{\rho}\,^{\prime})$. Therefore the size distribution function $W(\vec{\rho},\vec{r},\vec{r}^{\,\prime})$ depends on only three variables. We normalize this function to unity, and relate it to the $pp$ differential cross section of charmonium production, presented in Eqs. (B.2)-(B.16),

[TABLE]

This distribution also depends implicitly on $x_{2}$.

III.2 Nuclear effects

Now we are in a position to predict the nuclear effects,

[TABLE]

Here $\vec{B}$ is the impact parameter of the $pA$ collision; $z$ is the longitudinal coordinate of the incoherent color-exchange interaction, which leads to the production of a colorless $S$-wave $\bar{c}c$ dipole, projected to the $J/\psi$ wave function. The nuclear suppression factor $S_{A}^{(1N)}$ includes shadowing due to reduction of the $\bar{c}cg$ flux at $z^{\prime}<z$ and attenuation of the produced colorless $\bar{c}c$ dipole at $z^{\prime}>z$,

[TABLE]

where $\sigma_{4}(\rho,\alpha_{g})$ is given by Eq. (18), and $\Sigma_{1}(\vec{r},\vec{r}^{\,\prime})=[\sigma_{\bar{q}q}(r)+\sigma_{\bar{q}q}(r^{\prime})]/2$, by Eq. (A.11). The nuclear thickness functions, $T_{-}(B,z)$ and $T_{+}(B,z)$, which correspond to the propagation of the projectile $g\bar{c}c$ fluctuation up to the point $(\vec{B},z)$ and the propagation of the produced $\bar{c}c$ dipole afterwards, respectively.

[TABLE]

Apparently, $T_{-}(B,z)+T_{+}(B,z)=T_{A}(B)$, the full thickness function given by (16).

Now we can calculate the single-step term Eq. (23) and the results at $\sqrt{s}=200$ and $5000\,\mbox{GeV}$ are plotted in Figs. 4 and 5 by blue curves labelled as $R^{1N}$.

These results are close to the first simplified calculations done in nontrivial ; rhic-lhc , which agreed reasonably well with data phenix-psi at $\sqrt{s}=200\,\mbox{GeV}$, but grossly under-predicted the ratio $R_{pA}$ at the energy of LHC alice1 ; alice2 . This fact was already highlighted in andronic .

Notice that such a contradiction with the observed energy dependence of the nuclear ratio is not a simple failure of a concrete model, but discloses a deeper puzzle. The dipole cross section is well constrained by precise DIS data from HERA. It is known to steeply rise with $1/x$, therefore the magnitude of nuclear attenuation of dipoles must rise with energy. This expectation is beyond the details of a particular model, and cannot be easily changed. The observed similarity of nuclear suppression at both RHIC and LHC energies should be treated as an indication of a new mechanism of $J/\psi$ production in nuclei, for which a natural candidate is the second term in Eq. (17).

IV Double-step production

The second term in Eq. (17) is given by,

[TABLE]

where the double-step contribution to the numerator is illustrated in Fig. 2. Summing over final states one arrives at the cross section, expressed in terms of the density matrix, as is described in Appendix A.

The first color-exchange interaction, $g+N\to\bar{c}c+X$, can result in the production of a $\bar{c}c$ pair in three different states (at leading order): (i) antisymmetric relative to permutations of space and spin variables, color singlet $\{1^{-}\}$ or color octet $\{8^{-}\}$states; (ii) symmetric in spacial-spin variables, color octet state $\{8^{+}\}$. The notations used here are from kt-hf

Assuming that the finally produced state after the second interaction is a colorless $S$-wave $\bar{c}c$ dipole {$1^{+}$}, the intermediate $\bar{c}c$ pair, between the first and second collisions, must be a $P$-wave $\{8^{-}\}$ state. The first collision cross section has the form kth-psi ; kt-hf ,

[TABLE]

where $\Phi_{g}^{\mu\bar{\mu}}(\vec{r},\alpha)$ is the light-cone distribution function of $\bar{c}c$ incoming gluon, defined in (B.6);

[TABLE]

We fixed here $\alpha=\alpha^{\prime}=1/2$, because these values are strongly enhanced by the projection into the charmonium wave function kth-psi ; klss .

The second interaction, $\{\bar{c}c\}_{\{8^{-}\}}N\to J/\psi X$, is the time reversal of the usual inelastic (color exchange) interaction, $J/\psi N\to X$, which is related to the dipole cross section,

[TABLE]

Thus, we are in a position to calculate the numerator of the double-scattering term Eq. (26) as,

[TABLE]

where the gluon PDF in the beam proton, $g_{N}(x_{1})$, is taken at the scale $Q^{2}=4m_{c}^{2}$. The wave function of quarkonium is normalized according to

[TABLE]

For evaluations, we rely on the LC charmonium wave function obtained with the Cornell potential cornell1 ; cornell2 and boosted to another frame following the procedure developed in klps .

The $\bar{c}c$ light-cone distribution function is convoluted in (30) including the Clebsch-Gordan coefficient $\left\langle 1M\left|\frac{1}{2}\bar{\mu}\frac{1}{2}\mu\right.\right\rangle$, and $M$ is the spin $z$-projection. The nuclear suppression factor $S_{A}^{(2N)}$ is presented below.

IV.1 The nuclear suppression factor

This factor gets contributions from different parts of the dipole path through the nucleus (see Fig. 2) : (i) prior the first collision at longitudinal coordinate $z_{1}$ and production of the color-octet, P-wave $\bar{c}c$ pair $\{8^{-}\}$; (ii) attenuation of the produced $\{\bar{c}c\}_{8^{-}}$ pair on the path from $z_{1}$ up to the next color-exchange interaction at $z_{2}$; (iii) attenuation of the produced colorless dipole $\{\bar{c}c\}_{1^{+}}$on its way out of the nucleus. Correspondingly, the nuclear suppression can be presented as a product of three factors,

[TABLE]

The first factor $S_{1}(z<z_{1})$ has the meaning of shadowing, namely the competing probabilities of the process $g\to\bar{c}c$ to occur on different bound nucleons, which reduce the gluon flux piller ; kt-hf ,

[TABLE]

where $\Sigma_{3}=\left[\sigma_{3}(r,\alpha)+\sigma_{3}(r,^{\prime}\alpha^{\prime})\right]/2$, and

[TABLE]

The cross section $\sigma_{3}(r,\alpha)$ controlling the suppression, is the total cross section of a 3-body dipole ($g\bar{c}c$), responsible for the inclusive production process $gN\to\bar{c}cX$ piller ; kt-hf .

The second factor in (32) can be treated as the survival probability of the produced $(\bar{c}c)_{\{8^{-}\}}$ pair propagating through the medium. Its attenuation is controlled by only a part of the cross section $\Sigma_{8}(\vec{r},\vec{r}^{\,\prime})$ introduced in (A.11). While the diagonal transitions $\{\bar{c}c\}_{\{8^{-}\}}\to\{\bar{c}c\}_{\{8^{-}\}}$ do not affect the final result, the other channels, such as transitions of $\{\bar{c}c\}_{\{8^{-}\}}$ to a singlet $\{\bar{c}c\}_{1^{+}}$, or to a color octet S-wave $\{\bar{c}c\}_{\{8^{+}\}}$, eliminate further possibilities of production of $J/\psi$ at $z=z_{2}$. Summing up the cross sections of the last two channels, we arrive at the second suppression factor in (32),

[TABLE]

where $T_{12}(B,z_{1},z_{2})=T_{-}(B,z_{2})-T_{-}(B,z_{1})$, and

[TABLE]

In the last line we again employ the approximation $\alpha=\alpha^{\prime}=1/2$, for the sake of simplicity.

The last factor in (III.2) has a rather obvious form,

[TABLE]

where $\Sigma_{1}(\vec{r},\vec{r}^{\,\prime})$, is given by Eq. (A.11).

Notice that the $z$-dependent part of (23) can be integrated analytically,

[TABLE]

where we introduced the short-hand notations, $\Omega_{1}=\Sigma_{3}-\Sigma_{\{8^{-}\}}$; $\Omega_{2}=\Sigma_{\{8^{-}\}}-\Sigma_{1}$; $\Omega_{3}=\Sigma_{3}-\Sigma_{1}$.

IV.2 The $pp$ reference

In our calculation of $R^{(1N)}_{pA}$ for the single-step mechanism, we assumed that the same CSM model dominates both the numerator and denominator, and therefore they have nearly identical functional forms, except for the nuclear suppression factor and some corrections discussed below. So the reference $pp$ cross section nearly cancels.

The double-step term $R^{(2N)}_{pA}$ evaluation is more peculiar, because the numerator and denominator originate from different mechanisms and have distinct functional forms. While the former, given by Eq. (23), is calculated directly based on the well developed dipole phenomenology, the latter depends on the choice of a model for inclusive $J/\psi$ production (see Sect. III) and is assumed here to be dominated by CSM. Thus, the denominator of $R^{(2N)}_{pA}$ Eq. (26) has a rather wide theoretical uncertainty band, which is related to the accuracy of the CSM, and possibility of other missed contributions (such as three-gluon fusion motyka , certainly important at very forward/backward rapidities).

The least model dependent way to treat the denominator of (26) would be to take it directly from a fit to experimental data for $pp\to J/\psi X$, available within certain kinematic domains. We rely on our evaluations of the $pp\to J/\psi X$ cross section, performed within the dipole version of the CSM in Appendix B. The results, compared with data in Figs. 21 and 22, well reproduce the shape of the $y$-dependence of the cross section, however, slightly underestimate the normalization. At $\sqrt{s}=200\,\mbox{GeV}$ we employed the data from phenix-psi ; phenix-pp . Lacking experimental results at $\sqrt{s}=5\,\mbox{TeV}$, we interpolated between data at $\sqrt{s}=2.76\,\mbox{TeV}$ and $7\,\mbox{TeV}$ interpolation . The details are presented in Appendix B and the results are depicted in Fig. 22. Since, as we said, data is the most reliable source of information about the $pp$ cross section, we adjusted the normalization of the theoretical curves to fit the data, keeping the shape of the $y$-dependence unchanged.

Now we are in a position to calculate $R^{2N}$, the ratio of the cross section Eq. (30) to the chosen $pp$ reference, and the results at $\sqrt{s}=200$ and $5000\,\mbox{GeV}$ are plotted in Figs. 7 and 8 by dashed curves labelled as $R^{2N}$.

V Gluon shadowing

Leading twist gluon shadowing originates in the nuclear rest frame from coherent multiple interactions of the radiated gluons. It can also be treated as the contribution of higher Fock components in the projectile hadron, containing extra gluons, which have a coherence (radiation) time sufficiently long to experience multiple interactions in the nucleus kst2 . These gluons are complementary to the gluon radiated within the CSM mechanism (Fig. 3). Unlike quark shadowing, which is known to onset at $x_{2}\lesssim 0.1$ krt-shad1 , gluon shadowing needs an order of magnitude smaller $x_{2}$ to show up krt-shad2 . This is controlled by the coherence length of gluon radiation,

[TABLE]

which must be longer than the mean free path in nuclear matter. The factor $P_{g}\approx 0.1$, evaluated in krt-shad2 , makes the coherence time of gluon radiation significantly shorter than the Ioffe time for quarks. This happens due to the enhanced transverse momenta of gluons in hadrons kst2 ; spots , which make the fluctuations containing gluons much heavier. For the same reason, the mean quark-gluon separation is short, and the magnitude of the leading-twist gluon shadowing turns out to be rather small, even compared with the higher-twist quark shadowing. The weaknees of gluon shadowing, predicted in kst2 , was confirmed by the NLO analysis of DIS data florian1 ; florian2 .

The gluon shadowing suppression factor $R_{g}(x,Q^{2})$, calculated in kst2 , was applied to Drell-Yan process in jkrt-dy , and to heavy flavor production in kt-hf , where one can find the details of the calculations. This factor suppresses $J/\psi$ production on nuclei as well. In our case we include gluon shadowing by reducing the dipole cross section with the shadowing factor $R_{g}$, which also depends on the nuclear impact parameter $b$. Such a way of incorporation of gluon shadowing can be justified only at first order, which corresponds to radiation of a single gluon. In fact, radiation of two gluons lead to a quadratically short coherence time compared with Eq. (39) k-poetic , too short to cause shadowing at currently available energies.

The terms $R^{1N}$ and $R^{2N}$ in (17), with added gluon shadowing corrections at $\sqrt{s}=200\,\mbox{GeV}$ and $5000\,\mbox{GeV}$, are depicted in Figs. 4 and 5 respectively. The corrections are found to be rather small at the energy of RHIC (due to shortness of the coherence length), but significant at LHC. Nevertheless, even at the LHC energy gluon shadowing vanishes in the backward hemisphere, towards the minimal rapidity $y\sim-4$ in the kinematical range measured so far, because the coherence length Eq. (39) becomes shorter than the mean spacing between bound nucleons.

VI Energy loss

VI.1 Nonperturbative energy loss

Apparently, multiple soft interactions in the nuclear medium should lead to dissipation of energy by the projectile partons, reducing the production rate of $J/\psi$ at large Feynman $x_{F}$, where the restricted phase space of produced $J/\psi$ becomes an issue. Energy loss was first proposed in kn-psi as a mechanism of suppression of the $pA$-to-$pp$ ratio of $J/\psi$ production at large $x_{F}$, observed in serpukhov ; na3 . The rate of energy loss, treated within the string model, was independent of the incoming proton energy kn-psi . Perturbative calculations, performed in the approximation of soft gluon radiation, confirmed the string model result of energy independent parton energy loss feri ; bh . This, however, could not explain the observed $x_{F}$ scaling, i.e. similarity of the $x_{F}$-dependences of nuclear effects in $J/\psi$ production at different energies serpukhov ; na3 ; e866 .

Nonetheless, later, in e-loss ; boosting , it was found that the rate of energy loss, either in nonperturbative e-loss , or perturbative boosting regimes, rises proportional to the incoming energy. This is easily interpreted in terms of Fock-state representation for the light-cone wave function of the incoming hadron. The probability of giving a significant fraction of the hadron momentum to one parton (soft or hard) is more suppressed in the higher Fock states. Indeed, if one of the participating partons gets a large momentum fraction $x_{1}\to 1$, all other participants are pushed into a small phase space with $x<1-x_{1}$. The measured parton distribution function (PDF) is averaged over different Fock components, and the interaction of these Fock states with the nuclear target changes their weights, increasing the contribution of higher Fock components, so that the projectile parton distribution becomes softer, i.e. more suppressed at large $x_{1}\to 1$. Thus, the projectile proton PDF becomes target-dependent, violating QCD factorization at large $x_{1}$, where the energy sharing (energy loss) problem becomes important e-loss ; boosting . Such a beam-target correlation breaks factorization, because it occurs at a low scale. This explains why every process measured so far was found to be nuclear suppressed at large $x_{1}$ e-loss .

Glauber multiple hadron-nucleus soft inelastic interactions are not sequential (as is frequently naively believed), but correspond to multi-sheet configurations in the topological $1/N_{c}$ expansion of QCD for the inelastic amplitude, i.e. they are related to simultaneous propagation and interaction in the medium of different projectile partons from a high Fock component of the incoming hadron gribov1 ; gribov2 ; agk ; kaidalov ; capella . This leads to the problem of energy sharing between participating partons, which becomes especially severe at large fractional momentum $x_{1}$ carried by one of them. The associated nuclear suppression was calculated in e-loss using the Fock state expansion, weighted by the interaction with the target, corresponding to the Glauber model. The suppression factor $S(x_{1})$ for each additional topological sheet was evaluated in kaidalov ; capella relying on Regge phenomenology, and in e-loss treating it as a rapidity gap survival probability. Both approaches led to the same result: at $x_{1}\to 1$ suppression increases as $S\propto(1-x_{1})$. We apply here the model for energy loss developed in e-loss , in order to correct the nuclear ratio (17).

VI.2 Perturbative energy loss

Another source of nuclear modification of the projectile gluon distribution is an increased hard scale. Indeed, if in $pp$ collision the gluon distribution is taken at the scale $Q^{2}=4m_{c}^{2}$, a nuclear target generates another scale, known as saturation scale $Q_{s}^{2}$. So the effective scale of the process increases, $Q_{eff}^{2}=4m_{c}^{2}+Q_{s}^{2}$. This follows naturally from the interpretation of saturation in the rest frame of the nucleus, which is related to broadening of the transverse momentum of a gluon propagating through the nucleus broadening ,

[TABLE]

We employ the dipole description of broadening jkt , and for the saturation scale rely on the results of broadening . This result is based on the approximation of Bethe-Heitler regime of gluon radiation in multiple interactions, neglecting interferences of gluons radiated in collisions with different nucleons. Effects of coherence cause deviations from Eq. (40), however according to the discussion in Sect. V effects of coherence in gluon radiation are small even at the energies of LHC. Therefore, in what follows we employ the approximate effective scale $Q_{eff}^{2}=4m_{c}^{2}+Q_{s}^{2}$ for numerical evaluations.

Notice that broadening of the transverse momentum of a gluon propagating through the nucleus is equivalent to the effect of saturation in the $k_{T}$-dependent PDF of the nucleus in its infinite-momentum frame saturation .

With a larger scale the process resolves more partons in the incoming proton. Thus, via the effect of broadening the nuclear target activates higher Fock states in the incoming proton. The result is qualitatively similar to what we observed above, namely, parton density will be enhanced at small $x_{1}$, but suppressed at $x_{1}\to 1$. Such a nuclear modification of the gluon density in the incoming proton can be performed by evolving the projectile proton PDFs with DGLAP equations from the scale $Q^{2}=4m_{c}^{2}$ to $Q^{2}+Q_{s}^{2}$. Then the gluon PDF in the proton should be replaced $g_{N}(x_{1})\Rightarrow\tilde{g}_{N}(x_{1},B)$ in the numerators of $R^{(1N})$ and $R^{(2N})$. Some examples of modifications, $\tilde{g}_{N}(x_{1},B)/g_{N}(x_{1})$ are shown in Fig. 6 for $pPb$ collisions at $\sqrt{s}=5\,\mbox{TeV}$ vs $x_{1}$ and impact parameter $B$.

This modification of the $x_{1}$ dependence of the projectile gluon distribution can be treated as an effective energy loss, leading to nuclear suppression of heavy quark production at forward rapidities (large $x_{1}$). The results presented in Fig. 6 show that the effect is extremely weak, only a few percent suppression at very forward rapidities. The reason for this weakness can be easily traced in Fig. 1 of Ref. boosting . One can see that the effect of induced energy loss is controlled by the relation between the scale of the process, $Q^{2}$, and the saturation scale $Q_{s}^{2}$. The effect may be strong if $Q_{s}^{2}\gg Q^{2}$, but becomes vanishingly small at $Q^{2}\gg Q_{s}^{2}$.

Intuitively, this is pretty clear. It can be interpreted as a vacuum dead-cone effect hf-jets , namely a parton originating from a hard process at scale $Q^{2}$ is lacking gluon field with small transverse momenta $k_{T}^{2}<Q^{2}$. Gluon bremsstrahlung and medium induced energy loss of such a parton are significantly reduced compared to a nearly on-mass-shell parton. This is what we see in the above Fig. 6, where the characteristic scale of the process, $Q^{2}\approx 10\,\mbox{GeV}^{2}$, exceeds considerably the saturation scale.

Reduction of induced energy loss by a large genuine scale $Q^{2}$ of the process can be also interpreted in terms of the Landau-Pomeranchuk effect, which says that on a long length scale $l\gg R_{A}$, the radiation spectrum depends on the total accumulated kick acquired by the charge, rather than on the details of several kicks occurring on a short length scale (the nuclear radius $R_{A}$). The radiation spectrum $dk_{T}^{2}/k_{T}^{2}$ leads to a logarithmic scale dependence of the radiated energy. The induced energy loss is given by a difference between energies radiated in the processes with the effective scales $Q^{2}+Q_{s}^{2}$ (in pA) and $Q^{2}$ (in pp). Thus, the induced energy loss exposes the following scale dependence,

[TABLE]

if $Q^{2}\gg Q_{s}^{2}$, i.e. it turns out to be suppressed. This effect is of course included in the DGLAP analysis, whose results are presented in Fig. 6.

Notice that the suppressing effect of a large scale of the process was missed in the calculations ap of induced energy loss in charmonium production. As a result, the magnitude of energy loss was grossly overestimated compared with the DGLAP analysis.

VI.3 Numerical results for $J/\psi$

Now we are in a position to finalize the calculations of nuclear effects in $J/\psi$ production. The effects of energy loss, or modification of the projectile gluon distributions, have been already incorporated into our previous results corrected for gluon shadowing, as was plotted by the dashed curves in Fig. 4 and Fig. 5. The final results are compared with available data at $\sqrt{s}=200\,\mbox{GeV}$ in Fig. 7 and at $\sqrt{s}=5000\,\mbox{GeV}$ Fig. 8.

As was anticipated, the energy loss effects are strongest at the energies of RHIC. A substantial modification of nuclear effects due to energy loss has been already observed for other hard processes in high-pt ; e-conserv . Our results seem to agree reasonably well with data, especially taking into account the large uncertainties in the $pp$ reference, affecting the term $R^{(2N)}$ in (17).

In view of the forthcoming LHC measurements of $pA$ collisions at $\sqrt{s}=8000\,\mbox{GeV}$, we notice that our predictions are hardly different from those presented in Fig. 8 for $\sqrt{s}=5000\,\mbox{GeV}$.

VI.4 Nuclear modification of the $p_{T}$-distribution

Multiple interactions of the projectile partons in the nucleus are known to lead to broadening of the transverse momentum, the phenomenon also called saturation or color glass condensate. It can be effectively evaluated within the dipole phenomenology jkt , well adjusted to HERA data on small-$x$ DIS. The value of broadening at impact parameter $B$ is given by Eq. (40) derived in jkt .

Nuclear broadening of the $p_{T}$-distribution naturally leads to a ratio $R_{pA}(p_{T})$, rising with $p_{T}$, the effect, usually named after Cronin. The $p_{T}$ dependence of the $J/\psi$ production cross section in $pp$, $pA$ and $AA$ collisions is well described by the form, $d\sigma/dp_{T}^{2}\propto(1+p_{T}^{2}/6\langle p_{T}^{2}\rangle)^{6}$ e789 ; phenix-pt ; alice2 . Therefore, making a shift of $\langle p_{T}^{2}\rangle$ for $pA$ in comparison with $pp$ collisions, one arrives at a $p_{T}$-dependent nuclear ratio kps-psiAA ,

[TABLE]

where $R_{pA}$ in the r.h.s. of (42) is the ratio of the $p_{T}$-integrated cross sections (as was calculated above); $\xi=1+\Delta_{pA}(x_{2})/\langle p_{T}^{2}\rangle$; and $\Delta_{pA}(x_{2})=\langle p_{T}^{2}\rangle_{pA}-\langle p_{T}^{2}\rangle_{pp}$ is nuclear broadening of charmonium transverse momentum.

The magnitude of broadening was evaluated in broadening . At $\sqrt{s}=5.02\,\mbox{TeV}$ and the rapidity intervals of interest, $y\in(-4.46,-2.96)$, $y\in(-1.37,-0.43)$ and $y\in(2.03,3.53)$ the broadening magnitudes, averaged over impact parameters, are $0.35$, $0.73$ and $2.27\,\mbox{GeV}^{2}$ respectively. The $p_{T}$-dependent $R_{pA}(p_{T})$, given by Eq. (42), calculated with these values and $\langle p_{T}^{2}\rangle=7\,\mbox{GeV}^{2}$ alice2 are compared with data in figures 9 - 11, demonstrating good agreement.

VII Production of $\psi(2s)$

The first radial excitation $\psi(2S)$ has the mean radius squared about twice as large as that of $J/\psi$ buchmueller ; cornell1 ; cornell2 , and therefore comparison of nuclear effects for these two charmonium states offers a sensitive test of the production dynamics. Expectations are usually based on either of two popular ideas, both incorrect:

(i) The effect of color transparency makes the nuclear medium more transparent for smaller size state, $J/\psi$, which is expected to be considerably less suppressed than $\psi(2S)$. However, experiments at the SPS na3 and Fermilab e866 found similar magnitudes of nuclear suppression for the two charmonium states.

(ii) At first glance, the observed similarity of nuclear effects can be understood in line with the hierarchy of characteristic length scales discussed in Sect. II.1. Indeed, at high energies the formation length Eq. (2) substantially exceeds the nuclear dimension, so a perturbatively small $\bar{c}c$ dipole, rather than a formed charmonium of much larger larger size, propagates through the nucleus. Then one expects the dipole to evolve into either $J/\psi$ or $\psi(2S)$ outside of the nucleus, after experiencing an universal nuclear attenuation on the early perturbative stage. Naively, one might expect universal nuclear suppression for different charmonia. However, the dynamics, controlling the nuclear effects is more involved.

The second proposal (ii) explains why the first one, (i), is incorrect. Nonetheless, an universal nuclear attenuation of a $\bar{c}c$ dipole does not lead to an universal charmonium suppression, because the projection of the produced $\bar{c}c$ distribution function to the charmonium wave function depends on the latter. In particular, spectacular effects are expected for production of $\psi(2S)$, related to the specific shape of its wave function, which has a node and changes sign as function of the $\bar{c}c$ separation.

Unusual features of $\psi(2S)$ production were revealed in photoproduction of charmonia kz91 , the process of a similar, although simpler dynamics compared with hadro-production. It was found that in spite of its large size, the $\psi(2S)$ produced in nuclei may be less suppressed compared with $J/\psi$, sometimes even enhanced. This can be interpreted either in terms of the multi-channel generalised Glauber model hufner96 , or within the dipole description as a result of the specific nodal structure of the $\psi(2S)$ wave function kz91 ; knnz93 ; hikt2 . The $\bar{c}c$ distribution function, to be projected to the charmonium wave function, has a rather wide $r$-distribution, which peaks at $r\sim 2/m_{c}$ nontrivial ; knnz93 , close to the node position in the $\psi(2S)$ wave function. Therefore, a part of the overlap integral extends beyond the node and contributes with a negative sign, causing a significant compensation between dipole separations smaller and larger than the node position. This cancellation contributes to the observed suppression of $\psi(2S)$ production knnz93 in $pp$ collisions. A nuclear target serves as a color filter, which removes the large-size $\bar{c}c$ dipoles, and therefore the mean size of the $\bar{c}c$ wave packet is reduced and the overlap with the $\psi(2S)$ wave function increases.

The results of calculations of the nuclear ratio $R_{pA}$ for $\psi(2S)$ are compared with available data at RHIC and LHC in Figs. 12 and 13, respectively.

The double scattering term $R^{(2N)}$ turns out to be very small for $\psi(2S)$ at the energies of RHIC, but rises to a sizeable corrections at higher energies.

Again, we can conclude that our calculations do not contradict data, which has rather large errors. However, our results for the double ratio $R_{pA}^{\psi(2S)}/R_{pA}^{J/\psi}$, plotted in Fig. 14, show rather small values slowly rising with energy. These results contradict the precise data of E866 experiment e866 , which show that at small $x_{F}$ the double ratio is about $R_{pA}^{\psi(2S)}/R_{pA}^{J/\psi}=0.9$, with a small error.

The nuclear effects observed for the production of the first radial excitation $\psi(2s)$ demonstrate suppression, similar to $J/\psi$, in the energy range of fixed target experiments na3 ; e866 . However, in the energy range of RHIC-LHC, a stronger suppression of $\psi(2s)$ relative to $J/\psi$ was observed phenix-2s ; alice-2s .

VIII Upsilon production

The developed dipole description of charmonium production in $pA$ collisions can be naturally extended for bottomium production, replacing the charm quark mass by $m_{b}=4.5\,\mbox{GeV}$. In the Figures 15 and 16 we present the results at the energies of RHIC and LHC respectively. The term $R^{(1)}$ closely reproduces the earlier calculations in nontrivial , except for the added energy loss effect, which affects the results for RHIC, but not for LHC.

Due to larger $b$-quark mass and smaller dipole sizes, the two-nucleon term $R^{(2N)}$ in (17) is relatively smaller compared with $J/\psi$ production, as one can see in Figs. 15 and 16. As for charmonium, we calculate the $pp$ reference cross section used in the denominator of $R^{(2N)}$, within the CSM, and adjust its normalization to data.

The only available data point phenix-upsilon , plotted in fig. 15 has too large error bar to be considered as a support for our calculations.

We also performed calculations for the double ratio $R_{pA}(\Upsilon(2S))/R_{pA}(\Upsilon(1S))$, and plotted it as function of $\bar{b}b$ energy in Fig. 14. This ratio was measured with a good precision in the CMS experiment at $\sqrt{s}=5000GeV$ and $|y|<1.93$ cms-upsilon-double . This point, plotted in Fig. 14 at energy $E_{\bar{b}b}=e^{y}M_{\Upsilon}\sqrt{s}/2m_{N}$, agree well with our parameter-free calculations.

IX Summary and conclusions

The main objective of this work was to settle the problem of the energy independence of nuclear effects for $J/\psi$ production, observed in $pA$ collisions. This independence of energy is in striking contradiction with the steep energy dependence of the dipole cross section observed at HERA , which controls the nuclear effects. We revealed a novel mechanism enhancing charmonium production at high energies, which comes from the next order of the opacity expansion.

Crucial for the results was the choice of mechanism dominating the production of heavy flavor vector mesons in pp collisions. We favored the color-singlet model (CSM), which can dominate the small-$pT$ quarkonium production we are interested in. We developed a color-dipole formulation of CSM, which is crucial for the calculation of nuclear effects.

The second order term in the opacity expansion for the production cross section is dominated by a different mechanism, a double color-exchange interaction of the projectile heavy $\bar{Q}Q$ dipole. Its contribution helps to reach agreement with data for the nuclear suppression of $J/\psi$ production both at the energies of RHIC and LHC.

Other nuclear effects, gluon shadowing and energy loss, have also been included. Gluon shadowing corrections are found to be important at the energies of LHC, but very small at RHIC. On the contrary, energy loss effects substantially suppress quarkonium production rates at forward rapidities at RHIC, while have no influence at the energies of LHC. The main contribution comes from the nonperturbative mechanism of energy loss, related to the energy sharing problem at forward rapidities. The perturbative energy loss generated by $p_{T}$ broadening was found to be suppressed by the smallness of the saturated momentum relative the scale of the process. This suppression was missed in the previous calculations of the energy loss effect, which has been grossly overestimated.

Although we restricted these calculations with the $p_{T}$-integrated cross sections, the $p_{T}$-dependent ratio $R_{pA}(p_{T}$ was also evaluated, based on the known empirical shape of the $p_{T}$-distribution and the value of broadening, calculated in a parameter-free way (although not free of assumptions) within the dipole phenomenology. The results, obtained for several rapidity intervals, well agree with ALICE data.

Production of radial excitations, vector quarkonia in the $2S$ state, has always attracted interest, related to the nodal structure of the wave function. Differently from photoproduction, where $2S$ states are enhanced compared with the ground state, in hadroproduction we found a strong nuclear suppression of the $\psi(2S)$ to $J/\psi$ ratio, in good agreement with data. At the same time, for bottomia, the $2S$ to $1S$ ratio is nearly unaffected by the nuclear effects, what could be anticipated, because the $\bar{b}b$ dipoles are much smaller compared with $\bar{c}c$, so the convolution with the $\Upsilon$ wave function is less important.

Acknowledgements.

This work was supported in part by Fondecyt (Chile) grants 1170319, 1140842, 1140390 and 1140377, by Proyecto Basal FB 0821 (Chile), and by CONICYT grant PIA ACT1406 (Chile) . Powered@NLHPC: This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02). Also, we thank Yuri Ivanov for technical support of the USM HPC cluster where a part of evaluations has been done.

Appendix A Multiple color-exchange interactions of a high-energy dipole

At sufficiently high energy, when the length scales discussed in Sect. II.1 considerably exceed the nuclear dimensions, one can treat the transverse size of such a dipole as ”frozen” by Lorentz time dilation during propagation through the nucleus. The kinematic constraints for this regime can be found in Sect. II.2. This is a perturbative stage of interaction, so the one-gluon approximation for dipole-nucleon interaction is justified. However, multigluon-exchange interactions with different nucleons are enhances by powers of $A^{1/3}$ and cannot be neglected.

A.1 Evolution of the $\bar{c}c$ density matrix

Multiple soft color-exchange interactions of with the bound nucleons keep the dipole transverse separation $\vec{r}$ unchanged, but destroy the target,

[TABLE]

as is illustrated in Fig. 2. One cannot describe the dipole evolution in terms of the dipole-nucleus amplitude, because in the cross section the final states of each color-exchange collision must be summed-up, as is illustrated in Fig. 17.

Therefore, the dipole propagation in the medium is described in terms of density matrix ${}^{k}_{l}U^{i}_{j}(x_{1},x_{2};x_{1}^{\prime},x_{2}^{\prime})$, where $x_{1,2}$ and $x_{1,2}^{\prime}$ are the transverse coordinates of the quark and antiquark in the two conjugated amplitudes alesha ; k-zam ; kth-psi , which presented graphically in Fig. 18.

We will follow the evolution of the density matrix along the longitudinal coordinate $z$, which measures the propagation of the system through the nucleus.

Before the $\bar{c}c$ pair enters the nucleus, i.e. at $z\to-\infty$, it is in a pure colorless state, i.e.

[TABLE]

where $\Psi_{in}(r)$ is the distribution function of $\bar{c}c$ in the incoming beam, for instance a $\bar{c}c$ component of a projectile gluon.

At $z\to\infty$ the system leaves the nucleus and the density matrix can be projected directly to the final state wave function,

[TABLE]

Since for every interaction of the $\bar{c}c$ in the medium we sum up over the final states of nucleons, the density matrix is a colourless object, i.e. it is invariant under simultaneous rotations in all colour indices $i,j,k,l$. Therefore it can be conveniently decomposed into the irreducible parts corresponding to singlet and octet states of the pair,

[TABLE]

where $z$ is longitudinal coordinate of the target nucleon; $\vec{r}=\vec{x}_{1}-\vec{x}_{2}$, $\vec{r}^{\,\prime}=\vec{x}_{1}^{\,\prime}-\vec{x}_{2}^{\,\prime}$. We assume here that the impact parameters of the centers of gravity of the dipoles in the two amplitudes coincide, which is correct if the dipole-nucleon interaction radius can be neglected compared with the nuclear radius.

$P_{S}$ and $P_{O}$ in (A.4) are the singlet and octet projection operators,

[TABLE]

such that

[TABLE]

The elements $S(\vec{r}=\vec{r}^{\,\prime})$ and $O(\vec{r}=\vec{r}^{\,\prime})$ are the probabilities to find the quark-antiquark pair in color singlet or octet states respectively.

In the one-gluon-exchange model every interaction with a nucleon results in the change of the density matrix ${}^{k}_{l}U^{i}_{j}$, represented schematically in fig.3.

Explicit calculation of the diagrams gives the following variations of the density matrices as function of $z$,

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

If one is not interested in a particular spacial state of the outgoing $\bar{c}c$ pair and regarding only its colour state (e.g., one does not descriminate between different outgoing colourless states like $J/\psi$, $\eta_{c}$, $\chi$, etc.), only the elements diagonal in the space variables $\vec{x}_{1}=\vec{x}_{1}^{\prime}$ and $\vec{x}_{2}=\vec{x}_{2}^{\prime}$ of the density matrix are relevant. Then for $S(\vec{r};z)$ and $O(\vec{r};z)$, which also implicitly depend on $b$, one gets the following system of linear differential equations

[TABLE]

Here $S(r;z)$ and $O(r;z)$ are interpreted as the probabilities to find the $\bar{c}c$ pair in a color singlet or octet states respectively. Since the total probability is conserved,

[TABLE]

Assuming that the initial state is a pure singlet with distribution function $S_{in}(r)$, and solving Eqs. (A.12)-(A.13) one arrives at,

[TABLE]

Correspondingly, for a color-octet initial state one gets,

[TABLE]

We see that for large number of inelastic collisions of the $\bar{c}c$ dipole333One should not mix up this value with the number of collision usually used for normalization of hard reactions in $pA$ and $AA$ collisions. The latter is controlled by $\sigma_{in}^{NN}$, rather than by the small $\bar{c}c$ dipole cross section., $n^{\bar{c}c}_{coll}=\sigma_{\bar{c}c}(r)T_{A}(b,z)\gg 1$ the probability of production of color-singlet or octet states approach universal values, $1/9$ and $8/9$ respectively, independently of the color structure of the incoming $\bar{c}c$ pair. This could be anticipated, since after multiple rotations in the color space both quark become completely unpolarized in color. All of the possible $9$ ($3\times 3$) color states of the $\bar{c}c$ are produced with equal probabilities, and only one of them is a singlet, while the 8 others are octets.

Appendix B $J/\psi$ production in $pp$ collisions

The production of heavy quarks was described within the dipole approach in kt-hf . In the leading order of pQCD it is described by 15 Feynman graphs (Fig. 8 in kt-hf ). Only six of them, presented here in Fig. 20, contribute to the production of $J/\psi$ and its excitations.

B.1 Soft gluon approximation

The amplitude, corresponding to these graphs was derived in kt-hf , in the approximation of small fractional gluon LC momentum $\alpha_{g}\ll 1$.

[TABLE]

Here $\alpha$ and $\bar{\alpha}=1-\alpha$ are the fractional light-cone momenta of the $\psi$, carried by the charm quark and antiquark, respectively. The relative transverse momentum and separation of $c$ and $\bar{c}$ are $\vec{k}$ and $\vec{r}$ respectively. We employ here the result of kt-hf for the production of a colorless $\bar{c}c$ pair in $S$-wave, but projecting it to the charmonium light-cone wave function, $\Psi_{\psi}\left(\alpha,\vec{r}\right)$, normalized as

[TABLE]

The transverse momentum of $\psi$ as a whole, $p_{\psi}$, is related to the transverse momentum transfer to the target, $\vec{k}_{T}$, and the transverse momentum $\vec{k}_{g}$ of the radiated gluon as,

[TABLE]

Further notations in (B.1) are the transverse distances $\vec{b}$ between the target and the center of gravity of $\psi$-$g$, and $\rho$ between the $\psi$ and radiated gluon.

The light-cone distribution function for a quark, radiating a transversely polarized gluon with fractional momentum $\alpha_{cg}$, was derived in kst1 ,

[TABLE]

where $\vec{\rho}_{cg}$ is the transverse separation between the final gluon and quark, and $\tau^{2}=(1-\alpha_{cg})m_{g}^{2}+\alpha_{cg}^{2}m_{c}^{2}$. Notice that the non-perturbative effects strongly affect this distribution function, leading to a significant reduction of the mean quark-gluon separation. The magnitude of this reduction is constrained by the observed suppression of diffractive gluon radiation kst2 , as well by many other processes spots . Here we rely on the perturbative form Eq. (B.4 of the distribution function, but introduce an effective gluon mass $m_{g}\approx 0.7\,\mbox{GeV}$, which can be treated as a transverse mass of the gluon, which has a transverse motion enhanced by the non-perturbative effects.

The indices $\mu$ and $\bar{\mu}$ in B.4 are quark helicities before and after the gluon emission, and $\xi_{\bar{\mu}}$ and $\xi_{\mu}^{\dagger}$ are the spinors of the initial and final quarks respectively. The operator $\hat{\mathcal{Q}}_{cg}$ has the form kst1

[TABLE]

The light-cone distribution function for the $g\to\bar{c}c$ transition is given by

[TABLE]

where

[TABLE]

and

[TABLE]

The fractional momentum $\beta$ of the $c$ quark emerging from the incoming gluon (see Fig. 20) is different from that in the final state, due to gluon radiation by either $c$, or $\bar{c}$ quarks. Correspondingly, $\beta=\alpha/(1+\alpha_{g})$, or $\beta=(\alpha+\alpha_{g})/(1+\alpha_{g})$, as one can see in (B.1). Gluon radiation also changes the $\bar{c}c$ separation $\vec{R}$, which is different from the final $\vec{r}$, as one can see in the argument of $\Phi_{\bar{c}c}$ in Eq. (B.1). The $\bar{c}c$ distribution function contains proper convolution with a Clebsch-Gordan coefficient $\left\langle 1M\left|\frac{1}{2}\bar{\mu}\frac{1}{2}\mu\right.\right\rangle$, where $M$ is the spin projection.

Following the definitions of kt-hf , the function $\gamma(b)$ in (B.1) corresponds to the Fourier image of the dipole destruction amplitude, which can also be treated as an “elastic” (color-exchange) gluon-nucleon scattering amplitude. It is related to the dipole cross-section as

[TABLE]

B.2 The general case of arbitrary $\alpha_{g}$

A gluon, as a vector particle, is usually radiated at high energies with a small fractional momentum $\alpha_{g}\sim 1/\ln(s)$. However, in the process under consideration, the transition of a $\bar{c}c$ pair from color-octet to singlet states, small $\alpha_{g}$ values are suppressed by color screening, and one should go beyond this approximation, Eq. (B.1), and rely on the general form of the amplitude, where

[TABLE]

The functions under the trace operation are here 2$\times$2 matrices in quark helicity space (helicity indices are dropped). The matrix $\Lambda_{M}$ contains the convolution of spinors with the Clebsch-Gordan coefficients from the wave function,

[TABLE]

where $\sigma_{i}$ are the Pauli matrices in helicity space. The multiplier

[TABLE]

takes into account the ordering of $t_{a}$ matrices and a numerical pre-factor.

The functions $\Phi_{cg}(\tau_{n},\rho_{n})$ and $\Phi_{\bar{c}c}(\epsilon_{n},\,r_{n})$ are defined in (B.4) and (B.6) respectively. The contributions of different graphs depicted in Fig. 20 to the amplitude are summed in Eq. (B.2). The fractional momenta $\alpha_{n}$ and $\beta_{n}$, as well as the transverse separations $\vec{\rho}_{n}$ and $\vec{r}_{n}$, depend on the number of the corresponding graph in Fig. 20.

It is assumed that at least one of the quarks is on-shell. The parameters $\epsilon_{n},\,\tau_{n}$ as well as arguments $r_{n},\,r_{G,n}$ for different diagrams (1-6) in the Figure 20 are given by

[TABLE]

where $\bar{\alpha}_{g}=1-\alpha_{g}$.

The $p_{T}$-integrated differential cross-section of the inclusive charmonium production, which describes the distribution over $\rho$ and $r$, can be expressed in terms of the dipole cross-section,

[TABLE]

where $y$ is the charmonium rapidity, and

[TABLE]

In the difference $\vec{b}_{n}-\vec{b}_{n^{\prime}}^{\,\prime}$ in (B.15) the $b$ dependence cancels, so the dipole cross-section $\sigma_{\bar{q}q}$ in (B.2) is function of $\vec{\rho}$ and $\vec{r}$.

The integrated cross section,

[TABLE]

can be compared directly with data. Comparison with available data from RHIC and LHC, is shown in Figs. 21 and 22.

Although the calculations contain no free parameter adjusted to the data to be explained, there are theoretical uncertainties related to the different approximations that had been used. In particular, while the phenomenological dipole cross section takes into account the effects of gluon saturation, important at small $x_{2}$ in one of the protons, we rely on a single gluon approximation in the projectile gluon distribution, which is justified only at large $x_{1}$. Therefore the dipole description is ”asymmetric”, it is reliable only at sufficiently small $x_{2}$, but large $x_{1}$, and vice versa, i.e. at forward-backward rapidities, and in the central rapidity region is least reliable.