$J/\psi$ polarization in the CGC+NRQCD approach

We compute the $J/\psi$ polarization observables $\lambda_\theta$, $\lambda_\phi$, $\lambda_{\theta\phi}$ in a Color Glass Condensate (CGC) + nonrelativistic QCQ (NRQCD) formalism that includes contributions from both color singlet and color octet intermediate states. Our results are compared to low $p_T$ data on $J/\psi$ polarization from the LHCb and ALICE experiments on proton-proton collisions at center-of-mass energies of $\sqrt{s}=7$ TeV and 8 TeV. Our CGC+NRQCD computation provides a better description of data for $p_T \leq 15$ GeV relative to extant next-to-leading (NLO) calculations within the collinear factorization framework. These results suggest that higher order computations in the CGC+NRQCD framework have the potential to greatly improve the accuracy of extracted values of the NRQCD universal long distance matrix elements.


Introduction
The study of heavy quarkonium states in QCD is an essential ingredient in developing our understanding of the subtle interplay of short and long distance physics in QCD. However even though the simplest J/ψ meson was discovered more than 40 years ago, key features of how this state is produced in high energy collisions continue to elude us. A prominent example is the polarization of the J/ψ, which appeared to differ significantly from theoretical expectations.
The theoretical models describing the production of heavy flavors in QCD rely on the factorization between the hard process governing the production of the heavy quarkantiquark pair and the soft processes governing the hadronization of this pair into quarkonium states such as the J/ψ. The former can be computed in perturbative QCD while the latter is intrinsically nonperturbative and can be determined only from models or effective field theory approaches. For instance, in the color evaporation model [1], the perturbative production of the heavy quark-antiquark (QQ) pair with mass M is followed by its nonperturbative hadronization to the final state meson with a universal transition probability for all M below the mass threshold of producing two open flavor heavy mesons. In the color singlet model, the QQ pair is produced in a color singlet state before hadronizing into the quarkonium state. The QQ wave function in this approach is computed at zero separation between the quark-antiquark pair. The most sophisticated approach to describe the hadronization of heavy quarkonia is nonrelativistic QCD [2], an effective field theory valid in the limit of very heavy quark masses. NRQCD employs systematic power counting in the relative velocity of the QQ to determine the long distance matrix elements (LDMEs) of the dominant nonpeturbative operators contributing to the formation of quarkonia. For the J/ψ, the LDME of the color singlet channel is dominant in the power counting of the matrix elements, with significant contributions also arising from several color octet channels [2][3][4].
There has been a large amount of work in recent years developing the NRQCD formalism and applying it to quarkonium measurements in a wide range of experiments [5]. Our focus will be on quarkonium production in proton-proton collisions. In this case, the matrix elements for the production of QQ pairs in both color singlet and color octet states were calculated within the collinear factorization formalism including next-to-leading (NLO) perturbative corrections. These results were employed by several groups to extract the nonperturbative LDMEs from comparisons of the cross-sections to data [6][7][8][9][10].
For the production of J/ψ's with large transverse momentum p T , the most important contribution in NRQCD at leading order (LO) in α s comes from the 3 S [8] 1 channel, which suggests that the produced J/ψ's should be transversely polarized [11]. The polarization of the J/ψ is extracted from the angular distribution of positively charged (by convention) leptons in the decay of J/ψ into muons (J/ψ → µ + µ − ), that is parametrized by the coefficients λ θ , λ φ and λ θφ . Transversally (longitudinally) polarized J/ψ's have λ θ = 1(−1), λ φ , λ θφ = 0 while all coefficients are zero for unpolarized J/ψ's. Measurements of the J/ψ polarization by the CDF Collaboration at the Tevatron [12,13], as well as the ALICE [14,15], LHCb [16] and CMS [17] experiments at the LHC, showed that the J/ψ has weak or no polarization. This stark disagreement between the leading NRQCD expectation and collider data has been dubbed the "J/ψ polarization puzzle".
Extensions of the LO NRQCD J/ψ polarization studies in hadronic collisions to next-toleading order (NLO) in collinearly factorized perturbative QCD (pQCD) approaches have been discussed in [18,19]; further studies including feeddown contributions from higher states were also discussed in [20,21]. The latter computations are important because there are no available experimental data for direct J/ψ production; the polarization is measured either for "prompt" (including feeddown from the higher excited charmonium states ψ(2s), χ cJ . . .) [12,13,16,17] or "inclusive" production (prompt plus additional contributions from bottom meson decays) [14,15]. The conclusion of these NLO computations was that the 3 S J channels have a large cancellation between their transverse and longitudinal polarization components [19], thereby providing a possible explanation for the lack of J/ψ polarization at high p T .
However the J/ψ produced in proton-proton collisions are also weakly polarized at low values of p T where the collinear factorization formalism may not be applicable. At low p T 's at collider energies, large α s ln(1/x) contributions arise at higher orders that may not be fully accounted for in collinear factorization frameworks. Another source of O (1) contributions are higher twist multiparton matrix elements that are large at low p T . Small x kinematics is also accessed in either the projectile or target at forward rapidities. For instance, the LHCb [16] and ALICE [14,15] experiments measure J/ψ's at the forward rapidities of 2 < y < 4.5 (LHCb) and 2.5 < y < 4 (ALICE) for p T < 15 GeV, providing access to x values down to x ∼ 10 −4 in one of the protons.
The contribution of large α s ln(1/x) contributions as well as leading higher twist contributions to quarkonium production can be computed systematically in the Color Glass Condensate (CGC) effective field theory (EFT) [22][23][24][25][26]. This EFT treats large x degrees of freedom in the two hadrons as static color sources that are coupled to dynamical gauge fields at small x. Physical quantities such as heavy quark pair cross-sections are computed in a two-step procedure; they are first computed for a fixed distribution of color sources, in the gauge field background, and subsequently averaged over a gauge invariant distribution of color sources. The separation scale in x between large x sources and small x fields is arbitrary. However the requirement that physical quantities do not depend on this separation scale leads to a renormalization group equation, the JIMWLK equation [27][28][29][30], describing the change in the distribution of color sources with decreasing x. A key feature of this approach is a dynamically generated saturation scale Q s (x) [31][32][33][34] that grows both with decreasing x and with increasing nuclear size. When Q 2 s Λ 2 QCD , one can employ weak coupling methods to compute cross-sections even at low p T since α S (Q 2 s ) 1. The heavy quark pair production cross-section for proton-proton and proton-nucleus collisions was computed in a "dilute-dense" approximation of the CGC in [35,36]. This approximation corresponds to systematically keeping lowest order terms in an expansion of the smaller of the color charge densities of the two colliding hadrons, and terms to all orders in the larger of the two color charge densities -hence the moniker dilute-dense. It is strictly valid for forward proton-proton collisions or in proton-nucleus collisions in p T and rapidity windows that are consistent with this expansion 1 . The dilute-dense results for heavy quark pair production were later used to compute the short distance cross-section (SDC) for Onium production in a CGC+NRQCD approach [37]. Numerical results for J/ψ production in p + p at RHIC and LHC collisions were presented in [38] and likewise for p + A in [39]. More recently, results for J/ψ production in high multiplicity p + p and p + A collisions were obtained in [40]. The CGC+NRQCD approach describes quite well the systematics of the p T and rapidity dependence of the J/ψ yields in both proton-proton and proton-nucleus collisions within theory uncertainties. We note that the CGC has been applied to compute the SDC in quarkonium production in hadron-hadron collisions 2 approximating the LDME with the color evaporation model [45,46] and variants thereof [47,48].
In this paper, we shall extend the CGC+NRQCD analysis of [37] and [38] to address the J/ψ polarization puzzle. We will begin by first relating the coefficients of the angular distribution of the positively charged leptons produced in J/ψ leptonic decays to helicity dependent quarkonium cross-sections. These coefficients are frame dependent and are typically presented as such by the experiments though frame independent combinations of these can also be extracted. In Section 3, we will write down the explicit expressions for the helicity dependent SDCs in the CGC+NRQCD framework. Numerical results for the coefficients of the angular distribution are presented in Section 4. We observe that the agreement of the theoretical computations with the data is quite good for the low p T LHC data. In particular, we observe that there is also a cancellation between the polarizations of the 3 S [8] 1 channel with that of the 3 P [8] J channel in the CGC+NRQCD framework, even for the LO in α s impact factor. We will end with a summary and outlook on further work. Several details of the computation, and additional results, are provided in two appendices.
2 Angular distribution of J/ψ leptonic decay and helicity dependent cross-sections In order to extract the polarization of the J/ψ, first consider the leptonic decay of J/ψ in its rest frame [49]. The angular (θ, φ) distribution of the positive lepton is obtained by fixing a frame X, Y, Z with respect to which the angles θ and φ are measured: Here p l + is the four-momentum of the positively charged lepton created in the J/ψ's decay and | p l + | is a length of its three-momentum in the J/ψ's rest frame. The unit four-vectors X, Y, Z span the subspace perpendicular to J/ψ's momentum p µ and are perpendicular to each other: The orientation of the vectors X, Y, Z with respect to the momenta P 1 , P 2 of the incoming hadrons depends on a choice of frame. In this paper, we will consider two frames which are often used in the literature: the Collins-Soper [50] frame and the recoil (or helicity) frame [51]. In both frames, the Y four-vector is chosen to be perpendicular to the hadron plane, In the J/ψ center-of-mass frame, the corresponding three-vector is chosen to be 3 , In order to define the X and Z four-vectors, we first introduce projections of the hadron momenta on to the J/ψ's four-momentum p µ : where p µ p µ = M 2 . Then X and Z are linear combinations ofÃ andB, where the recoil and Collins-Soper frames are defined by the values of the coefficients α x,z , β x,z . Explicit expressions for these coefficients were given in [52]. For completeness, they are listed in Appendix A.
Since the X and Z vectors lie in the plane of the incoming hadrons, the positively charged lepton's angular distribution in the J/ψ rest frame can be parameterized by three coefficients as [53] dσ J/ψ(→l + l − ) dΩ ∝ 1 + λ θ cos 2 θ + λ φ sin 2 θ cos 2φ + λ θφ sin 2θ cos φ, (2.6) where Ω = (θ, φ) denotes the solid angle of the positive lepton in Eq. (2.1). The coefficients of this angular distribution computed using the spin density matrix elements for J/ψ production can be expressed as [53] λ θ = dσ 11 − dσ 00 dσ 11 + dσ 00 , The cross-section dσ ij corresponds to the product of the amplitude for inclusive production of a J/ψ with helicity i in the amplitude and helicity j in the complex conjugate amplitude. Hence dσ 11 (dσ 00 ) can be interpreted as the cross-section for the production of J/ψ with helicity h = +1 (h = 0). The unpolarized cross-section is given by a sum of contribution of three helicity states: h = +1, h = 0 and h = −1: which is the trace of the spin density matrix. The values of the coefficients λ θ , λ φ and λ θφ depend on the choice of frame-the choice of the X and Z axes. One can construct out of these coefficients frame-independent quantities as well, as discussed in [54][55][56][57][58]. We will present results for two of these invariants, which are defined as (2.9)

