Production and polarization of direct $J/\psi$ to ${\mathcal O}(\alpha_s^3)$ in the improved color evaporation model in collinear factorization

We calculate the production and polarization of direct $J/\psi$ in the improved color evaporation model at $\mathcal{O}(\alpha_s^3)$ in the collinear factorization approach. We present the first calculation of polarization parameters $\lambda_\vartheta$, $\lambda_{\phi}$, and $\lambda_{\vartheta \phi}$ in the helicity and the Collins-Soper frames, as well as the frame-invariant polarization parameter $\tilde{\lambda}$ as a function of transverse momentum. We find agreement with both $J/\psi$ cross sections and the invariant polarization parameters at small and moderate $p_T$.


I. INTRODUCTION
Quarkonium production is important to understand both the long-and short-distance aspects of QCD. Both the perturbative and nonperturbative natures of QCD are needed to model the production of quarkonium from heavy quark production in hard processes to the hadronization of the final state. Nonrelativistic QCD (NRQCD) [1], the most commonly employed model of quarkonium production cannot describe the J/ψ production and polarization while respecting the universality of the long-distance matrix elements for p T cuts less than twice the mass of the quarkonium state [2,3]. It also has difficulty describing the J/ψ polarization and the LHCb η c production [4,5] while using heavy quark spin symmetry [6][7][8]. On the other hand, the color evaporation model (CEM) [9][10][11] and the improved CEM (ICEM) [12] have only been employed extensively to hadroproduction on S-state quarkonia.
Our previous charmonium and bottomonium polarization calculations performed in the k T -factorization approach [13,14] describe both polarization and production at most p T . However, the p T dependence of the production calculation has a strong dependence on the factorization scale chosen. Also, the discrepancies between the ICEM polarization calculation and the measured data are not visualized in a frame-independent way. We address these issues in this paper by performing a polarized production calculation in the collinear factorization approach and also computing a frame-invariant polarization parameter to compare with the data. In this calculation, only the production and polarization of direct J/ψ is presented. We will address the effects of feed-down production on J/ψ in a later publication.
In this paper, we present both the yield and the polarization parameters of direct J/ψ production as a function of p T in the ICEM [12] using the collinear factorization approach. The polarized cross section calculation is presented in Sec. II. The results, along with comparison of unpolarized p T distributions and polarization parameters to data, are shown in Sec. III. Our conclusions are presented in Sec. IV.

