Production and polarization of prompt $J/\psi$ in the improved color evaporation model using the $k_T$-factorization approach

We calculate the polarization of prompt $J/\psi$ production in the improved color evaporation model at leading order employing the $k_T$-factorization approach. In this paper, we present the polarization parameter $\lambda_\vartheta$ of prompt $J/\psi$ as a function of transverse momentum in $p+p$ and $p$ + A collisions to compare with data in the helicity, Collins-Soper and Gottfried-Jackson frames. We also present calculations of the charmonium production cross sections as a function of rapidity and transverse momentum. This is the first $p_T$-dependent calculation of charmonium polarization in the improved color evaporation model. We find agreement with both charmonium cross sections and polarization measurements.


I. INTRODUCTION
The production mechanism of quarkonium remains uncertain even more than 40 years after the discovery of the J/ψ. Nonrelativistic QCD (NRQCD) [1], the most widely employed model of quarkonium production encounters serious challenges in both the universality of the long distance matrix elements (LDMEs) and the prediction of quarkonium polarization [2]. The production cross sections in NRQCD, based on an expansion in the strong coupling constant and the QQ velocity [3], is factorized into hard and soft contributions and divided into different color and spin states, including color octet contributions. The LDMEs, which weight the contributions from each color and spin state, are fit to the data above some minimum transverse momentum, p T . These LDMEs, which are conjectured to be universal, fail to describe the yields and polarization simultaneously for p T cuts less than twice the mass of the quarkonium state [4,5]. They also depend on the collision system [6][7][8][9]. Moreover, the polarization predicted by NRQCD is sensitive to the p T cut. The η c p T distributions calculated with LDMEs obtained from J/ψ yields using heavy quark spin symmetry [10][11][12], generally overestimates the high p T LHCb η c results [13]. NRQCD also consistently underestimates the Υ(nS) cross section ratio for 8 to 7 TeV as a function of p T [14].
On the other hand, the color evaporation model (CEM) [17][18][19][20], which considers all QQ (Q = c, b) production regardless of the quark color, spin, and momentum, is able to predict both the total yields and the rapidity distributions with only a single normalization parameter per state [21]. Both the CEM and NRQCD can describe production yields rather well but spin-related measurements like the polarization are strong tests of production models. However, polarization is not the only test of models. The CEM was also used recently to calculate transverse single spin asymmetries in J/ψ production [15,16].
We have previously presented the first polarization results in the CEM [22], which only considered charmonium and bottomonium production in general, followed by polarization results of the prompt J/ψ and Υ(1S) [23]. The later also took the feed-down production into account using the recently-developed improved CEM (ICEM) [20]. However, those results were at LO assuming collinear factorization and were thus p T -independent. This paper serves as a continuation of the previous work by presenting a p T -dependent leading order (LO) ICEM calculation of the polarization in prompt J/ψ production using the k T -factorization approach. This is a p T -dependent result because the transverse momenta of the incoming gluons and their off-shell properties are not neglected in the k Tfactorization approach. Our calculation provides the first p T -dependent ICEM polarization result and represents a step toward a full NLO ICEM polarization result. We will begin to address the p T dependence at NLO in a later publication.
The polarization 4-vector is where k µ T = (0, k T , 0). In the traditional CEM, all quarkonium states are treated the same as QQ below the HH threshold. The invariant mass of the heavy quark-antiquark pair is restricted to be less than twice the mass of the lowest mass meson (H) that can be formed with the heavy quark as a constituent. The distributions for all quarkonium family members are assumed to be identical.
In the ICEM the invariant mass of the intermediate heavy quark-antiquark pair is constrained to be larger than the mass of produced quarkonium state, M Q , instead of twice the quark mass, 2m q , the lower limit in the traditional CEM [17,22]. Because the charmonium momentum and integration range depend on the mass of the state, the kinematic distributions of the charmonium states are no longer identical in the ICEM and, for example the ψ to J/ψ ratio depends on p T . Using the k T -factorization approach, in a p + p collision, the ICEM production cross section for a directly-produced quarkonium state Q is where the square of the heavy quark pair invariant mass isŝ while the square of the center-of-mass energy in the p + p collision is s. Here Φ(x, k T , µ 2 F ) is the unintegrated parton distribution function (uPDF) for a parton with momentum fraction x and transverse momentum k T interacting with factorization scale µ F . The angles φ 1,2 in Eq. (3) are between the k T 1,2 of the partons and the p T of the final state quarkonium Q. The parton-level cross section is σ(R + R → QQ). Finally, F Q is a universal factor for the directly-produced quarkonium state Q, and is independent of the projectile, target, and energy. In this approach, the cross section is where the sum k 1T is over the roots of k 2 The momentum fractions x 10 and x 20 are Here, φ is the relative azimuthal angle between two incident Reggeons (φ = φ 1 − φ 2 ) and p T is the transverse momentum of the produced QQ.
The cross section may also be defined in terms of x F instead of rapidity as, where x 10 and x 20 are now Thus the transverse momentum distribution dσ/dp T in the ICEM is The two expressions are equivalent when calculating the transverse momentum without any longitudinal kinematic cuts. Eq. (12) is used to compare to collider data with defined rapidity cuts while Eq. (13) is used to compare to fixed-target data with x F cuts. Similarly, the rapidity distribution dσ/dy in the ICEM is dσ dy = dp T dŝdφ d 4 σ dp T dydŝdφ .
We take the renormalization and factorization scales to be µ F = µ R = m T , where m T is the transverse mass of the QQ. We will study the effect of varying these scales on the p T distributions and the polarization.