Computation of the spin density matrix in CGC+NRQCD
In what follows, we will write down the spin density matrix elements dσ ij in the CGC+NRQCD formalism. The spin density matrix elements can be expressed as [2] where O κ are the NRQCD long distance matrix elements (LDMEs). The SDCs dσ κ ij describe the production of the cc pair in a given quantum state κ = 2S+1 L denotes either the singlet [1] or the octet [8] color state. The LDMEs describe the nonperturbative transition of the cc pair into the J/ψ state; these are process independent and can be determined by fitting experimental data. For J/ψ production, the leading contribution to the sum in Eq. (3.1) comes from the states 3 S Based on NRQCD velocity scaling rules [2], the spin of the J/ψ is the same as that of the intermediate cc pair if it is produced via 3 S J states. See also [52,59] for further discussion.
In the CGC effective field theory, the SDC's are given by the expressions [37,38], for the color octet channels and for the color singlet channels. In these expressions, N Y denotes the forward scattering amplitude corresponding to the Fourier transform of the "dipole" correlator of lightlike Wilson lines in the fundamental representation [24], πR 2 p is the effective transverse area of the proton [38] and ϕ p is an unintegrated gluon distribution inside the proton: whereÑ Y is the Fourier transform of a dipole correlator, but in this case with lightlike Wilson lines that live in the adjoint representation. These dipole forward scattering amplitudes, N Y andÑ Y , are obtained by solving the running coupling Balitsky-Kovchegov (rcBK) equation [60,61] in momentum space, as a function of x, with McLerran-Venugopalan (MV) initial conditions [33,34] specified at an initial large scale x 0 = 0.01 [62]. For x > 0.01, we employ an extrapolation of the solutions of the rcBK equation [38] which is constrained by requiring that the corresponding integrated gluon distribution matches that in the collinear factorization framework. We refer the reader to [37,38], and the references therein, for details of the derivation of these expressions. The novel feature here is the unwrapping (so to speak) of the helicity integrated expressions derived in [37] to extract the helicity dependent functions Γ κ ij and G κ ij . The procedure is outlined in Appendix B, where we provide the detailed expressions for these functions as well.