II. POLARIZED CROSS SECTION
The ICEM assumes the J/ψ production cross section takes a constant fraction of the open cc cross section with invariant mass above the mass of the J/ψ but below the hadron threshold, the DD pair mass. A distinction is also made between the cc momentum and the J/ψ momentum in the ICEM compared to the traditional CEM. The unpolarized direct J/ψ production cross section in p + p collision in the ICEM is given by where i and j are q,q and g such that ij = qq, qg,qg or gg, F J/ψ is a universal factor at fixed order for direct J/ψ production in the ICEM independent of projectile and energy, x is the momentum fraction of the parton, and f (x, µ F ) is the parton distribution function (PDF) for a parton in the proton as a function of x and the factorization scale µ F . Finally,σ ij→cc+k are the partonlevel cross sections for initial states ij to produce a cc pair with a light final-state parton k. The invariant mass of the cc pair, M , is integrated from the physical mass of J/ψ (M J/ψ = 3.10 GeV) to two times the mass of the D 0 hadron (2m D 0 = 3.72 GeV). Because the O(α 3 s ) contribution diverges when the light parton is soft, in order to describe the p T distribution at low J/ψ p T , the initial-state partons are each given a small transverse momentum, k T , kick of k 2 are then multiplied by the Gaussian function g(k T ), assuming the x and k T dependences completely factorize. The same Gaussian smearing is applied in Refs. [12,15,16]. Note that in the traditional CEM, the lower invariant mass threshold for all charmonium states is set to the production threshold, which makes the kinematic distributions of the charmonium states identical except for the choice of F Q . The distinction between the J/ψ and cc momenta also helps describe the p T distributions at high p T . We consider diagrams with the projection operators, , applied to the cc [17,18] to calculate the partonic cross sections. We denote the momenta of i, j, c,c, and k in the partonic process i + j → c +c + k are denoted as k 1 , k 2 , p c , pc, and k 3 , respectively, where k is the emitted parton, with 1 (S 1z ), 2 (S 2z ), and 3 (S 3z ) denoting the polarization of the light partons. When calculating the 2 → 3 cross section, we transformed the momenta of the charm quark (p c ) and the anticharm quark (pc) into the momentum of the proto-J/ψ (p ψ ) and the relative momentum of the heavy quarks (k r ), and denote the polarization vector of the proto-J/ψ as ψ (J z ), where J z is the spin projection onto the polarization axis. Instead of taking the limit k r → 0, we note that, since the mass of the proto-J/ψ is integrated from the physical mass of J/ψ to the hadronic threshold, the relative momentum k r depends on the mass of the proto-J/ψ. We factorized the amplitudes into a product of the color factors, C, and the colorless amplitude A, so that M = CA. The color factors, C, in the squared amplitudes are calculated separately by summing over all colors and averaging over the initial-state colors.
There are 16 diagrams for the gg → ccg process. The factorized amplitudes, A, arranged by the number of three-gluon vertices, are The diagram with a four-gluon vertex is factorized according to with A gg4,2 = ig 3 s tr γ ν / * There are five gq → ccq diagrams, which, written in terms of Dirac spinors, are There are five gq → ccq diagrams, obtained by replacing the spinors in the above five gq → ccq diagrams: Finally, five diagrams contribute to qq → ccg, We assume that the angular momentum of the proto-J/ψ is unchanged by the transition from the parton level to the hadron level. We then convolute the partonic cross sections with the CT14 PDFs [19] in the domain where p ψ · k = 0. We restrict the partonic cross section calculations within the perturbative domain by introducing a regularization parameter such that all propagators are at a minimum distance of Q 2 reg = M 2 from their poles, as employed in Ref. [17]. We take the factorization and renormalizaton scales to be µ F /m T = 2.1 +2.55 −0.85 and µ F /m T = 1.6 +0.11 −0.12 respectively, where m T is the transverse mass of the charm quark produced . We also vary the charm quark mass around 1.27 ± 0.09 GeV. These variations were determined in Ref. [15], in which the uncertainties on the total charm cross section were considered.

III. POLARIZED PRODUCTION OF DIRECT J/ψ
We factor the polarization vector, ψ (J z ), from the unsquared amplitudes for all sub processes, giving us the form for each subprocess denoted by the initial states, n = gg, gq, gq, qq. The polarization vectors for J z = 0, ±1 in the rest frame of the proto-J/ψ are using the convention that the fourth component is the z component. While the unpolarized cross section does not depend on the choice of z axis, the polarized cross sections does depend on the orientation of the z axis. In this calculation, the y axis is chosen to be the normal vector of the plane formed by the two beams with momenta P 1 and P 2 ,ŷ In the helicity frame, the z HX axis is the flight direction of the cc pair in the center of mass of the colliding beams.
In the Collins-Soper frame [20], the z CS axis is the angle bisector between one beam and the opposite direction of the other beam. The x axis is then determined by the right-handed convention. We compute the polarized cross section matrix element, M n , in the rest frame of the cc pair by first taking the product of the unsquared amplitude with polarization vector of J z = i z and the unsquared amplitude with polarization vector of J z = j z in each subprocess (n), then adding them, and finally calculating the components of the polarized cross section matrix according to Eq. (1), where i z , j z = {−1, 0, +1} and the integral is over all variables explicitly shown in Eq. (1) as well as the Lorentzinvariant phase space in 2 → 3 scatterings. The unpolarized cross section is the trace of the polarized cross section matrix The polarization parameters are calculated using the matrix elements. The polar anisotropy (λ ϑ ), the azimuthal anisotropy (λ ϕ ), and polar-azimuthal correlation (λ ϑϕ ) are given by [21] These parameters depend on the frame (helicity or Collins-Soper) in which they are calculated and measured. Since the angular distribution itself is rotationally invariant, there are ways to construct invariant polarization parameters from Eqs. (32)-(34). One of the combinations to form a frame-invariant polarization parameter (λ) is [21] The choice ofλ is the same as the polar anisotropy parameter (λ ϑ ) in a frame where the distribution is azimuthally isotropic (λ ϕ = 0). We can remove the frameinduced kinematic dependencies when comparing theoretical predictions to data by also considering the frameinvariant polarization parameter,λ.