II. POLARIZATION OF DIRECTLY PRODUCED QQ
We define the polarization axis (z axis) in the helicity frame where z HX is the flight direction of the quarkonium in the center of mass frame of the colliding beams, as shown in Fig. 1. In this section we outline the kinematics required to compute the polarized scattering cross sections in the helicity frame as well as the procedure to relate them to the polarized scattering cross sections in the Gottfried-Jackson frame [27] and the Collins-Soper frame [28].
In the lab frame, using Eqs. (1) and (2) the momenta of the initial state Reggeons can be written as with polarization vectors We then boost the momenta along the beam direction to the frame where the total momentum of the Reggeons along the beam direction, k 1z + k 2z , is zero The polarization vectors are unchanged. We then apply a rotation such that the three momentum of the sum k µ 1 + k µ 2 is aligned with a new z-axis We then boost to the quarkonium rest frame where In this frame (helicity frame), the momenta of the initial state Reggeons are The orientation of polarization axis (z axis) in the helicity frame is indicated by the dashed arrow. The proton arrows indicate the incoming beam directions. The polarization axis is defined to be the direction of the produced (Q) travels in the center-of-mass frame of the colliding beams. If the quarks in the QQ pair with total angular momentum J = 1, they can either have the same angular momentum along the z-axis, Jz, or opposite resulting in Jz = 0 (longitudinal) or Jz = 1 (transverse), respectively.
The scattering amplitude of the process R + R → QQ is related to that of g + g → QQ by [25,26] where µ (k) is defined in Eq. (2). Evaluating A µν (g+g → Q +Q) using the conventional Feynman rules of QCD, there are three gg → QQ Feynman diagrams to consider at O(α 2 s ). Each diagram includes a color factor C and a scattering amplitude A. The generic matrix element for the gluon fusion process can be written as [29] M gg = C gg,ŝ A gg,ŝ + C gg,t A gg,t + C gg,û A gg,û . (31) In terms of the Dirac spinors u and v, the individual amplitudes are Here g s is the gauge coupling, m c is the charm quark mass, represents the gluon polarization vectors, γ µ are the gamma matrices, k (k) is the momentum of the initial state light quark (antiquark) or gluon, and p (p) is the momentum of final state heavy quark (antiquark).
In the process of evaluating the scattering amplitudes, we take advantage of the fact, that at O(α 2 s ), the final state QQ is produced with no dependence on the azimuthal angle φ (and thus L z = 0) in a rotated frame (primed frame) where the z axis is defined as the direction of one of the incoming Reggeons. Since the Reggeons are head to head in this frame, the scattering amplitudes are independent of the azimuthal angle φ . We first rotate the initial state momenta p from the helicity frame to the primed frame by an Euler rotation: The scattering amplitudes in the primed frame for S = 1, sorted S z , are Theŝ-channel amplitudes are: Thet-channel amplitudes are: Finally, theû-channel amplitudes are where χ = 1 − 4m 2 c /ŝ. The final state total spin is determined from the heavy quarks helicities. Two helicity combinations that result in S z = 0 are added and normalized to give the contribution to the spin triplet state (S = 1) in Eqs. (36), (38), and (40).
In this primed frame, to extract the projection on a state with orbital-angular-momentum quantum number L, we obtain the corresponding Legendre component A L in the amplitudes by Having obtained the amplitudes for S = 1 with S z = 0, ±1, and L = 0, 1 with L z = 0, we calculate the amplitudes for J = 0, 1, 2. The amplitudes for J = 1, found by adding S = 1 and L = 0, are Employing angular momentum algebra, the amplitudes for J = 0, 1, 2, obtained by adding S = 1 and L = 1, are Using a Wigner representation of the inverse rotation defined in Eq. (35), TABLE I. The mass MQ, the feed-down contribution ratio cQ, and the squared feed-down transition Clebsch-Gordan coefficients S Jz Q for all quarkonium states contributing to prompt J/ψ production.
3.56 0.14 2/3 1/2 Next, the amplitudes sorted by final state J and J z are squared for calculations in the helicity frame. For calculations in the other frames, the unsquared amplitudes can be further rotated to the Collins-Soper (CS) or the Gottfried-Jackson (GJ) frame. In the CS frame, the zaxis is defined as the angle bisector of the angle between one proton beam and the opposite of the other proton beam. In the GJ frame, the z-axis is defined as the direction of the momentum of one of the two colliding proton beams.
The squared matrix elements, |M| 2 , are calculated for each J, J z combination. The color factors, C, are calculated from the SU(3) color algebra and are independent of final state angular momentum [29]. They are Finally, the total squared amplitudes for a given J, J z state, are then employed to calculate the partonic cross sections The pT dependence of inclusive J/ψ production at √ s = 7 TeV in the ICEM obtained by varying the renormalization scale (blue solid), the factorization scale in the range 0.5 < µF /mT < 2 (magenta dashed), and the renormalization scale in the range 0.5 < µR/mT < 2 (green dot-dashed). The LHCb data [34] assuming the J/ψ polarization is totally transverse, λ ϑ = +1 (red), and totally longitudinal, λ ϑ = −1 (blue), are shown. The LHCb data assuming λ ϑ = 0 lie between the red and blue points and are not shown. by integrating over solid anglê The sum of the polarized partonic cross sections results for each final state total angular momentum J, is equal to the unpolarized partonic cross section, Having computed the polarized QQ production cross section at the parton level, we then convolute the partonic cross sections with the unintegrated parton distribution functions (uPDFs) to obtain the hadron-level cross section σ as a function of p T using Eq. (12) or (13) and as a function of y using Eq. (14). The quarkonium masses which appear as the lower limit of the QQ invariant mass are listed in Table I. We employ the ccfm-JH-2013-set1 [30] uPDFs in this calculation.