Numerical results
We will now explicitly compute the expressions in Eqs. (3.3) and (3.4) and use these to determine the angular distribution coefficients specifying J/ψ polarization. In our parameter set for the numerical computations, we will set the charm mass to be m c = 1.5 GeV-nearly one half of the J/ψ mass. The value of the color singlet LDME is estimated using the value of the wavefunction at the origin in a potential model [63]: For the color octet LDMEs, we employ the values obtained in Ref. [19] by fitting NLO collinear factorized pQCD + NRQCD results to the Tevatron high p T prompt J/ψ yields data: We will not use other sets of LDMEs extracted at NLO [20,64,65] as they contain negative values for some of the LDMEs 4 .
The solution of the rcBK equation employs the code of Albacete et al. [62] with MV initial conditions and the initial input parameters γ = 1, Q 2 s0,proton = 0.2 GeV 2 , α f r = 0.5 and C = 1; these were determined from fits to the HERA DIS data [62]. We have checked that our results for the angular coefficients, being ratios of cross-sections, are insensitive to the values of these parameters. The theoretical errors we quote therefore are for the angular coefficients (collectively denoted henceforth as λ) and are obtained by varying the LDME values by their statistical uncertainties and by taking the minimal/maximal value of the obtained set.

Spin density matrix elements in specific color channels
We begin with a comparison of the contributions from different channels to the spin density matrix elements. In Figure 1, we plot the the matrix elements with ij = {00, 11} for different quantum states. The rapidity interval is chosen to be 2.5 ≤ Y ≤ 4 and the centerof mass energy is √ S = 7 TeV. One sees immediately that the 1 S [8] 0 state is the dominant channel for both matrix elements, as was also seen in the NLO collinear pQCD calculations [18]. The color singlet state 3 S [1] 1 contribution is similar to that from the 3 P [8] J state in the Collins-Soper frame; they are both larger than the 3 S [8] 1 state at low p T and then decrease rapidly with p T such that 3 S [8] 1 starts to become important at higher p T . This is also the case for σ 11 in the recoil frame but σ 00 for the 3 S 1 state decreases as fast as the other channels. A similar behavior can be seen in Figure 2 of [18]. This explains why in the in the recoil frame at high p T we have strong transverse polarization in the 3 S [8] 1 channel. Even though the 1 S [8] 0 state is numerically dominant by far in the unpolarized crosssection in Eq. (2.8), it gives vanishing values for the polarization coefficients λ θ , λ φ and λ θφ because the produced quark-antiquark pair has no spin and orbital momentum -for further discussion, see [66]. We can write the helicity SDCs in this state as in other cases ,

