$J/\psi$ polarization in high multiplicity $pp$ and $pA$ collisions: CGC+NRQCD approach

Quarkonium production mechanism in high multiplicity small collision systems has recently been pursued in the color-glass-condensate (CGC) effective theory combined with non-relativistic QCD (NRQCD) factorization, allowing to study initial state interactions. Quarkonium polarization, potentially measured in future experiments, would help elucidate the quarkonium production mechanism at high multiplicities. In this paper, we provide predictions on $J/\psi$ polarization parameters in high multiplicity proton-proton ($pp$) and proton-nucleus ($pA$) collisions within the CGC+NRQCD framework. Theoretical predictions are given for $J/\psi$ rapidity $2.5<y_{J/\psi}<4$, charged-particle multiplicity pseudorapidity $|\eta_{ch} |<1$ and energies $\sqrt{S}=13\mathrm{~TeV}$ for $pp$, $\sqrt{S}=8.16\mathrm{~TeV}$ for $pA$ collisions. Considering two leptonic frame choices (Collins - Soper and helicity) we find a weak polarization of $J/\psi$ that additionally decreases with growing event activities. No significant system size dependence between $pp$ and $pA$ collisions is obtained - this could be a new constraint to initial state interactions in small collision systems.

Within the CGC framework, the cross-section for heavy quark pair production in a dilute-dense system (i.e. pA and forward pp) was computed in Refs. [46][47][48], and short distance coefficients (SDCs) in the CGC framework combined with the NRQCD factorization were calculated in Refs. [49,50]. Thus far, the CGC+NRQCD framework has been successful in describing the p ⊥ spectra of J/ψ in minimum bias pp and pA collisions [51,52] and polarization observables [53]; in those analyses, a specific set of the LDMEs was applied -they had been obtained in Ref. [16].
Furthermore, the event activity dependence of J/ψ production yield was explored in the CGC framework in Ref. [54]. In the CGC picture, one can attribute high multiplicity events to rare parton configurations in hadrons and nuclei, participating in scattering processes. Such rare parton configurations can be set up by changing the Q 2 s at large x [55][56][57], giving stronger multiple-scattering effects to the heavy quark pair. As to other approaches, Refs. [58,59] discussed the saturation effect on quarkonium production, while Ref. [60] studied quarkonium production in a dilute system where the saturation effect is less important. Percolation approach [61] shares the concept with the saturation approach. In addition, PYTHIA8 event generator has been used to study multi-parton interactions [62], which are partially included in the CGC, and another event generator EPOS3 can also examine the saturation effect [63]. Of particular importance is that recent LHC data [10] have shown that the initial state effects without final state flow effect can describe the event activity dependence of J/ψ yield quantitatively.
The polarization of the J/ψ meson can be analyzed using its decay into lepton's pair. Three polarization parameters λ θ , λ φ , and λ θφ can be extracted from data by measuring an angular distribution of the leptons. Measurements of the polarization parameters at the Tevatron [64,65] and later at RHIC [66][67][68] and the LHC [69][70][71][72] showed that J/ψ is produced with almost no polarization from low to high p ⊥ in hadronic collisions. On the other hand, the first theoretical predictions based on collinear factorized color singlet production channel at leading order (LO) suggested substantially nonzero polarization at high p ⊥ [73]. The inclusion of the octet states in the NRQCD expansion with SDCs calculated using NLO pQCD [15][16][17][18] reduced this discrepancy (known as "J/ψ polarization puzzle"), but no full agreement with polarization measurements has been reached yet. The polarization of hadronic inclusive J/ψ production was also studied within the CGC+NRQCD framework [53]; good overall description of minimum bias data on J/ψ polarization in pp collisions at forward rapidities from LHCb and ALICE [69,71,72] was obtained. Moreover, the STAR experiment has recently reported [68] nice agreement of the CGC+NRQCD predictions with data in minimum bias pp collisions.
In this paper, we study the polarization of J/ψ in high multiplicity small collision systems, which should provide us with useful information on the fundamental quarkonium production mechanism. For this purpose, we use the CGC+NRQCD framework following an approach from Ref. [54] and clarifying how much the initial state saturation effect could impact the polarization of J/ψ as charged multiplicity increases. Up to now, there are no experimental data for J/ψ's polarization as a function of charged particle multiplicity, so this paper aims to provide predictions for the future LHC experiments including High-Luminosity LHC [74] with the kinematical conditions used in similar measurements for unpolarized J/ψ production.
The paper is organized as follows. In section II, we briefly present calculations of the polarization parameters in the CGC+NRQCD framework. In section III, we discuss the charged hadron production within the CGC framework. Section IV demonstrates numerical results about the polarization of J/ψ in small collision systems. We then conclude in section V.