III. POLARIZATION OF PROMPT J/ψ
We assume that the angular momentum of each directly-produced quarkonium state is unchanged by the transition from the parton level to the hadron level, consistent with the CEM expectation that the linear momentum is unchanged by hadronization. This is similar to the assumption made in NRQCD that once a cc is produced in a given spin state, it retains that spin state when it becomes a J/ψ.
We calculate the J z = 0, ±1 to unpolarized ratios for each directly produced quarkonium state Q that has a contribution to prompt J/ψ production: J/ψ, ψ(2S), χ c1 (1P) and χ c2 (1P). These ratios, R Jz Q , are then independent of F Q . We assume the feed-down production of J/ψ from the higher mass bound states follows the angular momentum algebra. Their contributions to the J z = 0 to unpolarized ratios of prompt J/ψ are added and weighed by the feed-down contribution ratios c Q [31], where S Jz Q is the transition probability from a given state Q produced in a J z state to a J/ψ with J z = 0 in a single decay. We assume two pions are emitted for S state feed down, ψ(2S) → J/ψππ, and a photon is emitted for a P state feed down, χ c → Jψγ. S Jz Q is then 1 (if J z = 0) or 0 (if J z = 1) for Q = ψ(2S) since the transition, ψ(2S) → J/ψππ, does not change the angular momentum of the quarkonium state. For directly produced J/ψ, S Jz Q is 1 for J z = 0 and 0 for J z = 1. The S Jz Q for the χ states are the squares of the Clebsch-Gordan coefficients for the feed-down production via χ → J/ψ + γ. The values of M Q , c Q , and S Jz Q for all quarkonium states contributing to prompt J/ψ production are collected in Table I.