Results for the λ polarization coefficients in CGC+NRQCD
In Figure 2, we show results for all three angular distribution coefficients λ compared to data in the Collins-Soper frame (left column) and in the recoil frame (right column). The following data were used 5 : LHCb at 7 TeV [16], ALICE at 7 TeV [14] and ALICE at 8 TeV [15]. Both LHCb and ALICE measured angular coefficients in the rapidity window 2.5 < y < 4. Note that ALICE data are obtained for inclusive J/ψ production, so they also contain contributions from B-meson decays. However the contribution from B-meson  Figure 2. The angular distribution coefficients λ θ (first row), λ φ (second row) and λ θφ (third row) in the Collins-Soper frame (left column) and in the recoil frame (right column) as functions of the J/ψ transverse momentum p T . Data are from the LHCb experiment at 7 TeV [16], the ALICE experiment at 7 TeV [14] and ALICE at 8 TeV [15], all in the 2.5 < y < 4 rapidity window. Note that the ALICE data are for inclusive J/ψ production, containing contributions from B-meson decays.
decays are on the order of a few percent at low p T [67], so we can neglect them here 6 .
The polarization parameter λ θ is measured to be close to zero, indicating that the J/ψ are mostly unpolarized. At small p T , our results prefer a small transverse polarization (λ θ > 0). The data on the other hand seem to prefer a weak longitudinal polarization (λ θ < 0) albeit it should be noted that there is considerable variation between the experiments, with LHCb and ALICE 7 TeV data showing negative central values and ALICE 8 TeV data showing positive values. The data are consistent with each other to 1σ accuracy. At higher transverse momentum (p T 6 GeV), our results agree with data within two standard deviations. We note that the agreement of our theory results with the data in this kinematic region is significantly better than two of the three the collinearly factorized NLO pQCD+NRQCD computations; it is however close to the results of [19]. The compilations shown in [16] and [15] comparing NLO pQCD+NRQCD and color singlet model results to data demonstrate that there is considerable variation between the data for λ θ and those NLO pQCD+NRQCD theory results.
For the λ φ coefficient, we obtain very good agreement with data; the data are within the CGC+NRQCD theory band for both frames. For λ θφ , we obtain very good agreement with the LHCb data in the Collins-Soper frame. Our results are higher than the ALICE 8 TeV data though within 1σ accuracy; we note that there is tension between the LHCb and ALICE data at this level of accuracy. In the recoil frame, the ALICE data are well described except for high p T , where data are systematically above the CGC+NRQCD predictions. In contrast, the agreement with the LHCb data is good at the higher p T while we are slightly below at low p T .
The take away message here is that the experimental values, as well as our theory results, for λ φ and λ θφ , are consistent with zero. Our description of data for λ φ and λ θφ is significantly better than the NLO pQCD+NRQCD calculations of [18]: in fact, both the color singlet model and NLO pQCD+NRQCD approaches predict polarization coefficients with absolute values significantly larger than measured at LHCb and ALICE. A similar conclusion can be drawn for results obtained by a second group performing NLO pQCD+NRQCD analyses [20]. While, as noted, the third group [19,21] obtain a good description of data for λ θ in the recoil frame for p T > 5 GeV, no results are provided for the λ φ and λ θφ coefficients in this frame, or for any of the angular coefficients in the Collins-Soper frame.
In Figure 3, we show our results for the frame invariant quantities λ inv defined in Eq. (2.9), which are compared to the ALICE data [15] at 8 TeV. Once again, one observes that CGC+NRQCD provides a good description of λ (1) inv though indeed the error bars in the data are considerable. We also provide predictions for the λ (2) inv coefficient that was proposed in [57]. Since this coefficient is sensitive to all three of the frame-dependent polarization coefficients λ θ , λ φ , λ θφ , we believe it may provide an additional constraint both for theoretical studies and experimental measurements.
We can however compute the angular polarization variables as a function of rapidity 6 The good agreement between the LHCb data for prompt production and the ALICE data for inclusive production affirms this statement.   and p T and compare these to data to check the quality of agreement with varying rapidity. In Figure 4, we show the rapidity dependence of all three coefficients plotted in both Collins-Soper and recoil frames. They are plotted for three transverse momentum values, p T = 2.5, 4.5, 8.5 GeV represented by three bands on each plot (two of these are shifted by 0.5 and 1 unit for better visibility.) One observes that the CGC+NRQCD results are almost completely independent of rapidity in the range plotted. While this also appears to be the case for λ φ and λ θφ , the experimental values of λ θ seem to have some rapidity dependence (within uncertainties) for p T = 2.5 and 4.5 GeV. In general, one may conclude that at higher rapidities our dilute-dense CGC+NRQCD predictions are closer to data as they should be. In Figure 5, we show LHCb data [16] collected for the highest rapidity window, 4 < y < 4.5 and compare them with our predictions. One sees slightly better agreement for λ θ than that seen in the wider rapidity window of 2.5 < y < 4 plotted in Figure 2. Though tempting, it would be premature however to conclude that this better agreement is primarily due to the dilute-dense approximation being better satisfied.
Note that our dilute-dense approximation in the CGC EFT assumes asymmetric treatment of the two colliding protons: the target is treated as a dense parton system for the resolved transverse momentum in the target. Similarly, the projectile is assumed to be dilute. This assumption is best satisfied for forward J/ψ production at low p T . In Fig. 6, we show the kinematic range of x values probed in the projectile and target. For the projectile, the x values are quite large; however, as discussed in [38], the unintegrated gluon distributions can be constrained by a smooth matching to the collinear pQCD gluon distribution. In the case of the target, the very small values of x ∼ 10 −5 − 10 −4 motivate its treatment as a dense system. The kinematics of projectile and target suggest therefore that the application of the CGC dilute-dense formalism is appropriate. For p T ≥ 8 GeV, one is starting to probe x values where replacing the unintegrated k T distribution with the gluon parton 5 . distribution will begin to receive large corrections. This situation will be exacerbated for 4 < y < 4.5, where the hybrid formalism may be more appropriate; computations in such  Figure 5. The angular coefficients λ θ (first row), λ φ (second row) and λ θφ (third row) in the Collins-Soper frame (left column) and recoil frame (right column) as functions of J/ψ transverse momentum p T . Data are from the LHCb experiment at 7 TeV [16], for the highest rapidity window 4 < y < 4.5.
a framework matched to NRQCD are not available at present.

