Are there any challenges in the charmonia production and polarization at the LHC?

We analyze the inclusive prompt production of $\psi^\prime$, $\chi_c$, $J/\psi$ and $\eta_c$ mesons at the LHC using the $k_T$-factorization approach. Our consideration is based on the off-shell production amplitudes for hard partonic subprocesses, nonrelativistic QCD formalism for the formation of bound states and transverse momentum dependent (or unintegrated) gluon densities in a proton derived from Catani-Ciafaloni-Fiorani-Marchesini equation. The nonperturbative color octet transitions are treated in terms of the multipole radiation theory. We extract the corresponding long-distance matrix elements from a combined fit to transverse momentum distributions measured at various LHC experiments. We make predictions for the polarization parameters $\lambda_\theta$, $\lambda_\phi$, $\lambda_{\theta \phi}$ and the frame-independent parameter $\lambda^*$ and compare them to the available $\psi^\prime$ and $J/\psi$ data. Finally, we present a universal set of parameters that provides a reasonable simultaneous description for the whole body of the LHC data (on the $p_T$ distributions, relative production rates and polarization observables) for the whole charmonium family.


Introduction
Since it was first observed, the prompt charmonia production in hadronic collisions remains a topic of considerable theoretical and experimental interest. A commonly accepted theoretical framework for the description of charmonia production and decay is provided by the nonrelativistic QCD (NRQCD) factorization [1,2]. In this formalism, the perturbatively calculated cross sections for the short-distance production of a cc pair in an intermediate Fock state 2S+1 L (a) J with spin S, orbital angular momentum L, total angular momentum J, and color representation a are accompanied with long-distance matrix elements (LDMEs), which describe the transition from that intermediate cc state into a physical meson via soft gluon radiation. The LDMEs are assumed to be universal (process-and energy-independent) and obeying certain hierarchy in powers of the relative charmed quarks velocity v.
At present, the cross sections of ψ , χ c , J/ψ and η c production in pp collisions are known at the next-to-leading order (NLO) accuracy [3][4][5][6][7][8][9][10][11][12][13][14][15]. The dominant tree-level next-to-next-toleading order (NNLO * ) corrections to the color-singlet (CS) mechanism have been calculated [16]. The long-distance matrix elements (LDMEs) are not calculable within the theory and can only be extracted from fits to the data. Then, with properly adjusted LDMEs values, a reasonably good description of the ψ , χ c and J/ψ transverse momentum distributions measured at the Tevatron and LHC is achieved [3][4][5][6][7][8][9][10]. However, the extracted LDMEs dramatically depend on the minimal charmonium transverse momentum used in the fits and are incompatible with one another when obtained from fitting different data sets. Moreover, none of the fits is able to reasonably accommodate the polarization measurements (the so called "polarization puzzle").
The fits involving low p T data lead to the conclusion that ψ and J/ψ production at large transverse momenta must be dominated by the 3 S [8] 1 contributions with strong transverse polarization, that contradicts to the unpolarized production seen by the CDF Collaboration at the Tevatron [17,18] and CMS [19] and LHCb [20,21] Collaborations at the LHC. To obtain an unpolarized J/ψ meson, it is necessary to assume that its production is dominated by the scalar 1 S [8] 0 intermediate state [4]. This comes to an immediate conflict with the recent LHCb data [22] on the η c meson production, as the respective η c and J/ψ LDMEs are related by one of the basic NRQCD principles, the heavy quark spin symmetry (HQSS). The impact of the η c data [22] on the general understanding of the charmonia production and polarization phenomena was investigated in [12]. The overall situation is found difficult and has been even called 'challenging' [13]. At present, the conventional NRQCD is yet far from understanding the data (see also discussions [23][24][25]). So, the further theoretical studies are still an urgent task.
Recently, a solution to the polarization puzzle has been proposed [26] in the framework of a model that interprets the soft final state gluon radiation as a series of color-electric dipole transitions. In this way the LDMEs are represented in a form ispired by the classical multipole radiation theory, so that the spin structure of the transition amplitudes is explicitly specified. The calculations made in this approach lead to weak final J/ψ polarization, either because of the cancellation between the 3 P [8] 1 and 3 P [8] 2 contributions, or as a result of two successive color-electric (E1) dipole transitions in the chain 3 S [8] 1 diluting 1 S [8] 0 contribution to J/ψ, we neither need its HQSS counterpart process, the 3 S contribution to η c , while the production of η c is saturated by the color singlet mechanism alone.
We follow this approach [26] in the present paper and carry out a global study of the production and polarization phenomena for the entire charmonium family (ψ , χ c , J/ψ and η c ) at the LHC. To describe the perturbative production of a cc pair in hard scattering subprocess we employ the k T -factorization formalism [31,32]. This formalism is based on the Balitsky-Fadin-Kuraev-Lipatov (BFKL) [33] or Catani-Ciafaloni-Fiorani-Marchesini (CCFM) [34] evolution equations and has certain technical advantages in the ease of including higher-order radiative corrections (namely, a part of NLO + NNLO +... terms corresponding to the initial-state real gluon emissions) in the form of transverse momentum dependent (TMD) gluon densities 1 . Then we perform a simultaneous fit for charmonia LDMEs using the latest LHC data collected by the ATLAS [37][38][39], CMS [40][41][42] and LHCb [43][44][45][46][47][48][49] Collaborations at √ s = 7, 8, and 13 TeV 2 . We also pay attention to the relative production rates, for example, σ(χ c2 )/σ(χ c1 ), as the latter are sensitive to the color siglet and color octet production mechanisms. A clear understanding of χ c (and, of course, ψ ) production is an important component of any general description of J/ψ production and polarization since the feed-down contributions from radiative decays constitute about 30% of the visible J/ψ cross section at the LHC. Using the fitted LDMEs, we make predictions for the polarization parameters λ θ , λ φ and λ θφ and compare them to the currently available data on ψ and J/ψ mesons. Our main goal is to show that the consistently used approach [26] meets no troubles with the available charmonia data (including transverse momentum distributions, relative production rates and polarization observables). We end up with presenting a universal set of parameters that provides a reasonable simultaneous description of everything.
The outline of our paper is the following. In Section 2 we describe the basic steps of our calculations. In Section 3 we perform a numerical fit and extract the charmonia LDMEs from the latest LHC data [37][38][39][40][41][42][43][44][45][46][47][48][49] on the transverse momentum distributions. Later in this section we check the compatibility of the extracted LDMEs with the available data [50] on charmonia polarization. The comparison is followed by a discussion. Our final conclusions are collected in Section 4.