II. QUARKONIUM PRODUCTION IN THE CGC+NRQCD FRAMEWORK
We begin with reviewing steps toward calculating quarkonium production (particularly J/ψ) in a dilute-dense system in the CGC+NRQCD approach.
The spin density matrix elements, depending on the quarkonium transverse momentum p ⊥ and rapidity y, are calculated in the NRQCD factorization framework [13]: where the matrix indices i and j represent helicity of the J/ψ in the amplitude and conjugated amplitude, respectively. In the above expression dσ κ /d 2 p ⊥ dy denotes SDCs for heavy quark pair production in a specific spin-color state represented by κ. O κ are LDMEs that describe the nonperturbative transition between a cc pair and J/ψ meson. As nonperturbative quantities they need to be determined by data fitting. Meanwhile, LDMEs should be processindependent in the NRQCD factorization. For the J/ψ production, most dominant intermediate states in the quark velocity expansion are where standard spectroscopic notation is used for state κ: 2S+1 L [C] J with [C] denoting a color state (singlet [1] or octet [8]).
In the CGC framework, the SDC for a color singlet cc production in a dilute-dense system (e.g. pA) is given by [53]: and for a color octet state by Matrices Γ κ ij and G κ ij describe couplings of the Wilson lines to the cc pair. For the details of calculations and explicit expressions, see Appendix B of Ref. [53]. Their traces were calculated earlier in Refs. [49,51]. In expressions above πR 2 p(A) denotes the effective transverse area of the projectile proton (target nucleus) [51]. In this paper, as we shall focus on ratios of spin density matrix elements, the overall normalization factors are not considered below. Forward scattering amplitudes N x andÑ x correspond to the Fourier transform of the dipole correlator of light-like Wilson lines in the fundamental and adjoint representation respectively at For the LDMEs we will apply values obtained in Ref. [16] by fitting NLO collinear factorized pQCD+NRQCD results to the Tevatron high p ⊥ prompt J/ψ data. This set of LDMEs was also used in the previous CGC+NRQCD studies, where a good description of data was obtained for the p ⊥ spectra of J/ψ in pp and pA collisions [51,52] and J/ψ polarization in pp collisions [53]. The color singlet LDME is estimated using the value of the quarkonium wavefunction at the origin in a potential model [75]: where m c is the charm quark mass for which we take the value 1.5 GeV. One practical reason for choosing this set of the LDMEs is that the numerical values for each of the color octet states are positive, so our computations are not conflicted with cancellations among the color octet contributions.
From the spin density matrix elements (1) one obtains the polarization parameters: Those three parameters can be measured experimentally by considering the leptonic decay of J/ψ in its rest frame.
By Ω = (θ, φ) we denote the solid angle of the positive lepton w.r.t. arbitrary chosen vectors X, Y, Z. If Y is chosen to be perpendicular to the hadronic plane, the angular distribution of the measured positive lepton is given by [76,77]: The angular distribution (6) depends on the orientation of vectors X, Z w.r.t. hadrons' momenta. There are several choices of frames used in literature. We will use the most popular two: the Collins-Soper [78] frame and the recoil (also called helicity) frame [79]. For the explicit definition of these frames see for example Refs. [53,80].