IV. RESULTS
Although the matrix element in this calculation is LO in α s , by convoluting the polarized partonic cross sections with the transverse momentum dependent uPDFs using the k T -factorization approach, we can calculate the yield as well as the polarization parameter λ ϑ as a function of p T . The full NLO polarization including qq and (q + q)g contributions, requiring us to go to O(α 3 s ), will be discussed in a future publication.
The traditional CEM can describe the unpolarized yield of charm and J/ψ production at both LO and NLO assuming collinear factorization [21,33]. The ICEM can also describe the ψ(2S) to J/ψ ratio at NLO while, in the traditional CEM, this ratio is independent of p T [20]. Since this is the first calculation in the ICEM using the k T -factorization approach, it is important to check if the unpolarized yield is also in agreement with the data.
In the remainder of this section, we first present how our approach describes the transverse momentum and rapidity distribution of the charmonium states in collider experiments. We then discuss the transverse momentum and rapidity dependence of the polarization parameter λ ϑ for prompt J/ψ production and direct production of quarkonium states that contribute to the feeddown production. We compare our results to the polarization measured in fixed-target experiments as well as collider experiments in the helicity, Collins-Soper, and Gottfried-Jackson frames to discuss the frame dependence of the polarization parameter. Finally, we discuss the sensitivity of our results to the factorization and renormalization scales, the weight of each diagram, and the feed-down ratios considered. In our calculations, we construct the uncertainty bands by varying the charm quark mass, around its base value of 1.27 GeV in the interval 1.2 < m c < 1.5 GeV, and the renormalization scale around its base value of m T in the interval 0.5 < µ R /m T < 2 while keeping the factorization scale fixed at µ F = m T . The total uncertainty band is constructed by adding the two uncertainties in quadrature.

A. Unpolarized charmonium production
In this section, we present the p T and rapidity distributions of charmonium states in our approach. In the spirit of the traditional CEM, F Q in Eq.(3) has to be independent of the projectile, target, and energy for each quarkonium state Q. Even though the focus of this paper is on polarization, which is F Q independent, the unpolarized yield in the ICEM using the k T -factorization approach was not considered before. Therefore, it is important to first confirm that this approach can indeed describe the charmonium yields as a function of p T and rapidity before discussing polarization predictions. We first obtain F J/ψ and F ψ(2S) by comparing our results with the experimental data measured by the LHCb Collaboration and the CDF Collaboration respectively. Using the same F J/ψ and F ψ(2S) , we compare our results with the experimental data measured at CDF and ALICE. We can only obtain F χc1 and F χc2 for the χ c states by comparing the unpolarized yield with the data measured by the AT-LAS Collaboration at √ s = 7 TeV because these are the only measurements. We instead give predictions of χ c1 and χ c1 production at √ s = 13 TeV. We also compare The direct χc1 (left) and χc2 (middle) pT dependence multiplied by the branching ratios for χc → J/ψγ and for J/ψ → µ + µ − , and the ratio of χc2 to χc1 (right) at √ s = 7 TeV (top panels) and at √ s = 13 TeV (bottom panels) in the ICEM with combined mass and renormalization scale uncertainties. The ATLAS data for prompt χc production are also shown [41]. and predict the ratio of χ c2 to χ c1 at √ s = 7 TeV and √ s = 13 TeV. Note that we cannot expect that our LO values of F Q to be equal to those found for Jψ and ψ(2S) in Ref. [20]. Those calculations are NLO in the total cross section assuming collinear factorization, and include the qq and (q +q)g channels where the contribution of the later is non-negligible.   [34] are shown as in Fig. 2. The LHCb data assuming λ ϑ = 0 are not shown.