Theoretical framework
We start with recalling the essential calculation steps. Our consideration is based on the following leading-order off-shell gluon-gluon fusion subprocesses for ψ and J/ψ mesons: for χ c mesons (with J = 0, 1, 2): 1 See reviews [35,36] for more information. 2 In our previous studies [27][28][29] such fits were performed using the LHC data at √ s = 7 TeV only. and for η c mesons: Additionally, we took into account the feed-down contribution to η c production from the decays h c → η c γ. The leading contributions to h c come from the off-shell partonic subprocesses where the four-momenta of all particles are indicated in the parentheses. The corresponding production amplitudes contain spin projection operators which discriminate the spin-singlet and spin-triplet cc states [51]: where m c is the charmed quark mass. States with various projections of the spin momentum onto the z axis are represented by the polarization four-vector (S z ). Here p c and pc are the four-momenta of the charmed quark and anti-quark: The relative momentum q of the quarks in a bound state is associated with the orbital angular momentum L. According to the general formalism [52,53], the terms showing no dependence on q are identified with the contributions to the L = 0 states while the terms linear (quadratic) in q µ are related to the L = 1 (L = 2) states with the proper polarization vector µ (L z ) (resp., polarization tensor µν (L z )). The hard scattering amplitude A(q) has to be multiplied by the bound state wave finction Ψ [a] (q) and integrated over q. The integration is done after expanding the amplitude around q = 0: A term-by-term integration of this series employs the identities [51] where R [a] (x) is the radial wave function in the coordinate representation. The first term in (10) only contributes to S-waves, but vanishes for P -waves. On the contrary, the second term only contributes only to P -waves, but vanishes for S-waves. The corresponding LDMEs are related to the wave functions R [a] (x) and their derivatives [1,2] as for S-waves and for P -waves. All algebraic calculations are straightforward and have been done in our previous papers [27][28][29][30]. The resulting expressions have been tested for gauge invariance by substituting the gluon momenta for corresponding polarization vectors. We have observed gauge invariance even with off-shell initial gluons. Now, let us turn to non-perturbative ingredients of the theory. As it is motivated by the HQSS relations, the LDMEs should be identical for transitions in both directions (from vectors to scalars and vice versa) and can only differ by an overall normalizing factor representing the averaging over spin degrees of freedom. Thus, we strictly have from this property [1,2]: O ηc 1 P O hc 1 P O hc 1 S for all S-and P -wave bound states Q and color states a. The relations between the different LDMEs require that the fit be done simultaneously for the entire charmonium family. Following the ideas of [26], we employ the classical multipole radiation theory to describe nonperturbative transformations of the color-octet cc pairs produced in hard subprocesses into observed final state charmonia. Only a single E1 transition is needed to transform a P -wave state into an S-wave state, and the structure of the respective 3 P 1 + g amplitudes is taken as [54] A( 3 P where k and l = p − k are the four-momenta of the emitted gluon and the produced meson, µ (k), µ (l), µ (p) and µν (p) are the polarization vectors (tensor) of the respective particles and e µναβ is the fully antisymmetric Levi-Civita tensor. The transformation of an S-wave state into another S-wave state (such as J/ψ meson) is treated as two successive E1 transi- [8] , or 3 P 2 [8] . For each of the two transitions we exploit the same effective coupling vertices (23) - (25). Note that the expressions describing E1 transitions are the same for gluons and photons (up to an overall color factor) and therefore can also be used to calculate the polarization variables in radiative decays in feed-down processes ψ → χ cJ + γ and χ cJ → J/ψ + γ. The polarization of the outgoing mesons can then be calculated without any ambiguity. The production cross section for a charmonium Q is calculated as a convolution of the off-shell partonic cross section and the TMD gluon densities in a proton. We have for the 2 → 1 and 2 → 2 subprocesses, respectively: where f g (x, k 2 T ), µ 2 ) is the transverse momentum dependent (TMD, or unintegrated) gluon density in a proton, p T and y are the transverse momentum and rapidity of produced charmonium Q, y g is the rapidity of outgoing gluon and √ s is the pp center-of-mass energy. The initial off-shell gluons have fractions x 1 and x 2 of the parent protons longitudinal momenta, non-zero transverse momenta k 1T and k 2T and azimuthal angles φ 1 and φ 2 . In accordance with the general definition [55], the off-shell gluon flux factor in (26) is taken as F = 2λ 1/2 (ŝ, k 2 1 , k 2 2 ), whereŝ = (k 1 + k 2 ) 2 . In the numerical analysis below, we have tried a few sets of TMD gluon densities in a proton, referred to as A0 [56], JH'2013 set 1 and JH'2013 set 2 [57]. These gluon densities were obtained from CCFM evolution equation where the input parametrization (used as boundary conditions) was fitted to the proton structure function F 2 (x, Q 2 ) and, in the case of JH'2013 set 2, to F c 2 (x, Q 2 ) also. The CCFM equation provides a suitable tool for our phenomenologycal study since it smoothly interpolates between the small-x BFKL gluon dynamics and high-x DGLAP one. The renormalization and factorization scales were set to where m Q and Q T are the produced charmonium Q mass and the transverse momentum of the initial off-shell gluon pair, respectively. The choice of µ R is rather standard for charmonia production, while the unusual choice of µ F is connected with the CCFM evolution (see [56,57] for more details). The multidimensional phase space integration has been performed by means of the Monte-Carlo technique using the routine vegas [58].