IV. RESULTS
We first show how our approach describes the transverse-momentum-distributions of J/ψ compared to ALICE [22] and ATLAS [23] measurements at √ s = 7 TeV and compare our results with previous calculations in the ICEM. We then discuss the transverse momentum and rapidity dependences of the frame-dependent polarization parameters λ ϑ , λ ϕ , and λ ϑϕ as well as the frameinvariant polarization parameterλ compared to the data measured by the LHCb Collaboration [24], the ALICE Collaboration [25], and the CMS Collaboration [26]. In our calculations, we consider theoretical uncertainties by varying the charm quark mass, the renormalization scale, and the factorization scale as discussed in Sec. II. We also estimate the uncertainty on the fit parameter, δF Jψ . The total uncertainty band is then constructed by adding the mass and scale uncertainties in quadrature and varying F J/ψ from F J/ψ − δF Jψ to F J/ψ + δF Jψ .
A. Unpolarized J/ψ pT distribution We calculate the p T distribution of direct J/ψ production at √ s = 7 TeV with |y| < 0.9. We assume direct production is a constant fraction, 0.62, of the inclusive production [27] to obtain the inclusive J/ψ p T distribution. We compare our ICEM inclusive J/ψ p T distribution with the data measured by the ALICE Collaboration [22]. The comparison is presented in Fig. 1. By comparing the total cross section for p T < 7 GeV, we find F J/ψ = 0.0363, consistent with previous CEM [15] and ICEM calculations [13]. We add the statistical and systematic uncertainties of the ALICE total cross section in quadrature. This results in a proportional uncertainty of 17% on the total cross section. Using the same proportional uncertainty, we estimate the uncertainty on F J/ψ , δF Jψ to be ±0.0062. We note that since this calculation is done using collinear factorization, the variation in the factorization scale does not result in a large uncertainty band as seen in our previous calculation using the k T -factorization approach [13]. Overall, we have good agreement with the data over the p T range measured.
We compare the ICEM inclusive J/ψ p T distributions at √ s = 7 TeV in the forward rapidity region covered by the ALICE detector, 2.5 < y < 4 [22], and in the  [22]. The ALICE data are measured while assuming J/ψ production is unpolarized. The pT dependence of inclusive J/ψ production at √ s = 7 TeV with 2.5 < y < 4 (blue region) and |y| < 2 in the ICEM (red region). The combined mass, renormalization scale, factorization scale, and F J/ψ uncertainties are shown in the band. They are compared to the ALICE data [22] (blue squares) and the ATLAS data [23] (red squares) respectively. central rapidity region covered by the ATLAS detector, |y| < 2 [22]. These comparisons are presented in Fig. 2. Both ALICE and ATLAS data are shown assuming J/ψ production is unpolarized. Since the detector acceptance depends on the polarization assumption, the ATLAS Collaboration reports a set of correction factors for each polarization assumptions. These factors ranges from 0.7 to 1.7 for the lowest p T bin, and from 0.90 to 1.06 for the highest p T bin. Except for the highest p T bin reported by the ATLAS Collaboration, we find our ICEM results in good agreement with the data over the p T range measured in both kinematic regions using the same F J/ψ found by comparing to the ALICE central rapidity data.
We compare the p T distribution of inclusive J/ψ production at √ s = 7 TeV with 2 < y < 4.5 in this calculation with that from previous calculations in the ICEM shown in Fig. 3. We compare these calculations with the data measured by the LHCb Collaboration [28]. Although this distribution is for unpolarized J/ψ, the calculation selects only cc with the same spin as the J/ψ. We thus refer to this calculation as polarized collinear ICEM. We compare this distribution with the polarized ICEM in the k T -factorization approach [13] and the unpolarized ICEM in the collinear factorization approach [12]. The former also selects cc with the same spin as the J/ψ but the latter is a spin-averaged calculation. The uncertainty band of the unpolarized collinear ICEM is constructed in the same way as the polarized collinear ICEM. The uncertainty band of the k T -factorized ICEM is constructed by varying the renormalization scale in the interval 0.5 < µ R /m T < 2 and varying the charm mass in the interval 1.2 < m c < 1.5. We find our polarized collinear ICEM p T distribution behaves very similarly to the unpolarized collinear ICEM as p T increases and agrees with the data and other calculations in the ICEM in general.
B. pT and y dependences of λ ϑ , λϕ, and λ ϑϕ We calculate the p T dependence of the framedependent polarization parameters λ ϑ , λ ϕ , and λ ϑϕ at √ s = 7 TeV in the helicity frame and in the Collins-Soper frame. We compare the polarization parameters at low and moderate p T (2 < p T < 15 GeV) with the data measured by the LHCb Collaboration [24] and the ALICE Collaboration [25], where the data are collected in rapidity ranges 2 < y < 4.5 and 2.5 < y < 4. The comparisons in the helicity frame and in the Collins-Soper frame are presented in Figs. 4 and 5, respectively. We then compare the polarization parameters at high p T (14 < p T < 70 GeV) with the data measured by the CMS Collaboration [25], where the data are collected in rapidity ranges |y| < 0.6 and 0.6 < |y| < 1.2. These comparisons in the helicity frame and in the Collins-Soper frame are presented in Figs. 6 and 7, respectively.
The polar anisotropy parameter (λ ϑ ) reflects the proportion of the J/ψ in each spin projection state, with λ ϑ = 1 referring to completely transverse production of J z = ±1 and λ ϑ = −1 referring to completely longitudinal production of J z = 0. At low p T , λ ϑ is close to zero in both the helicity frame and the Collins-Soper frame, indicating equal J/ψ yields in each spin projection state (J z = 0, ±1). However, as p T grows larger, the difference between the λ ϑ calculated in the two frames increases with p T . In the helicity frame, the transverse component (J z = ±1) falls off more slowly than the longitudinal √ s = 7 TeV in the polarized collinear ICEM (this calculation) (blue region), in the polarized ICEM using the kT -factorization [13] (magenta region), in the unpolarized collinear ICEM [12] (green region). They are compared to the LHCb data [28] assuming that the J/ψ polarization is totally transverse, λ ϑ = +1 (red squares), and totally longitudinal, λ ϑ = −1 (blue squares). The LHCb data assuming λ ϑ = 0 lie between the red and blue points and are not shown.  4. (a) The polar anisotropy parameter (λ ϑ ), (b) the azimuthal anisotropy parameter (λϕ), and (c) the polar-azimuthal correlation parameter (λ ϑϕ ) (c) in the helicity frame at √ s = 7 TeV in the ICEM. The combined mass, renormalization scale, and factorization scale uncertainties are shown in the band and compared to the LHCb data [24] (blue) and the ALICE data [25] (red).  [24] (blue) and the ALICE data [25] (red). component (J z = 0). As a result, λ ϑ becomes positive as p T increases. This is consistent with the (color glass condensate) CGC+NRQCD approach at low and moderate p T [29] and a NRQCD calculation at high p T [30]. On the other hand, in the Collins-Soper frame, the longitudinal component dominates with increasing p T . Thus, λ ϑ becomes negative as p T increases. This relative behavior of λ ϑ in the two frames is expected because the polarization z axes are parallel at p T = 0 and become orthogonal in the limit p T → ∞.
The azimuthal anisotropy parameter (λ ϕ ) reflects the azimuthal symmetry of J/ψ production. When λ ϕ = 0, the production is azimuthally symmetric. When λ ϕ = ±1, the azimuthal distribution is maximally asymmetric. We note that this parameter strongly depends on the production mechanism as well as the frame the distribution is measured in. In the Collins-Soper frame, this parameter is close to zero over all p T as the matrix element σ +1,−1 is small relative to σ +1,+1 and σ 0,0 . This means that the z CS axis is approximately the azimuthal symmetry axis. In the helicity frame, the matrix element σ +1,−1 is negative and become more negative. As a result, λ ϕ becomes negative as p T grows larger, showing that the z HX axis is not the symmetry axis of the distribution. However, the distribution itself is rotationally invariant. The discrepancy between λ ϕ in these two frames is a combination of two factors: z CS and z HX becomes approximately orthogonal as p T increases and production is not spherically symmetric.
The polar-azimuthal correlation parameter (λ ϑϕ ) describes the angular correlation between 2ϑ and ϕ. When λ ϑϕ = 0, the two angles are uncorrelated and as λ ϑϕ departs from 0, the behavior of the distribution becomes similar at locations where 2ϑ = ϕ. In both the helicity frame and the Collins-Soper frame, λ ϑϕ is consistent with 0 at forward rapidity, which agrees with the data. At central rapidity, λ ϑϕ pulls away from 0, but in general within the uncertainties of the data.
We observe a substantial rapidity dependence in the polarization parameters. To illustrate the rapidity dependence, we present the polarization parameters as a function of p T in Fig. 8 in the most central rapidity and the forward rapidity regions covered by the CMS detector and the LHCb detector respectively. In both the Collins-Soper frame and in the helicity frame, the polar anisotropy parameter, λ ϑ switches sign in the high p T limit as rapidity increases. We observe a difference in the curvature of the polarized cross sections as a function of p T in the forward rapidity region compared to the central rapidity region. At central rapidity, λ ϑ becomes constant for p T 20 GeV as the longitudinal and transverse cross sections fall off at similar rates while keeping  [26] in |y| < 0.6 (blue) and 0.6 < |y| < 1.2 (red).  [26] in |y| < 0.6 (blue) and 0.6 < |y| < 1.2 (red).
an approximately fixed ratio. At forward rapidity, in the Collins-Soper frame, the longitudinal cross section first falls off faster than the transverse cross section at moderate p T , but at high p T , the decrease of the longitudinal cross section slows down until it eventually dominates production. We observe the opposite behavior in the helicity frame, where the transverse cross section falls off more slowly than the longitudinal cross section and ultimately dominates production.
We also observe a rapidity dependence in λ ϕ . At central rapidity, the polarized cross section matrix element, σ +1,−1 , from the gg contribution in the helicity frame is negative and relatively constant for p T 20 GeV. As rapidity increases, the matrix element becomes smaller at moderate p T and finally becomes small and positive at high p T . Thus, λ ϕ follows the gg contribution, and the production becomes more azimuthally symmetric with increasing rapidity. In the Collins-Soper frame, the azimuthal dependence of the distribution shows the opposite behavior as rapidity increases. Compared to at central rapidity, the magnitude of λ ϕ at forward rapidity is larger at both moderate p T and high p T . As a result, the production becomes less azimuthally symmetric as rapidity increases. λ ϑϕ is also rapidity dependent because the magnitude of the matrix element, σ 0,+1 in both frames is smaller at forward rapidity than at central rapidity. Note that when the polarization parameters are close to zero, such as λ ϕ and λ ϑϕ , the variation with rapidity is harder to visualize because a given fractional change in the size of a matrix element corresponds to a smaller change in the polarization parameter.
The resulting angular distributions in the helicity frame and in the Collins-Soper frame at p T = 12 GeV are shown in Figs. 9(a) and 9(b), respectively. Note that two distributions are almost identical except that one appears to be rotated 90 • with respect to the other. Thus, the two distributions can be interpreted as a top view and a side view of the production distributions. The angular distribution in the Collins-Soper frame, based on the data collected in the 10 < p T < 15 GeV bin by the LHCb Collaboration [24], is presented in Fig. 9(c).
We compare the polar anisotropy parameter, λ ϑ , in this calculation with the one calculated in the polarized ICEM using the k T -factorization approach [13] in Fig. 10. Since the k T -factorized ICEM considers only the contribution from off-shell Reggeized gluons at O(α 2 s ), we compare the p T dependence of λ ϑ in the k T -factorized ICEM with that in this calculation using the contribution from gg → ccg at O(α 3 s ) only. The polarization calculated in the gg channel alone is approximately equivalent to that summed over all channels in the collinear factorization approach because the gg component dominates both GeV. They represent the same angular distribution separated by one rotation. We constructed the angular distribution in the Collins-Soper frame based on the data collected in the 10 < pT < 15 GeV bin by the LHCb Collaboration [24], and it is shown in (c) for comparison.
the polarized and unpolarized cross sections. Although the kinematics in the two calculations are different, they both predict longitudinal polarization at high p T , where the experimental results are nearly unpolarized. At small and moderate p T , the difference in λ ϑ is within the uncertainty of the data.
C. pT and y dependence ofλ In Figs. 6 and 7, we observe that, at high p T , the ICEM is in better agreement with the measured data in the helicity frame than in the Collins-Soper frame. However, even though we are switching from one frame to another, we are still comparing the same angular distributions. The difference between the ICEM polarization results and the data is then best quantified by a frameindependent polarization parameter. Thus, we compute the frame-invariant polarization parameterλ as a function of p T using λ ϑ and λ ϕ . We compareλ as a function of p T with the data from LHCb and CMS in the helicity and the Collins-Soper frames in Fig. 11.
Since the azimuthal anisotropy parameter λ ϕ in the Collins-Soper frame is close to zero in all p T and y ranges considered, the p T and y dependences of the invariant polarization parameterλ are very similar to those in λ ϑ of the Collins-Soper frame. We find the curves are generally within 1σ of the low and moderate p T data, as shown in Fig. 11 (a). We also find the invariant polarization of direct J/ψ is relatively constant at high p T . We observe a rapidity dependence similar to that of the λ ϑ in the high p T limit, where the invariant polarization is less longitudinal in the forward rapidity region than in the central rapidity region. We note that the invariant polarization parameter was devised as a way to treat each frame equally rather than to be the final arbiter between √ s = 7 TeV with 2 < y < 4.5 in the Collins-Soper frame calculated in the collinear ICEM in this work (blue) and in the kT -factorized ICEM in Ref. [13]. Both calculations include only the gg contributions and are compared to the LHCb data [24]. theoretical model calculations. Our direct J/ψ-invariant polarization results are in reasonable agreement with the measured data over all p T even though our calculations show a decrease inλ with p T . We note that including feed down from higher mass states could change the slope in this p T range. Our results are in agreement with the J/ψ-invariant polarization found in the CGC+NRQCD approach at low and moderate p T [29].

V. CONCLUSIONS
We have presented the transverse-momentum dependence of the direct J/ψ cross section as well as the polarization in p + p collisions in the improved color evaporation model in the collinear factorization approach. We compare the p T and y dependences to data for inclusive J/ψ production measured by the LHCb, ALICE, and CMS Collaboration. We also present the frame-invariant parameter,λ, as a function of p T and compare our results with the data. We find direct J/ψ production is consistent with the unpolarized data at small and moderate p T and becomes longitudinal in the high p T limit. In the near future we will see if the feed-down contribution will have a positive influence on the parameter to match the data in that region. We will study the effects of feed-down production in this approach in a future publication.