J/ψ pT distribution
We first discuss why we fix the factorization scale at µ F = m T instead of including a factor of two variation, as usual in most other approaches. In Fig. 2, we show the p T distributions of inclusive J/ψ production at √ s = 7 TeV found by fixing m c = 1.27 GeV, and varying the factorization scale over the range 0.5 < µ F /m T < 2 and the renormalization scale over the range 0.5 < µ R /m T < 2 separately. We also fix µ F /m T = µ R /m T = 1 and vary the charm quark mass over the range 1.2 < m c < 1.5 GeV. The direct production cross section is calculated using Eq. (12) by integrating the pair invariant mass from M J/ψ to 2m D 0 (m D 0 = 1.86 GeV) over the rapidity range 2.0 < y < 4.5. We assume the direct production is a constant fraction, 0.62 of the inclusive production [31]. We then are able to compare the inclusive p T distribution in the ICEM with the LHCb data [34]. The result has a significant dependence on the factorization scale for p T > 5 GeV. This is because the uPDFs have a sharp cutoff for k T > µ F and are thus very sensitive to the chosen factorization scale. The yield varies more as p T approaches m T at high p T . At low p T , m T ∼ M Q and the cross section is independent of the factorization scale since k T << µ F . At moderate p T , the variation with µ F is similar to or smaller than that due to the charm quark mass. At p T ∼ 10 GeV, m T ∼ p T . Thus the lower limit on the factorization scale, m T /2, is on the order of k T and the yield drops off at this cutoff limit, while the upper limit on the factorization scale, 2m T , is still greater than k T , enhancing the yield. Since at LO, only the QQ pair carries the transverse momentum, the  predictive power of the yield is limited by the uPDFs. Therefore, to construct a meaningful uncertainty band, we fix the factorization scale at µ F = m T . As we push toward the limit of the k T -factorization approach with uPDFs at high p T at LO, we can only improve the high p T limit by a full NLO calculation.
After fixing the factorization scale, the variation in renormalization scale then gives the largest uncertainty, followed by the variation in charm mass. When µ R is reduced, the strong coupling constant is larger, increasing the yield. On the other hand, when m c is reduced, the yield increases. In the remainder of this section, we present our results by adding the uncertainties due to variations of the charm mass and renormalization scale in quadrature.
The inclusive J/ψ p T distribution at √ s = 7 TeV with combined uncertainty is shown in Fig. 3. The ICEM result has a peak at p T ∼2 GeV, in agreement with the experimental results but slightly overestimates the data at high p T . The ICEM p T distribution is within reasonable agreement with the data for all p T . The experimental prompt production cross section depends on the polarization of J/ψ since the polarization affects the acceptance and reconstruction efficiencies. LHCb checked the yields for the three polarization assumptions: λ ϑ = −1, 0, +1. The experimental p T distribution for all polarization assumptions is within the uncertainty band constructed in the ICEM. By matching to the experimental unpolarized yield λ ϑ = 0, we find that the ICEM can describe the J/ψ p T distribution with F J/ψ = 0.0216. This is the fraction of cc pairs produced in the invariant mass range from M J/ψ to 2m D 0 that result in direct J/ψ, defined in Eq. (3). We test the universality of F J/ψ by comparing the inclusive J/ψ p T distribution in the ICEM at √ s = 1.96 TeV and |y| < 0.6 with the CDF data [35] in Fig. 3. We again assume the direct production takes a constant fraction of 0.62 of the inclusive production [31] to obtain the inclusive J/ψ cross section. The ICEM results slightly overshoot the data at high p T because both the direct and non-prompt contributions to J/ψ production are p T dependent [34,36]. The direct-to-prompt J/ψ ratio decreases as p T grows and the contribution from b decay to inclusive production is measured to be larger at high p T than at low p T . Combining the effects of both, using a constant direct-to-inclusive ratio of 0.62 gives an overestimate of the yields at high p T . The calculated cross section differs from the measurements more as p T increases. We note that if we fix F J/ψ from the CDF data alone, it agrees within 1.5% of that extracted from comparison to the LHCb data.

ψ(2S) pT distribution
The inclusive ψ(2S) p T distribution at √ s = 1.96 TeV is shown in Fig. 5. Here, the direct production cross section is calculated using Eq. (12) by integrating the pair invariant mass from M ψ(2S) to 2m D 0 over the rapidity range |y| < 0.6. We assume the direct production is the same as the prompt production as there are no quarkonium states that feed down to ψ(2S) since its mass is just below 2m D 0 . Therefore, we compare the p T -integrated yield of direct ψ(2S) with the CDF measurement [37]. We find F ψ(2S) = 0.117. We note that F ψ(2S) > F J/ψ , primarily because the mass range is much smaller for ψ(2S) than J/ψ. In the traditional CEM, F ψ(2S) is smaller than F J/ψ because the integration over the pair invariant mass is the same for both J/ψ and ψ(2S). We add the contribution from non-prompt production reported by the CDF Collaboration to our prompt production yield to give the inclusive ψ(2S) yield shown in Fig. 5. We find agreement with the data within the combined uncertainty band constructed by varying the charm mass and the renormal-  ization scale in the ICEM.

χc1 and χc2 pT distribution
We now turn to the p T dependence of χ c production. The p T distributions of direct χ c1 , direct χ c2 , and the ratio of χ c2 to χ c1 at √ s = 7 TeV and 13 TeV are presented in Fig. 6. The direct production is calculated using Eq. (12) by integrating the pair invariant mass from M χc to 2m D 0 (m D 0 = 1.86 GeV) over the rapidity range |y| < 0.75. We assume the prompt production of χ c is approximately the same as the direct production. Thus, by comparing the direct χ c1 and χ c2 yields in the ICEM with the experimental yield of prompt χ c1 and χ c2 at √ s = 7 TeV measured by the ATLAS Collaboration [41], we obtain F χc1 = 0.180 and F χc2 = 0.20. As is the case for F ψ(2S) and F J/ψ , F χc2 > F χc1 is because the integration range over the pair invariant mass is smaller for χ c2 than for χ c1 . In the tradition CEM, F χc2 is smaller than F χc1 . The direct production in the ICEM describes prompt production of both χ c1 and χ c2 at √ s = 7 TeV within the uncertainty bands constructed by varying the charm quark mass and renormalization scale. The ratio of the cross sections is also described by the ICEM. We calculate the χ c2 to χ c1 ratio to be ∼ 0.5, almost independent of p T . The ratios disagree with a recent NRQCD calculation [40], which the ratio decreases as p T increases and is above the data. We assume that p T χc ≈ p T J/ψ , not unreasonable since the mass difference is ∼500 MeV and the decay photon is soft. We anticipate the direct χ c1 and χ c2 yields will be increased by 51% (at p T = 10 GeV) to 120% (at p T = 30 GeV) when √ s is increased from 7 TeV to 13 TeV. However, the ratio of χ c2 to χ c1 should remain approximately the same.