Global fit of charmonia LDMEs based on the LHC data
We have performed a global fit to the charmonium production data at the LHC for the entire cc family and determined the corresponding LDMEs. Specifically, for ψ mesons, we included in the fitting procedure the transverse momentum distributions measured by ATLAS [38,39] and CMS [41,42] Collaborations at moderate and large transverse momenta 8 < p T < 130 GeV at √ s = 7, 8 and 13 TeV, where the NRQCD formalism is believed to be most reliable. We have excluded from our fit the LHCb data [45] since they mainly lie in the low p T region, where a more accurate treatment of large logarithms ln m 2 ψ /p 2 T and other nonperturbative effects becomes necessary 4 . In the case of χ c mesons, we considered the χ c1 and χ c2 transverse momentum distributions measured by ATLAS Collaboration [37] at √ s = 7 TeV and also include in the fitting procedure the ratio of the production rates σ(χ c2 )/σ(χ c1 ) measured by CMS [40], ATLAS [37] and LHCb [43,44] Collaborations at the same energy. Note that most of the theoretical uncertainties cancel out in the ratio; in particular, the uncertainties related to the behavior of the TMD gluon densities in the low-p T region.
Following the suggestion [61], we consider the CS wave functions of χ c1 and χ c2 mesons as independent (not necessarily identical) parameters. The reasoning for such a suggestion is that treating charmed quarks as spinless particles (as in the potential models [62][63][64][65]) might be an oversimplification, and that radiative corrections to the wave functions may be large 5 .
To determine the LDMEs of J/ψ mesons (as well as their η c counterparts) we performed a simultaneous fit of J/ψ and η c transverse momentum distributions using the latest CMS [41,42], ATLAS [38] and LHCb [22] data taken at 7, 8 and 13 TeV. Here, the NRQCD factorization principle seems to be on solid theoretical grounds again because of not too low p T values for both J/ψ and η c (at least, p T > 8 GeV for η c mesons). Of course, we took into account the feed-down contributions to J/ψ and η c production from radiative decays of χ c , ψ and h c mesons using the corresponding branching fractions as listed above.
Nowhere we impose any kinematic restrictions but the experimental acceptance. The fitting procedure was separately done in each of the rapidity subdivisions (using the fitting algorithm as implemented in the gnuplot package [66]) under the requirement that the LDMEs be strictly positive, and then the mean-square average of the fitted values was taken. The relevant uncertainties are estimated in the conventional way using Student's t-distribution at the confidence level P = 95%.
We observe in Figs. 1 -9 quite a nice agreement between our calculations and the LHC data for the entire charmonium family at different energies and in a wide p T range for all of the considered TMD gluons (with the LDMEs values shown in Table 1). In particular, we have achieved good simultaneous description of the prompt η c and J/ψ production, see Figs. 6 -9. Such an agreement turned out to be impossible in the traditional NRQCD scheme, where the calcullated cross sections for η c are either at odds with the measurements or at odds with theoretical principles [11,13]. Further on, we have achieved a good agreement with the LHCb data [45][46][47][48][49], originally not included in the fitting procedure (see Figs. 3 and 9). The extracted LDMEs values strongly depend on the TMD gluon density (see Table 1), that reflects the different x and k T behavior of the latter (see discussion [57]). The estimated theoretical uncertainties of our predictions are rather small and comparable with the uncertainties of the NLO NRQCD calculations.
Our fits show unequal values for the χ c1 and χ c2 wave functions at the origin |R [1] (0)| 2 . We present these values in Table 1. The difference in the values of the wave functions mainly follows from the prompt measurements of the ratio σ(χ c1 )/σ(χ c2 ). For each of the considered gluon densities, our extracted values of |R [1]χ c2 (0)| 2 (but not |R [1]χ c1 (0)| 2 ) are close to the estimations based on the potential models [62][63][64][65] and two-photon decay width [59]; namely, |R [1] (0)| 2 = 7.5·10 −2 GeV 5 (that is a widely adopted choice). However, it differs significantly from |R [1] (0)| 2 = 3.5 · 10 −1 GeV 5 obtained from a combined fit [9] to the Tevatron and LHC data. Note that the fit [9] was performed under the assumption of equal χ c1 and χ c2 wave functions. We interpret the available LHC data [37,40,43,44] as supporting their unequal values, that qualitatively agrees with the previous results [28,61]. In such an interpretation, the data on the σ(χ c2 )/σ(χ c1 ) ratio lie almost inside the theoretical uncertainty bands, as one can see in Fig. 5. Moreover, the ratio |R [1]χ c2 (0)| 2 /|R [1]χ c1 (0)| 2 2.5 is practically independent on the TMD gluon density. Finally, we find that χ c production is dominated by the CS contributions in the considered p T range, that agrees with some earlier conclusions [9]. The obtained LDMEs for ψ and χ c mesons were further used to calculate the feeddown contributions to J/ψ production. The results of our fits for J/ψ and ψ polarization parameters are discussed in the next Section.