III. INCLUSIVE HADRON PRODUCTION IN THE CGC FRAMEWORK
This section summarizes the calculation of charged hadron multiplicity in the CGC framework; for more details, see Ref. [54].
The pseudo-rapidity dependence of charged hadron multiplicity can be obtained in our approach as: where we have skipped the overall normalization constant as we are interested only in ratios of this quantity. dσ g /d 2 p g⊥ dy g denotes the gluon production cross section calculated in the CGC framework. For the fragmentation function D h (z) we take the NLO parametrization of the Kniehl-Kramer-Potter fragmentation function at µ = 2 GeV [81]. The hadrons carry z fraction of the gluon's transverse momentum: p h⊥ = zp g⊥ . The Jacobian J(y h → η) accounts for the transformation between rapidity y h and pseudo-rapidity η of the hadron. In addition we assume that the rapidity of the hadron is the same as the rapidity of the parent gluon: y h = y g . For the typical mass of a light hadron we take m h = 0.3 GeV. The integration over z in (7) is constrained from below by the kinematical conditions x 1,2 ≤ 1, where x 1,2 = p g⊥ e ±yg / √ S. In the CGC framework, the cross section for the inclusive gluon production reads [82,83]: For the numerical purposes we restrict values of p h⊥ in (7) to 0.1 GeV ≤ |p h⊥ | ≤ 10 GeV. This cut has negligible effect on the results since the cross section decreases rapidly with p g⊥ and, in addition, we shall consider only ratios of dN ch /dη. The dipole correlatorsÑ x1 (k ⊥ ) and N x2 (k ⊥ ), which appear in (3), (4), and (8), can be obtained by solving the BK equation. In this paper we use numerical solutions of the running coupling BK (rcBK) equation [84] in momentum space [85]. Initial conditions at x 0 = 0.01 were chosen according to McLerran-Venugopalan (MV) model [86,87] which in the position space reads This initial condition is parametrized by saturation scale Q 2 s0 in the proton/nuclei at x 0 = 0.01. Q 2 s0 is proportional to Q 2 0 -the saturation scale of the proton for minimum bias events. Values Λ = 0.241 GeV, γ = 1.119 and Q 2 0 = 0.168 GeV 2 were obtained in Refs. [88][89][90] by fitting to HERA DIS data with x ≤ x 0 .
In the dense saturation regime at high energies, hadron multiplicity per unit of rapidity and transverse area in impact parameter space in the central rapidity region for symmetric pp collisions can scale as follows [91,92]: dN ch /dη ∝ k 2 ⊥ S ⊥ /α s ∼ Q 2 s S ⊥ /α s with S ⊥ being an effective transverse area for impact parameter, k ⊥ the average transverse momentum of produced hadrons. The above scaling is a consequence of Eq. (8). High multiplicities can be achieved due to the hot spots where gluons are highly occupied, giving larger saturation scales; such rare parton configurations are less likely to happen compared to minimum bias events. It has been demonstrated numerically in [54] that increasing the initial saturation scales gives rise to high multiplicities. Therefore, we use the value of Q 2 s0 in (9) as a parameter controlling charged multiplicity of a given event. We assume that for pp collisions the initial saturation scales in both protons are the same: Q 2 s0,proton 1 = Q 2 s0,proton 2 ≡ Q 2 s0,proton 2 . It is equal to Q 2 0 for minimum bias events and ξQ 2 0 , ξ > 1 for high multiplicity events. So we have: for minimum bias events, where l.h.s. is calculated using Eq. (7) and . For pA collisions, the initial saturation scale in the nucleus, which is embedded in the MV initial condition, is expected to be larger than in the proton, Q 2 s0,nucleus = νQ 2 s0,proton with ν > 1. In this case the scaling of dN ch /dη is expected to be more complicated than symmetric pp collisions because there can be the dependence 2 Note that Q 2 s0 is the initial saturation scale, hence at x 0 = 0.01. The "evolved" saturation scales probed in the given event are determined by the rcBK equation and depend on the x values. By increasing the initial saturation scale, we also increase the saturation scales at smaller x, which lead to higher event activities.
on ln(Q s,proton /Q s,nucleus ) [93]. Nevertheless, we drop such a explicit logarithmic correction for simplicity and instead tune the initial saturation scales by hand [52,54,94,95]. We then have for minimum bias pA collisions: For higher multiplicity events, we assume that the initial saturation scales grow by the same factor ξ in the proton and nucleus: with ξ > 1. The parameter ν, being the ratio of initial saturation scales in the nucleus and proton, is not derived from first principles, i.e., QCD. Some previous studies [52,[94][95][96][97] suggest that it should be between 2 and 3 for a heavier target, i.e., Pb. We will calculate results both for ν = 2 and ν = 3, and treat this factor as an additional source of uncertainty. As we described above, we use the rcBK solution when x < x 0 = 0.01. For x > 0.01, we employ an extrapolation of the solutions of the rcBK equation by requiring that the corresponding integrated gluon distribution matches that in the collinear factorization framework, see [51,54] for more details.
In what follows, we will calculate the ratios between the yield in high-multiplicity pp (pA) collisions and that in minimum bias events: where X = p, A and multiplicities are defined in Eqs. (10)- (13). Both multiplicities are integrated over the same pseudorapidity range. The polarization parameter λ θ (5) in pA collisions can be written as: where all density matrix elements (1) are evaluated using the same initial saturation scales. Just like for the light hadrons multiplicity, the parameter ξ 'controls' the activity of the events. Polarization parameters λ φ and λ θφ (5), are calculated similarly to (15). The calculation for pp collisions is analogous, with ν = 1, as in Eq. (11). As we have shown above, the J/ψ production cross-section in the CGC+NRQCD framework (1) does not include any medium effect but the saturation effect. In this paper, we assume that the bound state formation happens at a later stage without modifying the LDMEs, which could be acceptable only when considering J/ψ production [97]. Thus any modifications in the J/ψ's polarization in pA collisions and high multiplicity events can be attributed to the saturation effect at a short distance. In the following section, we shall clarify the impact of the saturation effect on the polarization of J/ψ.

IV. RESULTS
Since experimental data for polarization parameters in high multiplicity events are not available at present, we will give our predictions with kinematical conditions probed in the related measurements. We choose energies from Run 2 of the LHC, √ S = 13 TeV for pp collisions and √ S = 8.16 TeV for pA collisions. We constrain the charged hadrons' pseudorapidity to the window |η ch | < 1, which was used in several ALICE measurements of J/ψ's yield as a function of charged-particle multiplicity [5,9,10]. The J/ψ rapidity we restrict to the window 2.5 < y J/ψ < 4, where positive rapidity is defined by the proton-going direction (for pA collisions). This setup provides the dilute projectile -dense target scattering, where the CGC formalism is particularly effective. Moreover, this range of rapidity has already been used for the measurements of J/ψ polarization [69,72]   In order to estimate uncertainties of our predictions, we use errors of the LDMEs (see section II) -we vary each LDME by the value of its error and calculate maximal and minimal value of the polarization parameters obtaining a band of uncertainty. For the pA collisions also the uncertainty from unknown ratio ν = Q 2 s0,nucleus /Q 2 s0,proton was included (see (13) and its discussion).
Our results for three polarization parameters are displayed in Figure 1. We provide predictions for two leptonic frames: Collins-Soper frame (left column) and helicity frame (right column). For each polarization parameter, the predictions for pp (red hatched band) and pA (green filled band) collisions are shown. We have checked that results do not depend significantly on energy √ S. This is consistent with the previous result that the correlation between the yield of J/ψ and charged-hadron may be characterized by Q 2 s , so that the √ S dependence is weak [54]. The difference between pp and pA results in Figure 1 comes entirely from the different type of target. Values of all parameters do not exceed (in the module) 0.1, which means almost unpolarized production of J/ψ. The parameter λ θ in both frames decreases with the event activity -this behavior is particularly well seen in the helicity frame where λ θ is equal to 0.1 at minimum bias events and decrease with the growing event activity to negative values. In the Collins-Soper frame, the λ θ in pA collisions is slightly smaller than in pp collisions. In the helicity frame, the opposite is true. The λ φ parameter is negative and very close to 0 in both frames. The λ θφ parameter is larger than λ φ , has a different sign in both frames, and very slightly decreases with the event activity. We note that the small polarization of produced J/ψ in the CGC+NRQCD approach was also predicted for the inclusive events (see [53]).
Using our approach, we predict a slight difference between the J/ψ polarization in pp collisions and that in pA collisions. All the differences are within the uncertainty bands. This characteristic behavior can be traced back to the feature of our model. Firstly, the saturation effect participates in partonic scatterings at short distance, where a produced cc pair does not experience the hadronization effects. Secondly, the same values of the LDMEs were used both in pp and pA collisions -this universality is suggested by the good description of the p ⊥ spectra of J/ψ in pp and pA collisions with the same set of the LDMEs [51,52]. Finally, in our model, there is no final state interaction of J/ψ. One should note that the similarity between pp and pA collisions is not present in data for J/ψ's yield and mean p ⊥ at high multiplicity events. If we use the same set up shown in this paper, the CGC predictions for those observables differentiate between pp and pA collisions [54] because the J/ψ's yield and p ⊥ are driven by combinations of the saturation scales, which are very different in the proton and nucleus. On the other hand, the polarization observables are not sensitive to those differences, as we have shown here. Therefore, we conclude that J/ψ's polarization data will provide a complementary test of the CGC theory at the high multiplicity events.

V. SUMMARY
In this paper, we have demonstrated the polarization parameters of J/ψ as a function of charged-particle multiplicity in pp and pA collisions in the CGC framework. Our work is a continuation of the previous analyses, where the CGC framework was applied to high multiplicity events [54] and (joined with the NRQCD) to J/ψ polarization studies [53]. In this approach, the event activity is controlled by the value of the saturation scale encoded in the solution of the rcBK evolution equation.
We predict a very small polarization of the produced J/ψ in minimum bias pA collisions, which is almost consistent with the J/ψ polarization in pp collisions. According to our model, the polarization parameters decrease slowly as the event activity increases. We found that the J/ψ's polarization has a weak dependence on the system size and √ S, albeit the saturation scale in the target is different between pp and pA collisions. This is contrary to the other observables including J/ψ yield and mean p ⊥ [54].
The polarization measurements of quarkonium in high multiplicity events provide an independent insight into the QCD dynamics behind high multiplicity events. Given the data which are currently available, such observable should be accessible experimentally at the LHC. Together with other measurements, like quarkonium's yield and mean p ⊥ , polarization will help us understand the role of initial state interactions, hadronization, and final state interactions in quarkonium production processes.