Discussion
The computation of the helicity SDCs in our paper are performed in the CGC weak coupling framework which, in principle, differs significantly from those in the NLO collinear factorization framework which has a different kinematic window of applicability. However there is also a significant regime of overlap between the two approaches. As it was shown in [68], the dilute-dilute limit of the dilute-dense limit we have considered here is equivalent to the k T factorization formalism. Further, it was shown in [68] that for k T → 0 one can express the heavy-quark pair cross-section as a convolution of the LO collinear pQCD matrix element for gg → qq scattering and the product of small x gluon parton distributions. Thus in including the NLO BK evolution equation with running coupling [69][70][71] in our computation of the cross-sections, we are including important pieces of the leading NLO, NNLO,· · · collinear pQCD contributions at small x, as previously also emphasized in [72,73]. The matching between the two frameworks will fail when p T becomes sufficiently large that p T dependent contributions that are subleading in x begin to play a role. It is at present not known analytically where in p T this mismatch occurs. This would require higher order computations in both frameworks than currently available. However one can see how good the matching is phenomenologically and where they begin to differ. Such a comparison was performed for J/Ψ production in [38]; good agreement was obtained for p T 10GeV. For the case of J/Ψ polarization, to illustrate the fact that our CGC computation includes important NLO collinear contributions, we will compare our result with the NLO pQCD+NRQCD calculation by Chao et al. [19,21]. We define the coefficient which is a polarization parameter λ θ calculated for the given channel κ 7 . Note that it 7 We introduce absolute value in the denominator of (4.2) because NLO pQCD dσ 3 P does not depend on the values of the LDMEs. In Fig. 7, we showλ κ θ 's in the recoil frame calculated using the two frameworks, CGC+NRQCD (left plot) and NLO pQCD+NRQCD (right plot). A characteristic property of the NLO pQCD+NRQCD computations is that the 3 S [1] 1 and 3 P [8] J channels giveλ κ θ with a sign opposite to that for the 3 S [8] 1 channel; this cancellation leads to λ θ ∼ 0. As can be seen in Fig. 7, our CGC+NRQCD results also have a similar behavior confirming our expectation that the latter includes key physics of the NLO pQCD+NRQCD framework. We should emphasize that the divergence seen in Fig. 7 for the 3 P