J/ψ and ψ polarization
It is known that the polarization of ψ or J/ψ mesons can be described with three parameters λ θ , λ φ and λ θφ , which determine the spin density matrix of a charmonium decaying into a lepton pair. In general, the double differential angular distribution of the decay leptons in the charmonium rest frame can be written as dσ d cos θ * dφ * ∼ 1 3 + λ θ 1 + λ θ cos 2 θ * + λ φ sin 2 θ * cos 2φ * + λ θφ sin 2θ * cos φ * , where θ * and φ * are the polar and azimuthal angles of the decay lepton. So, the angular parameters λ θ , λ φ and λ θφ can be measured experimentally. The case of (λ θ , λ φ , λ θφ ) = (0, 0, 0) corresponds to unpolarized state, while (λ θ , λ φ , λ θφ ) = (1, 0, 0) and (λ θ , λ φ , λ θφ ) = (−1, 0, 0) refer to fully transverse and longitudinal polarizations. The CMS Collaboration has measured [50] all these parameters as functions of J/ψ and ψ transverse momentum in the complementary frames: the Collins-Soper, helicity and perpendicular helicity ones 6 . In the Collins-Soper frame the polarization axis z bisects the two beam directions whereas the polarization axis in the helicity frame coincides with the direction of the charmonium momentum in the laboratory frame. In the perpendicular helicity frame the z axis is orthogonal to that in the Collins-Soper frame and lies in the plane spanned by the two beam (P 1 and P 2 ) momenta. In all cases, the y axis is taken to be in the direction of the vector product of the two beam directions in the charmonium rest frame, P 1 × P 2 and P 2 × P 1 for positive and negative rapidities, respectively. Additionally, the frame-independent polarization parameter λ * = (λ θ + 3λ φ )/(1 − λ φ ) was investigated [50].
To estimate the polarization parameters λ θ , λ φ , λ θφ and λ * we generally follow the experimental procedure. We collect the simulated events in the kinematical region defined by the CMS measurement [50], generate the decay lepton angular distributions according to the production and decay matrix elements and then apply a three-parametric fit based on (27). Of course, in the case of J/ψ production we took into account the polarization of J/ψ mesons originated from radiative χ c and ψ decays, that is in full agreement with the experimental case. Since the ψ → J/ψ+X decay matrix elements are unknown, these events were generated according to the phase space. In Figs. 10 -12 we confront our predictions for all polarization parameters with the CMS data [50]. For both J/ψ and ψ mesons we find only weak polarization (λ θ −0.2) at p T ∼ 15 GeV in the Collins-Soper and helicity frames and practically zero polarization (λ θ −0.1 or even close to zero) at large transverse momenta p T ∼ 50 GeV. In the perpendicular helicity frame our simulation shows practically unpolarized J/ψ and ψ production with λ θ ∼ 0 in the whole p T range. The λ φ and λ θφ parameters are close to zero everywhere, as one can see in Figs. 10 -12. Moreover, these results are practically independent of the J/ψ and/or ψ rapidity. Thus, we demonstrate that treating the soft gluon emission within the NRQCD as a series of explicit color-electric dipole transitions leads to unpolarized charmonia production, that is in agreement with available LHC data. The absense of strong polarization is not connected with parameter tuning, but seems to be a natural and rather general feature of the scenario [26]. We would like to point out here that the conventional NLO CS calculations predict large longitudinal charmonia polarization at high transverse momenta, while the NLO NRQCD predicts large transverse polarization. None of these predictions is supported by the LHC measurements.
The obtained unpolarized J/ψ and ψ production at the LHC is our main result. The qualitative predictions for the λ θ , λ φ , λ θφ and λ * are stable with respect to variations in the model parameters. In fact, there is no dependence on the strong coupling constant and TMD gluon densities, i.e. two of important sources of theoretical uncertainties cancel out. Despite large experimental uncertainties (especially for λ * parameter), the agreement between our predictions and the data is rather satisfactory and shows no fundamental problems in describing the data. So, the proposed way, in our opinion, can provide an easy and natural solution to the long-standing polarization puzzle.