Prompt J/ψ pT distribution
After fixing F J/ψ , F ψ(2S) , F χc1 and F χc2 , we calculate the prompt J/ψ p T distribution at √ s = 7 TeV in the rapidity range 2.0 < y < 4.5 using the direct J/ψ, ψ(2S), χ c1 and χ c2 yields and their branching ratios to J/ψ. The prompt J/ψ p T distribution is shown in Fig. 7. The ICEM p T distribution describes the data for most p T but overshoots the data slightly at the highest p T bin. The ICEM p T distribution is within reasonable agreement with the data for all p T . We extract the p T dependent feed-down ratios c ψ 's by taking the direct to prompt ratio in this distribution. We find the feed-down ratios are very similar to those listed in Table I. Additionally, we find c J/ψ decreases as p T increases, in agreement with Ref. [36].

J/ψ rapidity distribution
We now turn to the rapidity dependence of J/ψ production. The rapidity distribution of inclusive of J/ψ at √ s = 7 TeV is shown in Fig. 8. The direct production is calculated using Eq. (14) by integrating over the p T range 0 < p T < 7 GeV (|y| < 0.9) and 0 < p T < 8 GeV (2.5 < y < 4). We again assume the direct production is a constant 62% [31] of the inclusive production. We use the same F J/ψ again to compare the rapidity distribution in the ICEM with the measurement made by the ALICE Collaboration [38]. The difference in the integrated p T range has a negligible on the rapidity distribution because the p T dependence has already dropped by an order of magnitude by p T ∼ 7 − 8 GeV. We find the ICEM can describe the ALICE rapidity distribution at √ s = 7 TeV using the F J/ψ obtained at the same energy by LHCb in the forward rapidity region.

ψ(2S) rapidity distribution
The rapidity distribution of direct ψ(2S) at √ s = 7 TeV is shown in Fig. (9). Here, the rapidity distribution is calculated in the interval p T < 12 GeV at forward rapidity (2.5 < y < 4). We use the same F ψ(2S) compare with inclusive ψ(2S) data from ALICE [39]. While the lower bound of our uncertainty band should still be lower than the data when the contribution from B decays are added, our baseline should slightly overshoot the inclusive ψ(2S) data. Our results also agree with the direct ψ(2S) rapidity distribution from a recent NRQCD calculation at LO using the k T -factorization approach [40].
Here, we present the p T dependence of the polarization parameter λ ϑ in p + p and p+A collisions. Because the polarization parameter is defined as the ratio of polarized to unpolarized cross sections in Eq. (71) and these cross sections depend on µ R and µ F in the same way, the polarization parameter is independent of the scale choice. However, the amplitudes themselves are mass dependent so that the polarized to unpolarized ratio in λ ϑ depends on the charm quark mass. Thus the only uncertainty on λ ϑ in our calculation is due to the variation of m c in the range 1.2 < m c < 1.5 GeV. In this section, the uncertainty band is only due to the mass variation and therefore the uncertainty is reduced relative to the yield calculations.

Charmonium polarization in p+A collisions at fixed-target energies
The polarization results for prompt production of J/ψ at √ s N N = 41.6 GeV are shown in Figs. 10 and 11. Although the HERA-B data are taken on nuclear targets, C and W, and there are known nuclear modifications of the parton densities in the nucleus, λ ϑ is independent of any modification. This is because the ratios of the polarized to unpolarized cross sections are in the same kinematic acceptance and any nuclear effects cancel in the ratio. Thus there is no difference in polarization between the two target nuclei. We compare our results with the C and W combined data measured by the HERA-B Collaboration in the region −0.34 < x F < 0.14 [43]. Prompt J/ψ polarization in the ICEM is close to unpolarized in both the CS and GJ frames for p T < 5 GeV. At p T = 0, the two z-axes z CS and z GJ , are in the same direction. Thus the polarization is the same in that limit. As p T increases, the two axes depart from each other. Thus the polarization is slightly less longitudinal in the GJ frame than in the CS frame. This behavior is also consistent with the experimental data showing that the J/ψ polarization at very low p T is not affected by switching from the CS frame to the GJ frame. At higher p T the polarization is slightly less longitudinal in the GJ frame than in the CS frame. The ICEM results are in fair agreement with the experimental data except at the lowest p T .