[8]
J channel is not physical, since the quantity we are plotting is merely illustrative and is not what is measured in the J/Ψ polarization studies. In fact, we see that the CGC+pQCD computation is more stable than the NLO pQCD+NRQCD computation over the entire kinematic region shown. Note that the cancellations seen at NLO is not present [19] for the LO pQCD+NRQCD computation. To summarize, the most important difference relative to the NLO pQCD computations is that our framework includes higher twist gluon saturation contributions that become comparable to the leading twist contributions at small x and low p T . However at high p T our framework reduces to k T factorization which has a significant regime of overlap in p T with higher order collinear pQCD computations.
The accuracy of the CGC+NRQCD computations should be significantly improved once NLO computations to the heavy quark-antiquark pair impact factor become available. There has been recent progress in this direction [74][75][76] but much work remains to be done. An interesting question is the validity of the eikonal approximation that is assumed in the computation; these may potentially impact polarization observables more than unpolarized quantities. While there have been some efforts towards computing non-eikonal corrections in the CGC framework [77][78][79], the application of these methods to quarkonium polarization is beyond the scope of this paper.
We should note further that the results presented in this paper are obtained for direct J/ψ production, whereas only prompt [16] or inclusive polarization data [14,15] are available. However our calculation is quite accurate because feeddown contributions to J/ψ production from decays of higher charmonium states and B-hadrons are smaller at the low p T 's that we consider [67,80]. Furthermore, an analysis within the collinear factorization pQCD+NRQCD framework suggests that these feeddown contributions do not change the result significantly [81]. Indeed, we have checked within the CGC framework itself that the feeddown corrections from χ cJ states has a small impact on J/ψ polarization parameters, and are within the estimated theoretical uncertainty band.
Finally, we note that the LDME set we are using was obtained from the fit to the data using SDCs calculated in the NLO collinear factorization framework [19]. An independent fit of LDMEs to the data using the SDCs computed in the CGC is possible; this program will however be most effective when the above mentioned NLO impact factors will become available. The good description of the data both for the yields [38] and for polarization variables suggests that such an exercise is very worthwhile and much needed. This is especially so because some of the LDME sets obtained in the NLO pQCD+NRQCD approach from being positive to negative as pT increases [19]. This sign change is the reason forλ   (for instance in [20] and [64]) have negative values for some of the LDMEs. Using these LDMEs in our approach would lead to significant discrepancies between theory and data. A concern with using cross-sections at high p T is that small variations in p T can lead to large uncertainties in the LDMEs; a robust framework at lower values of p T can therefore greatly improve the precision with which they are extracted.