Conclusion
We have considered the inclusive prompt production of ψ , χ c , J/ψ and η c mesons at the LHC in the framework of k T -factorization approach. Our consideration was based on the offshell production amplitudes for hard partonic subprocesses (including both color-singlet and color-octet contributions) and NRQCD formalism for the formation of bound states. Treating the nonperturbative color octet transitions in terms of multipole radiation theory and applying the TMD gluon densities in a proton derived from the CCFM evolution equation, we extracted charmonia LDMEs in a combined fit to transverse momentum distributions measured on various LHC experiments at √ s = 7, 8 and 13 TeV. Then, using the extracted LDMEs, we estimated polarization parameters λ θ , λ φ , λ θφ and frame-independent parameter λ * which determine the charmonia spin density matrix. We have demonstrated that treating the soft gluon emission as a series of explicit color-electric dipole transitions within the NRQCD leads to unpolarized charmonia production at moderate and large transverse momenta, that is in agreement with the recent LHC data on ψ and J/ψ mesons. Thus, we achieved a reasonable simultaneous description for all of the available data (transverse momentum distributions, relative production rates and polarization observables) on the entire charmonia family at the LHC.   Fig. 1. The experimental data are from CMS [41,42] and LHCb [45].    Figure 5: The relative production rate σ(χ c2 )/σ(χ c1 ) calculated as a function of decay J/ψ meson transverse momenta at √ s = 7 TeV. Notation of all curves is the same as in Fig. 1. The experimental data are from ATLAS [37], CMS [40] and LHCb [43,44].       |y| < 0.6 0.6 < |y| < 1.2 1.2 < |y| < 1.5 Figure 10: Polarization parameters λ θ , λ φ , λ θφ and λ * of prompt J/ψ (left panels) and ψ (right panels) mesons calculated as a function of transverse momentum in the Collins-Soper frame. The yellow, blue and green histograms correspond to the predictions obtained at |y| < 0.6, 0.6 < |y| < 1.2 and 1.2 < |y| < 1.5. The JH'2013 set 2 gluon density is used. The experimental data are from CMS [50]. |y| < 0.6 0.6 < |y| < 1.2 1.2 < |y| < 1.5 Figure 11: Polarization parameters λ θ , λ φ , λ θφ and λ * of prompt J/ψ (left panels) and ψ (right panels) mesons calculated as a function of transverse momentum in the helicity frame. Notation of all histograms is the same as in Fig. 8. The experimental data are from CMS [50]. |y| < 0.6 0.6 < |y| < 1.2 1.2 < |y| < 1.5 Figure 12: Polarization parameters λ θ , λ φ , λ θφ and λ * of prompt J/ψ (left panels) and ψ (right panels) mesons calculated as a function of transverse momentum in the perpendicular helicity frame. Notation of all histograms is the same as in Fig. 8. The experimental data are from CMS [50].