Charmonium polarization in p+p(p) collisions
We present the polarization parameters for prompt J/ψ in p+p collisions at √ s = 200 GeV in Fig. 12.
We compare our results with the data from the STAR Collaboration in the region |y| < 0.5 [44] in the he-  licity frame. The ICEM polarization of prompt J/ψ in the helicity frame is slightly transverse at low p T (p T < M J/ψ ). The result becomes unpolarized at moderate p T (M J/ψ < p T < 2M J/ψ ) before changing to slightly transverse at high p T . The ICEM polarization agrees fairly well with the data at small and moderate p T for inclusive J/ψ polarization at STAR.
We also compared the polarization parameters for prompt J/ψ in p+p collisions at √ s = 1.96 TeV with the data measured by the CDF Collaboration in the region |y| < 0.6 [45] in the helicity frame, shown in Fig 13. The ICEM prompt J/ψ polarization does not depend strongly on √ s or whether the collision is p+p or p+p. We find the trend in the p T dependence of the polarization is the same. At high p T , the prompt J/ψ polarization measured by the CDF Collaboration is slightly longitudinal to un- polarized while the ICEM polarization is slightly transverse. The polarization predicted by NRQCD also shows a similar behavior at this energy [46]. However, NRQCD predicts a stronger transverse polarization (λ ϑ ∼0.6) than ICEM in the high p T limit.
C. Rapidity dependence of λ ϑ Next we turn to the rapidity dependence of λ ϑ . We calculate the prompt J/ψ polarization in the helicity frame for p + p collisions at √ s = 7 TeV in two rapidity ranges, |y| < 0.6 and 0.6 < |y| < 1.2, shown in Figs. 14 and 15 respectively. We compare our results to the experimental data from the CMS Collaboration [47]. There is no difference in the polarization of prompt J/ψ in these two rapidity regions in the ICEM. In the ICEM, the polarization parameter λ ϑ of prompt J/ψ production increases very slowly in the high p T limit and reaches λ ϑ ∼ 0.12 at p T = 70 GeV. The ICEM polarization agrees with the the experimental results at central rapidity within uncertainty except the data in the 30 < p T < 35 bin. The experiment reports the polarization is less transverse in the forward rapidity region. Our results in the ICEM still agrees with the data even though the calculated polarization does not depend on rapidity in this range at 7 TeV.
We also do not observe variations in the polarization parameter λ ϑ at √ s = 7 TeV in the region of y < 4 using the same kinematics cut compared to the ALICE yield measurement in Fig. 8. We present the polarization as a function of rapidity in Fig. 16. The polarization parameter of prompt J/ψ for the p T -integrated results is λ ϑ = 0.26 ± 0.02.

D. Frame dependence of λ ϑ
We now turn to the frame dependence of our 7 TeV results. We calculate the polarization parameter in p + p collisions at √ s = 7 TeV in both the helicity frame and the Collins-Soper frame, shown in Figs. 17 and 18 respectively. The polarization in the Collins-Soper frame is opposite to that in the helicity frame in the ICEM. We expect this because, in these kinematics, at order α 2 s , the polarization axis in the Collins-Soper frame is always perpendicular to that in the helicity frame. Therefore, at low p T , where the J/ψ is predicted to be slightly transverse in the helicity frame, it is predicted to be slightly longitudinal in the Collins-Soper frame. Whereas, at moderate p T , where the J/ψ is predicted to be unpolarized, it is also predicted to be unpolarized in the Collins-Soper frame. This behavior, however, is not measured experimentally. As we compare our results with the ALICE data [48], the ICEM polarization agrees with the data in the Collins-Soper frame but does not agree with the data in the helicity frame, especially at low p T where the frame dependence is most significant.
We find similar results by comparing to the LHCb data in the Collins-Soper frame [49], show in in Figs. 19 and 20: the polarization in the ICEM agrees with the data in the Collins-Soper frame but not in the helicity frame. We expect that the difference in agreement of the calculations in different frames with the data may be resolved with a full α 3 s calculation of the ICEM cross section. Finally, we note that at low p T the polarization in the Gottfried-Jackson frame is similar to that in the Collins-Soper frame, as shown in Figs. 10 and 11 for fixed-target energies. However at high p T , the polarization in the Gottfried-Jackson frame is similar to that in the helicity frame. The differences are due to the definition of the polarization axes in the quarkonium rest frame. When p T << m T , the angle between the polarization axis in the Gottfried-Jackson frame and that in the Collins-Soper frame is small. As p T increases, the polarization axis in the Gottfried-Jackson frame becomes collinear with that in the helicity frame. Therefore, the polarization calculated in the Gottfried-Jackson frame is opposite to that in the helicity frame at low p T , and thus similar to that in the Collins-Soper frame. But as p T increases, the polarization in the Gottfried-Jackson frame should asymptotically approach the polarization in the helicity frame.