Summary and outlook
The CGC effective field theory has by now been used to successfully compute a large number of final states at collider energies. In this paper, we extended the CGC+NRQCD approach [37] which was previously used to compute J/ψ and ψ(2S) yields in protonproton collisions [38] to address the J/ψ polarization puzzle. We have computed the three polarization coefficients λ θ , λ φ , λ θφ for proton-proton collisions in this approach and have obtained on the whole quite good results describing the J/ψ polarization measured by the LHCb and ALICE experiments in both the Collins-Soper and recoil frames.
In Section 4.3, we discussed the differences between our approach and that of collinearly factorized pQCD+NRQCD approaches, as well as some future refinements of the extant CGC+NRQCD computations. In particular, our results provide strong motivation to compute the NLO impact factor in the CGC+NRQCD EFT for quarkonium production in proton-proton and proton-nucleus collisions. The NLO impact factor results, when available, will allow for significant improvements in the accuracy of the extracted LDMEs; these at present differ considerably between different NLO pQCD+NRQCD analyses.
Further, the increasing experimental precision of collider data opens up the possibility of computing the polarizations of the higher charmonium states in the CGC+NRQCD framework . In particular, the good description of data obtained for ψ(2s) yields [38,48] suggests that this framework may also describe the polarization of this meson. Measurements of the χ cJ states [82] give access to the 3 P [1] J channel [83]. What more, the production of quarkonia containing b quark pairs have been analyzed, in particular the polarization of Υ(ns) [84]. The analysis of these states requires that we extend the CGC+NRQCD computations to resum large log(p T /M ) terms that appear in the perturbative computations [85][86][87][88].
Finally, an analysis of J/ψ polarization in high multiplicity events in proton-proton and proton-nucleus collisions is possible. First computations of the dependence of J/ψ yields in proton-proton and proton-nucleus collisions were performed in [40] showing good agreement with data. The extension of this work studying the multiplicity dependence of the J/ψ polarization coefficients is in progress and will be reported separately.

B Computation of the CGC+NRQCD SDC
In this section, we sketch the procedure of obtaining the functions Γ κ ij and G κ ij that contribute to the short distance coefficients of the spin density matrix, in a full analogy to that computed in [37,38] for the CGC+NRQCD unpolarized quarkonium cross-sections.

B.1 Amplitude
We start by writing the amplitude obtained within the dilute-dense CGC formalism for cc pair production in a state κ with fixed spin and orbital momentum projections S z , L z and momentum p (see Eq. (3.7) in [37]). This has the form, These functions represent the action of the covariant spin projectors [89,90] on the amplitudes T QQ (p, q, k 1⊥ , k ⊥ ), T g (p, q, k 1⊥ , k ⊥ ) for the process g(k 1 ) → Q p 2 + q Q p 2 − q with Wilson lines attached to quarks or gluons respectively [35]. Explicit expressions for the latter can be found in Eq. (2.9) of [37]. Finally, the quarkonium polarization vectors µ can be expressed in terms of the unit axis vectors X, Y, Z in the J/ψ's rest frame defined in section 2 and are given by The functions Γ κ ij and G κ ij appearing in Eqs. (3.3) and (3.4) are obtained by taking the modulus squared of the amplitude in Eq. (B.1). The averaging over the color source densities in the projectile and target results in the unintegrated gluon densities and dipole amplitudes respectively, as described in [37,38]. The rest of the expression, containing all the information on the helicities then condensed into One needs to evaluate the Dirac traces of the form Tr Π 00 T QQ q=0 , Tr Π µ T QQ q=0 , Tr [Π µ T g ]| q=0 and ∂ ∂q α Tr Π ν T QQ q=0 and then contract them with polarization vectors µ . We do not show the explicit results for these traces because the expressions are lengthy but will provide them upon request.

C Contributions from different states to λ coefficients
To obtain further insight into how different channels contribute to the polarization coefficients λ, we define the variables λ κ ; these quantities have only contributions from given channel κ in the numerators while the denominator have contributions from all channels: Note that we have denoted here dσ ij = κ dσ κ ij to ensure that the denominators in Eq. (C.1) are the same as those in Eq. (2.7). These definitions then clearly satisfy Our results for λ κ are shown in Figure 8. We first observe that for p T → 0 all λ κ φ and λ κ θφ vanish. This is as anticipated because in this limit azimuthal symmetry is restored and no φ-dependence in Eq. (2.6) is possible. Note that the 1 S   J channels which add up with the same sign in most cases. At higher p T , the 3 S  J states. This is especially striking for λ θ , where the 3 S [8] 1 state contributes with a transverse polarization while 3 P [8] J state is longitudinally polarization, with the summation of the two giving a nearly unpolarized result. As noted in the introduction to this paper, this phenomenon is very similar to that  Figure 8. Parameters λ κ θ (first row), λ κ φ (second row) and λ κ θφ (third row) defined in Eq. (C.1) plotted for the Collins-Soper frame (left column) and recoil frame (right column). Different colors of bands represents different channels. "SUM" (denoted as black band) represents sum of all channels, so is equal to λ θ , λ φ , λ θφ from Figure 2 (note the difference in scales between these plots and those from Figure 2). observed in the NLO collinear factorization framework. So albeit the CGC computation is computed at LO in the impact factor, our result suggests that it may contain key features of the dynamics of the NLO collinear factorization framework. This was also seen in the close matching of the yields of the unpolarized cross-section observed in [38].