E. Sensitivity to scales and quark mass
We have already discussed the sensitivity of the charmonium yields to the factorization and the renormalization scales in section IV A 1. Here we note that the longitudinal to unpolarized fraction R Jz=0 J/ψ used in the calculation of λ ϑ , is insensitive to scale variations because the longitudinal and transverse change similarly as the scales are varied. Therefore, the polarization parameter λ ϑ for prompt J/ψ is independent of the scales for all energies considered. Similarly, while the unpolarized χ c1 and χ c2 cross section vary appreciably with the scale choice, the  χ c2 to χ c1 ratio is also independent of scales.
While the scale variations affect the polarized and unpolarized cross sections the same way, making λ ϑ scale independent, the J z components of the polarized cross section depend differently on quark mass. When p T ≤ M Q , the longitudinally polarized partonic cross section decreases faster with increasing m c than the transversely polarized partonic cross section in the helicity frame. Thus increasing the charm mass results in more transverse polarization. When p T > M Q , the longitudinally polarized partonic cross section decreases more slowly with increasing m c than the transversely polarized partonic cross section. Thus, here increasing the charm mass results in more longitudinal polarization. As p T ŝ, λ ϑ becomes insensitive to m c . Thus the uncertainty in λ θ is narrower. The pT dependence of the polarization parameter λ ϑ for prompt J/ψ production at √ s = 1.96 TeV in the ICEM with mass uncertainty when theŝ-channel contribution is excluded. The CDF data are also shown [45].

F. Sensitivity to feed-down ratios
We have tested the sensitivity of our results to the feeddown ratios used in our calculations [31]. Since prompt J/ψ production is dominated by direct J/ψ, we vary the feed-down ratio by changing the relative contribution of direct J/ψ and decays from excited states. Thus when the direct fraction, c J/ψ , increases, all other c ψ decrease and vice versa. Using the base values of c ψ in Table I and the reported uncertainty, we vary the feed-down ratios as given in Table II. Since the polarization of prompt J/ψ production does not vary at central rapidity, we study changes in the polarization by varying the feed-down ratios at y = 0. The p T -integrated polarization parameter for prompt J/ψ production at √ s = 7 TeV at y = 0 varies by 0.04 from 0.26 in the helicity frame. This variation is similar to that due to the charm quark mass and renormalization scale variations combined.

G. Sensitivity to diagram weights
We have tested the sensitivity of our results to diagram weights. As shown in Ref. [26], theŝ-channel diagram dominates color-octet production at high p T . Turning off the contribution from this diagram by setting A gg,ŝ = 0 in Eq. (31) makes a significant difference in polarization as well as the uncertainty band in the high p T limit. At 5 GeV, turning off the contribution from theŝ-channel diagram reduces the cross section by 70%. The difference is larger at higher p T . Thus the polarization is more sensitive to charm mass and gives a wider uncertainty band. The polarization parameter at √ s = 1.96 TeV in the rapidity region |y| < 0.6 in the helicity frame in this case is shown in Fig. 21. The polarization at low p T is more transverse compared to Fig. 13. Instead of becoming slightly transverse at high p T , prompt J/ψ production will remain approximately unpolarized with λ ϑ = +0.14 +0.04 −0.14 in the helicity frame when theŝ-channel amplitude is completely turned off.

V. CONCLUSIONS
We have presented the transverse momentum and rapidity dependence of the charmonium cross section as well as the the polarization of prompt J/ψ production in p + p and p+A collisions in the improved color evaporation model in the k T -factorization approach. We compare the p T dependence to data at both fixed-target energies and collider energies. We also present χ c predictions as a function of p T at √ s = 13 TeV. We find prompt J/ψ production to be unpolarized at moderate p T and slightly transverse in the high p T limit in the helicity frame. We do not observe any rapidity dependence in the polarization in the ranges considered. We report the p Tintegrated polarization parameter for prompt J/ψ production at √ s = 7 TeV to be λ ϑ = 0.26 ± 0.02 at y = 0 in the helicity frame. We will study the p T dependence of bottomonium states in this approach in a future publication. Since our calculation of the matrix elements is leading order in α s , the high p T cross section varies strongly with the choice of factorization scale due to the limitations on the uPDFs as x increases. We expect improvements at high p T when we calculate the cross section to O(α 3 s ) in a future publication.