1 Introduction

The production of quarkonium states in high energy hadronic collisions has been under intense theoretical and experimental study [13] since two decades ago, when the measurements of prompt \(J/\psi \) and \(\Upsilon \) production cross sections at the Tevatron revealed a more than one order-of-magnitude discrepancy with theoretical expectations obtained in the framework of the color singlet model [47]. This fact has induced extensive theoretical activity mainly connected with modeling the formation of quarkonium states from unbound heavy quark pairs produced in hard interactions. There exist two competing theoretical approaches known in the literature under the names of color-singlet (CS) and color-octet (CO) models. In general, a quark–antiquark pair is produced in a state \(^{2S+1}L_J^{(a)}\) with spin S, orbital angular momentum L, total angular momentum J, and color a, which can be either identical to the respective quantum numbers of the resulting quarkonium (as accepted in the CS model) or different from those. In the latter case, the heavy quark pair transforms into a physical quarkonium state by means of soft (non-perturbative) gluon radiation, as considered in the formalism of non-relativistic quantum chromodynamics (NRQCD) [810]. The quarkonium formation probability is then determined by the respective non-perturbative matrix elements (NMEs), which are assumed to be universal (process independent) and not depending on the quarkonium momentum. Though not strictly calculable within the theory, the NMEs are assumed to obey a certain hierarchy in powers of the relative quark velocity v. To the leading order in v, an S-wave vector meson such as \(\psi (2S)\) can be formed from a quark pair produced as a color singlet \(^3S_1^{(1)}\) or via one of the intermediate color octet states \(^1S_0^{(8)}\), \(^3S_1^{(8)}\) or \(^3P_J^{(8)}\) with \(J = 0\), 1 or 2.

We know already that the CS model with leading-order (LO) hard scattering matrix elements fails to describe the experimental data on \(J/\psi \) and \(\Upsilon \) production at the Tevatron and LHC. Including the next-to-leading order (NLO) [11] and dominant tree-level next-to-next-to-leading order (NNLO\(^*\)) [12] corrections to the CS mechanism significantly improves the description of the collider data [13]. In the NRQCD formalism, reasonably good agreement with the measured quarkonia cross sections can be achieved by adjusting the NME values, which play the role of free fitting parameters [1419]. This was already demonstrated by comparing the calculations with the ATLAS [20], CMS [21] and LHCb [22] experimental data taken at \(\sqrt{s} = 7\) TeV. However, the values of the extracted NMEs dramatically depend on the minimal quarkonium transverse momentum \(p_T\) used in the fits [23] and are incompatible with each other (both in size and even in sign!) when obtained from fitting the different sets of data. Furthermore, none of the fits is able to accommodate the polarization measurements [24, 25]. The fits involving low \(p_T\) measurements lead to the conclusion that the production of S-wave quarkonia at high \(p_T\) must be dominated by CO contributions with transverse polarization (namely, by the \(^3S_1^{(8)}\) channel). The latter contradicts the unpolarized production seen by the CDF Collaboration [26, 27] at the TevatronFootnote 1 and CMS [24, 28] and LHCb [25] Collaborations at the LHC. To obtain an unpolarized state it is necessary to either assume that quarkonium production is dominated by the scalar \(^1S_0^{(8)}\) channel [16] or restrict the NRQCD fits to very high \(p_T\) region [23]. This problem is known as the “quarkonium polarization puzzle” and is still far from understanding.

A new solution to the polarization puzzle has been guessed in the \(k_T\)-factorization approach of QCD [2932], where studies of quarkonia production and polarization have their own long history (see, for example, [3340] and references therein). The \(k_T\)-factorization approach is based on the Balitsky–Fadin–Kuraev–Lipatov (BFKL) [4143] or Ciafaloni–Catani–Fiorani–Marchesini (CCFM) [4447] gluon evolution equations and provides theoretical grounds for including the effects of initial gluon radiation and intrinsic gluon transverse momentum.Footnote 2 A combination of the usual CS scheme and the \(k_T\)-factorization formalism results in a reasonably good description of the data on \(J/\psi \), \(\chi _c\), and \(\Upsilon \) production at HERA [33, 35, 36, 39], Tevatron [34, 38], and LHC [40]. These results are also in good agreement with more complicated explicit NNLO\(^*\) calculations [1113] performed in the collinear approximation of QCD. The longitudinal polarization of directly produced quarkonia in the \(k_T\)-factorization is an immediate consequence of initial gluon off-shellness [33]. Adding the feed-down from P-wave states (\(\chi _c\) and \(\chi _b\)) leads to essentially unpolarized prompt \(J/\psi \) and \(\Upsilon \) mesons [38, 40]. In the \(k_T\)-factorization approach at LO, the P-wave states are produced in \(2\rightarrow 1\) gluon–gluon fusion and are expected to dominate at high transverse momenta. However, the latest LHC data show that the feed-down contributions \(\chi _c\rightarrow J/\psi {+}\gamma \) and \(\chi _b\rightarrow \Upsilon {+}\gamma \) do not constitute more than only 20–30 % of the visible cross section at large \(p_T\) values.

So, the production of \(\chi _c\) and \(\chi _b\) mesons requires a dedicated study which will be the subject of our forthcoming papers, while in the present analysis we concentrate on the direct mechanism and only restrict our study to \(\psi (2S)\) mesons having no contamination from higher states. Here we present a systematic analysis of ATLAS [20], CMS [21], and LHCb [22] data collected at \(\sqrt{s} = 7\) TeV regarding the transverse momentum distributions and polarization parameters \(\lambda _\theta \), \(\lambda _\phi \), and \(\lambda _{\theta \phi }\), which describe the spin density matrix of the produced \(\psi (2S)\) mesons.

The outline of our paper is as follows. In Sect. 2 we briefly recall the NRQCD formalism and the \(k_T\)-factorization approach. In Sect. 3 we perform a numerical fit to the latest CMS data and extract the color-octet NMEs using three different sets of transverse momentum dependent (TMD) gluon distributions. Later in this section we check the compatibility of the extracted parameters with ATLAS and LHCb data on the production and polarization of \(\psi (2S)\) mesons. The comparison is followed by a discussion. Our conclusions are collected in Sect. 4.

2 Theoretical framework

We start by briefly recalling the essential calculation steps. Our consideration is based on the following leading-order off-shell partonic subprocesses [810]:

$$\begin{aligned} g^*(k_1) + g^*(k_2)\rightarrow & {} \psi ^\prime \left[ ^3S_1^{(1)}\right] (p) + g(k), \end{aligned}$$
(1)
$$\begin{aligned} g^*(k_1) + g^*(k_2)\rightarrow & {} \psi ^\prime \left[ ^1S_0^{(8)}, \, ^3S_1^{(8)}, \, ^3P_J^{(8)}\right] (p), \end{aligned}$$
(2)

where \(J = 0, 1\) or 2, and the four-momenta of all particles are indicated in parentheses. The subprocesses (1) and (2) represent the CS and CO contributions, respectively. The corresponding production amplitudes can be obtained from the one for an unspecified \(c\bar{c}\) state by applying the appropriate projection operators, which guarantee the proper quantum numbers of the \(c\bar{c}\) state under consideration. These operators for the different spin and orbital angular momentum states read [47]

$$\begin{aligned} \Pi \left[ ^1S_0\right]= & {} \gamma _5 \left( \hat{p}_c + m_c\right) /m^{1/2},\end{aligned}$$
(3)
$$\begin{aligned} \Pi \left[ ^3S_1\right]= & {} \hat{\epsilon }(S_z) \left( \hat{p}_c + m_c\right) /m^{1/2},\end{aligned}$$
(4)
$$\begin{aligned} \Pi \left[ ^3P_J\right]= & {} \left( \hat{p}_{\bar{c}} - m_c\right) \hat{\epsilon }(S_z) \left( \hat{p}_c + m_c\right) /m^{3/2}, \end{aligned}$$
(5)

where \(m=2m_c\) is the mass of the considered \(c\bar{c}\) state, and \(p_c\) and \(p_{\bar{c}}\) are the four-momenta of the charmed quark and antiquark. States with various projections of the spin momentum onto the z axis are represented by the polarization four-vector \(\epsilon _\mu (S_z)\).

The probability for the two quarks to form a meson depends on the bound state CS and fictitious CO wave functions \(\Psi ^{(a)}(q)\), where the relative four-momentum q of the quarks in the bound state is treated as a small quantity in the non-relativistic approximation. So, we represent the quark momenta as

$$\begin{aligned} p_c = p/2 + q, \, p_{\bar{c}} = p/2 - q. \end{aligned}$$
(6)

Then we multiply the hard subprocess amplitude \({\mathcal A}\) (depending on q) by the meson wave function \(\Psi ^{(a)}(q)\) and perform integration with respect to q. The integration is done after expanding the amplitude \({\mathcal A}\) around \(q = 0\):

$$\begin{aligned}&{\mathcal A}(q)\Psi ^{(a)}(q) = {\mathcal A}|_{q = 0}\,\Psi ^{(a)}(q) \nonumber \\&\quad + q^\alpha (\partial {\mathcal A} / \partial q^\alpha )|_{q = 0}\,\Psi ^{(a)}(q) + \cdots \end{aligned}$$
(7)

Since the expressions for \({\mathcal A}|_{q = 0}\) and \(\partial {\mathcal A} / \partial q^\alpha |_{q = 0}\) are no longer dependent on q, they may be factored outside the integral sign. A term-by-term integration of this series employs the identities [47]:

$$\begin{aligned} \int {\text {d}^3 q\over (2 \pi )^3} \, \Psi ^{(a)}(q)= & {} {1\over \sqrt{4\pi }} {\mathcal R}^{(a)}(0),\end{aligned}$$
(8)
$$\begin{aligned} \int {\text {d}^3 q\over (2 \pi )^3} \, q^\alpha \Psi ^{(a)}(q)= & {} - i \epsilon ^\alpha (L_z){\sqrt{3}\over \sqrt{4\pi }} {\mathcal R}^{\prime \, (a)}(0), \end{aligned}$$
(9)

where \({\mathcal R}^{(a)}(x)\) are the radial wave functions in the coordinate representation. The first term in (7) contributes only to S-waves, but it vanishes for P-waves. On the contrary, the second term contributes only to P-waves, but it vanishes for S-waves. States with various projections of the orbital angular momentum onto the z axis are represented by the polarization four-vector \(\epsilon _\mu (L_z)\).

The NMEs of S-wave states are directly related to the CS and fictitious CO wave functions:

$$\begin{aligned} \left\langle {\mathcal O}^{\psi }\left[ {^{2S + 1}L_J^{(a)}}\right] \right\rangle = 2 N_c (2J + 1) |{\mathcal R}^{(a)}(0)|^2 /4\pi , \end{aligned}$$
(10)

where \(N_c=3\) and \(J=1\). A similar relation holds for \({\mathcal R}^{\prime \,(a)}\) if P-wave states are involved. The CS wave function at the origin of coordinate space is known from the measured \(\psi (2S)\) leptonic decay width. The color-octet NMEs are extracted from experimental data and obey the relation

$$\begin{aligned} \left\langle {\mathcal O}^{\psi }\left[ {^{3}P_J^{(8)}}\right] \right\rangle = (2J + 1) \, \left\langle {\mathcal O}^{\psi }\left[ {^{3}P_0^{(8)}}\right] \right\rangle , \end{aligned}$$
(11)

coming from heavy quark spin symmetry at LO in v. The polarization vectors \(\epsilon _\mu (S_z)\) and \(\epsilon _\mu (L_z)\) are defined as explicit four-vectors. In the frame where the z axis is oriented along the quarkonium momentum vector \(p_\mu = (E,0,0,|{\mathbf p}|)\), these polarization vectors read

$$\begin{aligned} \epsilon _\mu (\pm 1) = (0, \pm 1, i, 0)/\sqrt{2}, \, \epsilon _\mu (0) = (|{\mathbf p}|, 0, 0, E)/m. \end{aligned}$$
(12)

The states with definite \(S_z\) and \(L_z\) are translated into states with definite total momentum J and its projection \(J_z\) using the Clebsch–Gordan coefficients:

$$\begin{aligned} \epsilon _{\mu \nu }(J, J_z) = \sum _{S_z, \, L_z} \langle 1, L_z; 1, S_z | J, J_z \rangle \, \epsilon _\mu (S_z) \, \epsilon _\nu (L_z). \end{aligned}$$
(13)

Further evaluation of partonic amplitudes is straightforward and is done using the algebraic manipulation system form [51]. Our results for perturbative production amplitudes squared and summed over polarization states agree with the ones in [52, 53]. Here we only mention several technical points. First, according to the \(k_T\)-factorization prescription [2932], the summation over the incoming off-shell gluon polarizations is done using the gluon spin density matrix in the form \(\overline{ \epsilon ^\mu \epsilon ^{*\, \nu } } = {\mathbf k}_T^{\mu } {\mathbf k}_T^{\nu }/{\mathbf k}_T^2\), where \({\mathbf k}_T\) is the gluon transverse momentum orthogonal to the beam axis. In the collinear QCD limit, when \(|{\mathbf k}_T| \rightarrow 0\), this expression converges to the ordinary \(\overline{ \epsilon ^\mu \epsilon ^{*\, \nu } } = - g^{\mu \nu }/2\) after averaging over the azimuthal angle. Second, the \(\psi (2S)\) spin density matrix is expressed in terms of the momenta \(l_1\) and \(l_2\) of the decay leptons and is taken as

$$\begin{aligned} \sum \epsilon ^\mu \epsilon ^{*\, \nu } = 3 \left( l_1^\mu l_2^\nu + l_1^\nu l_2^\mu - {m^2\over 2} g^{\mu \nu } \right) /m^2. \end{aligned}$$
(14)

This expression is equivalent to the standard one \(\sum \epsilon ^\mu \epsilon ^{*\, \nu } = - g^{\mu \nu } + p^\mu p^\nu /m^2\) but is better suited for studying the polarization observables because it gives access to the kinematic variables describing the orientation of the decay plane. In all other respects the evaluation follows the standard QCD Feynman rules.

Table 1 The NMEs for \(\psi (2S)\) meson derived from the fit of the CMS data [21]. The NMEs obtained in the NLO NRQCD fits [15, 19] are shown for comparison

An important point in the NRQCD formalism is connected with the emission of soft gluons taking place after the hard interaction is over. It is usually assumed that the emitted soft gluons take away the unwanted color and change the other quantum numbers of the produced CO system but do not carry any energy, thus keeping the kinematics intact. However, such an emission contradicts the basic QCD property: soft gluons can never be radiated as they are confined. In order that the quantum numbers get changed, one needs to radiate a real gluon with some energy \(E \sim \Lambda _\mathrm{QCD}\), giving us the confidence that we do not enter into the confinement or perturbative domains. When considering the gluon radiation from the \(^3P_1^{(8)}\) and \(^3P_2^{(8)}\) states, we rely upon the dominance of the electric dipole E1 transitions, which is supported by the E835 experimental data [54]. The corresponding invariant amplitudes can be written as follows [55]:

$$\begin{aligned}&{\mathcal A}(^3P_1^{(8)} \rightarrow \psi ^\prime + g)=g_1 \, e^{\mu \nu \alpha \beta } k_\mu ^{(g)} \epsilon _\nu ^\mathrm{(CO)} \epsilon _\alpha ^{(\psi ^\prime )} \epsilon _\beta ^{(g)},\end{aligned}$$
(15)
$$\begin{aligned}&{\mathcal A}(^3P_2^{(8)} \rightarrow \psi ^\prime + g) \nonumber \\&\quad =g_2 \, p_\mu ^\mathrm{(CO)} \epsilon _{\alpha \beta }^\mathrm{(CO)} \epsilon _\alpha ^{(\psi ^\prime )} \left[ k_\mu ^{(g)} \epsilon _\beta ^{(g)} - k_\beta ^{(g)} \epsilon _\mu ^{(g)} \right] , \end{aligned}$$
(16)

where \(p_\mu ^\mathrm{(CO)}\), \(k_\mu ^{(g)}\), \(\epsilon _\mu ^{(\psi ^\prime )}\), \(\epsilon _\mu ^{(g)}\), \(\epsilon _\mu ^\mathrm{(CO)}\), and \(\epsilon _{\mu \nu }^\mathrm{(CO)}\) are the four-momenta and polarization four-vectors (tensor) of the corresponding particles and \(e^{\mu \nu \alpha \beta }\) is the fully antisymmetric Levi-Civita tensor. The gluon radiation from other CO states is generated according to the phase space.

The cross section of \(\psi (2S)\) production at high energies in the \(k_T\)-factorization approach is calculated as a convolution of the off-shell partonic cross sections and the TMD gluon densities in a proton. The contribution from the CS production mechanism can be presented in the following form:

$$\begin{aligned}&\displaystyle \sigma (p p \rightarrow \psi ^\prime + X) \nonumber \\&\quad = \int {1\over 16\pi (x_1 x_2 s)^2 } f_g(x_1,{\mathbf k}_{1T}^2,\mu ^2) f_g(x_2,{\mathbf k}_{2T}^2,\mu ^2) \nonumber \\&\qquad \displaystyle \times \, |\bar{\mathcal A}(g^* + g^* \rightarrow \psi ^\prime + g)|^2 \nonumber \\&\qquad \times \text {d}{\mathbf p}_{T}^2 \text {d}{\mathbf k}_{1T}^2 \text {d}{\mathbf k}_{2T}^2 \text {d}y \text {d}y_g \, {\text {d}\phi _1 \over 2\pi } {\text {d}\phi _2 \over 2\pi }, \end{aligned}$$
(17)

where \(f_g(x,{\mathbf k}_{T}^2,\mu ^2)\) is the TMD gluon density, \({\mathbf p}_T\) and y are the transverse momentum and rapidity of produced \(\psi (2S)\) meson, \(y_g\) is the rapidity of outgoing gluon, and \(\sqrt{s}\) is the pp center-of-mass energy. The initial off-shell gluons have a fraction \(x_1\) and \(x_2\) of the parent protons longitudinal momenta, non-zero transverse momenta \({\mathbf k}_{1T}\) and \({\mathbf k}_{2T}\) (\({\mathbf k}_{1T}^2 = - k_{1T}^2 \ne 0\), \({\mathbf k}_{2T}^2 = - k_{2T}^2 \ne 0\)), and azimuthal angles \(\phi _1\) and \(\phi _2\). For the CO production we have

$$\begin{aligned}&\displaystyle \sigma (p p \rightarrow \psi ^\prime + X) \nonumber \\&\quad =\int {2 \pi \over x_1 x_2 s \, F} \, f_g(x_1,{\mathbf k}_{1T}^2,\mu ^2) f_g(x_2,{\mathbf k}_{2T}^2,\mu ^2) \nonumber \\&\qquad \times \displaystyle |\bar{\mathcal A}(g^* + g^* \rightarrow \psi ^\prime )|^2 \, \text {d}{\mathbf k}_{1T}^2 \text {d}{\mathbf k}_{2T}^2 \text {d}y \nonumber \\&\qquad \times {\text {d}\phi _1 \over 2\pi } {\text {d}\phi _2 \over 2\pi }. \end{aligned}$$
(18)

According to the general definition [56], the off-shell gluon flux factor in (18) is definedFootnote 3 as \(F = 2 \lambda ^{1/2}(\hat{s},k_1^2,k_2^2)\), where \(\hat{s} = (k_1 + k_2)^2\). The squares of the corresponding off-shell partonic amplitudes, being too lengthy, are not presented, but the full C++ code is available on request.Footnote 4 The multidimensional integration has been performed by means of the Monte Carlo technique, using the routine vegas [57].

Numerically, we have tested several different sets of the TMD gluon densities. Two of them (namely, the JH and A0 sets) have been obtained [58, 59] from the CCFM equation where all input parameters have been fitted to describe the proton structure function \(F_2(x, Q^2)\). Besides the CCFM-evolved gluon densities, we applied the one obtained from the Kimber–Martin–Ryskin (KMR) prescription [60, 61]. The KMR approach is a formalism to construct the TMD quark and gluon distributions from well-known conventional ones. For the input, we have used the leading-order Martin–Stirling–Thorn–Watt (MSTW’2008) set [62].

3 Numerical results

We now are in a position to present our numerical results. First we describe our input and the kinematic conditions. Having the TMD gluon distributions chosen, the cross sections (17) and (18) depend on the renormalization and factorization scales \(\mu _R\) and \(\mu _F\). We set \(\mu _R^2 = m^2 + {\mathbf p}_{T}^2\) and \(\mu _F^2 = \hat{s} + {\mathbf Q}_T^2\), where \({\mathbf Q}_T\) is the transverse momentum of the initial off-shell gluon pair. The choice of \(\mu _R\) is the standard one for studying the charmonia production, whereas the special choice of \(\mu _F\) is connected with the CCFM evolution [58, 59]. Following [63], we set \(\psi (2S)\) mass \(m = 3.686\) GeV, branching fraction \(B(\psi ^\prime \rightarrow \mu ^+ \mu ^-) = 0.0077\) and use the LO formula for the coupling constant \(\alpha _s(\mu ^2)\) with \(n_f = 4\) quark flavors at \(\Lambda _\mathrm{QCD} = 200\) MeV, such that \(\alpha _s(M_Z^2) = 0.1232\).

Fig. 1
figure 1

The double differential cross sections of prompt \(\psi (2S)\) meson production at the LHC. Left panel: the dashed, dash-dotted, dotted, and short dash-dotted curves correspond to the color singlet \(^3S_1^{(1)}\) and color octet \(^1S_0^{(8)}\), \(^3S_1^{(8)}\), \(^3P_J^{(8)}\) contributions calculated with the KMR gluon density. The solid curve represent the sum of CS and CO terms. Right panel: the solid, dashed, and dash-dotted curves correspond to the predictions obtained with the A0, JH, and KMR gluon distributions, respectively. The experimental data are from ATLAS [20]

Fig. 2
figure 2

The transverse momentum distribution of prompt \(\psi (2S)\) meson production in pp collisions at the LHC. Notation of all curves is the same as in Fig. 1. The experimental data are from CMS [21]

Fig. 3
figure 3

The transverse momentum distribution of prompt \(\psi (2S)\) meson production in pp collisions at the LHC. Notation of all curves is the same as in Fig. 1. The experimental data are from LHCb [22]

Fig. 4
figure 4

Polarization parameters \(\lambda _\theta \), \(\lambda _\phi \), \(\lambda _{\theta \phi }\), and \(\lambda ^*\) of prompt \(\psi (2S)\) mesons calculated as a function of \(\psi (2S)\) transverse momentum in the Collins–Soper frame. The solid, dashed, and dash-dotted histograms correspond to the predictions obtained at \(|y| < 0.6\), \(0.6 < |y| < 1.2\), and \(1.2 < |y| < 1.5\). The KMR gluon distribution is used. The experimental data are from CMS [24]

Fig. 5
figure 5

Polarization parameters \(\lambda _\theta \), \(\lambda _\phi \), \(\lambda _{\theta \phi }\), and \(\lambda ^*\) of prompt \(\psi (2S)\) mesons calculated as a function of \(\psi (2S)\) transverse momentum in the helicity frame. The solid, dashed, and dash-dotted histograms correspond to the predictions obtained at \(|y| < 0.6\), \(0.6 < |y| < 1.2\), and \(1.2 < |y| < 1.5\). The KMR gluon distribution is used. The experimental data are from CMS [24]

Fig. 6
figure 6

Polarization parameters \(\lambda _\theta \), \(\lambda _\phi \), \(\lambda _{\theta \phi }\), and \(\lambda ^*\) of prompt \(\psi (2S)\) mesons calculated as a function of \(\psi (2S)\) transverse momentum in the perpendicular helicity frame. The solid, dashed, and dash-dotted histograms correspond to the predictions obtained at \(|y| < 0.6\), \(0.6 < |y| < 1.2\), and \(1.2 < |y| < 1.5\). The KMR gluon distribution is used. The experimental data are from CMS [24]

Fig. 7
figure 7

Polarization parameters \(\lambda _\theta \), \(\lambda _\phi \), \(\lambda _{\theta \phi }\), and \(\lambda ^*\) of prompt \(\psi (2S)\) mesons calculated as a function of \(\psi (2S)\) transverse momentum in the Collins–Soper frame. The solid, dashed, dash-dotted, dotted, and short dash-dotted histograms correspond to the predictions obtained at \(2 < y < 2.5\), \(2.5 < y < 3\), \(3 < y < 3.5\), \(3.5 < y < 4\), and \(4 < y < 4.5\). The KMR gluon distribution is used. The experimental data are from LHCb [25]

Fig. 8
figure 8

Polarization parameters \(\lambda _\theta \), \(\lambda _\phi \), \(\lambda _{\theta \phi }\), and \(\lambda ^*\) of prompt \(\psi (2S)\) mesons calculated as a function of \(\psi (2S)\) transverse momentum in the helicity frame. The solid, dashed, dash-dotted, dotted, and short dash-dotted histograms correspond to the predictions obtained at \(2 < y < 2.5\), \(2.5 < y < 3\), \(3 < y < 3.5\), \(3.5 < y < 4\), and \(4 < y < 4.5\). The KMR gluon distribution is used. The experimental data are from LHCb [25]

In Table 1 we list our results for the NMEs fits obtained for three different TMD gluon densities. We have fitted the transverse momentum distributions for prompt \(\psi (2S)\) mesons measured recently by the CMS Collaboration at \(\sqrt{s} = 7\) TeV [21]. These measurements were done at moderate and high transverse momenta \(10 < p_T < 100\) GeV, where the NRQCD formalism is believed to be most reliable. In contrast with [15, 19], we performed the fitting procedure under the requirement that the NME values be positive only. The color-singlet NMEs were not fitted, but just taken from the known \(\psi (2S) \rightarrow \mu ^+ \mu ^-\) partial decay width [63]. For comparison, we also present in Table 1 two sets of NMEs, obtained within the NLO NRQCD in [15, 19]. The main difference between [15, 19] is in that these fits were based on differently selected sets of data points.

In the \(k_T\)-factorization approach, the fitted NME values strongly depend on the choice of TMD gluon density. We find that the \(^1S_0^{(8)}\) contribution is compatible with zero if the CCFM-evolved gluon distributions are used, but it is non-negligible in the case of KMR gluons. For the latter, the extracted value of \(\left\langle {\mathcal O}^{\psi }\left[ ^1S_0^{(8)}\right] \right\rangle \) is very close to the one obtained in the NLO NRQCD analysis [15]. Both the NLO NRQCD fits [15, 19] significantly (by one order of magnitude) exceed the values of \(\left\langle {\mathcal O}^{\psi }\left[ ^3S_1^{(8)}\right] \right\rangle \), obtained with all of the TMD gluon densities. It is almost consistent with estimates performed by other authors [64, 65]. In contrast with the results [52, 53], our fit leads to non-zero \(\left\langle {\mathcal O}^{\psi }\left[ ^3P_0^{(8)}\right] \right\rangle \) values. Summing up, we can conclude that the NME values obtained by the different authors on the basis of different data sets or by the same authors using different gluon densities are widely spread, which spoils the belief in the universality of the matrix elements.

Now we turn to a comparison of our predictions with the data collected by the ATLAS [20], CMS [21] and LHCb [22] Collaborations. The ATLAS Collaboration has measured the prompt \(\psi (2S)\) transverse momentum distribution at \(10 < p_T < 100\) GeV at central rapidities \(|y| < 2\) [20]. The CMS Collaboration probed the transverse momentum in the range \(10 < p_T < 100\) GeV at \(|y| < 1.2\) [21], and the LHCb Collaboration worked in the kinematic range \(p_T < 16\) GeV and \(2 < y < 4.5\) [22]. In all cases the data were obtained at \(\sqrt{s} = 7\) TeV. The results of our calculations are shown in Figs. 1, 2 and 3. With our sets of NMEs, we achieve a reasonably good description of the data with any of the considered TMD gluon densities. We observe dominance of the CO contributions in the whole \(p_T\) range. In particular, the \(^3S_1^{(8)}\) contribution dominates at the large \(p_T > 25\) GeV, whereas the \(^1S_0^{(8)}\) channel is mostly important at low \(p_T\) values. Taken solely, the CS contributions (even incorporated with the \(k_T\)-factorization) are unable to describe the data. They are important at relatively low \(p_T\) only and are comparable there with the \(^3P_J^{(8)}\) contributions. At moderate and high transverse momenta the CS contributions are below the data by about one order of magnitude. The predictions obtained with the chosen TMD gluon densities are very close to each other at \(p_T > 6\) GeV, while the difference becomes only sizable at low \(p_T < 6\) GeV (see Fig. 3). Therefore, similar to the collinear QCD factorization, including the low \(p_T\) data to the fit procedure can change the relative weight of different NMEs in the \(k_T\)-factorization approach, which is beyond the scope of our present study.

Now we turn to the \(\psi (2S)\) polarization issue, which is the most interesting part of our study. In general, the spin density matrix of a vector particle decaying into a lepton pair depends on the three angular parameters \(\lambda _{\theta }\), \(\lambda _\phi \), and \(\lambda _{\theta \phi }\), which can be measured experimentally. The double differential angular distribution of the decay leptons can be written as [66]

$$\begin{aligned}&{\text {d}\sigma \over \text {d}\cos \theta ^* \text {d}\phi ^*} \sim 1 + \lambda _\theta \cos ^2 \theta ^* \nonumber \\&\quad + \lambda _\phi \sin ^2 \theta ^* \cos 2 \phi ^* + \lambda _{\theta \phi } \sin 2 \theta ^* \cos \phi ^*, \end{aligned}$$
(19)

where \(\theta ^*\) and \(\phi ^*\) are the polar and azimuthal angles of the decay lepton measured in the charmonium rest frame. The case of \((\lambda _{\theta }, \lambda _\phi , \lambda _{\theta \phi }) = (0,0,0)\) corresponds to the unpolarized state, while \((\lambda _{\theta }, \lambda _\phi , \lambda _{\theta \phi }) = (1,0,0)\) and \((\lambda _{\theta }, \lambda _\phi , \lambda _{\theta \phi }) = (-1,0,0)\) refer to the fully transverse and longitudinal polarizations. The CMS [24] and LHCb [25] Collaborations have measured all these parameters as functions of the \(\psi (2S)\) transverse momentum in two complementary frames: the Collins–Soper and helicity ones. In addition, the CMS Collaboration provided measurements in the perpendicular helicity frame. 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 \(\psi (2S)\) 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 momenta. Additionally, the frame-independent polarization parameter [6668] \(\lambda ^* = (\lambda _\theta + 3 \lambda _\phi )/(1 - \lambda _\phi )\) was investigated. Below we estimate the polarization parameters \(\lambda _\theta \), \(\lambda _\phi \), \(\lambda _{\theta \phi }\), and \(\lambda ^*\) for the CMS and LHCb conditions. Our calculation generally follows the experimental procedure. We collect the simulated events in the kinematical region defined by the CMS and LHCb experiments, generate the decay lepton angular distributions according to the production and decay matrix elements, and then apply a three-parametric fit based on (19).

In Figs. 4, 5, 6, 7 and 8 we confront our predictions for the polarization parameters \(\lambda _\theta \), \(\lambda _\phi \), \(\lambda _{\theta \phi }\), and \(\lambda ^*\) with the latest CMS [24] and LHCb [25] data. We find a slight transverse polarization (\(\lambda _\theta \sim 0.2\)) in the Collins–Soper frame and a slight longitudinal polarization (\(\lambda _\theta \sim - 0.2\)) in the helicity frame at low transverse momenta covered by the LHCb experiment. These results are practically independent on the \(\psi (2S)\) rapidity. At higher \(p_T\) the polarizations of \(\psi (2S)\) mesons, calculated in the Collins–Soper frame, go from slightly transverse (\(\lambda _\theta \sim 0.15\)) to almost zero values (\(\lambda _\theta \sim 0.05\)) as the transverse momentum increases from \(p_T \sim 10\) GeV to 50 GeV. In the helicity and perpendicular helicity frames, the \(\psi (2S)\) polarization changes from longitudinal (\(\lambda _\theta \sim - 0.2\)) to slightly longitudinal (\(\lambda _\theta \sim - 0.1\)). Here we arrive at the key point of our paper. Figures 4, 5, 6, 7 and 8 show that \(\psi (2S)\) production, calculated in the \(k_T\)-factorization approach, tends to be unpolarized at high \(p_T\), in agreement with the CMS data [24]. Indeed, as a strict consequence of the initial gluon off-shellness, a large fraction of the \(\psi (2S)\) mesons with zero helicity is produced in the partonic subprocesses, including both CS and CO contributions. Moreover, the fraction of such events increases when \(p_T\) increases (see, for example, [3340]). The only exception is in the \(^1S_0^{(8)}\) channel, which is produced unpolarized due to its spinless nature. It is a remarkable property of the \(k_T\)-factorization scheme that the gluon fragmentation to \(^3S_1^{(8)}\) states produces nearly unpolarized \(\psi (2S)\) mesons (in contrast with conventional collinear NRQCD where the mesons carry strong transverse polarization). Thus, we can conclude that the problem of \(\psi (2S)\) spin alignment can be solved if the initial gluon off-shellness is taken into account.

A comparison of our predictions with the LHC data [24, 25] shows that the latter seem to support the trend observed in the \(k_T\)-factorization formalism. However, while our predictions for the \(\lambda _\phi \) and \(\lambda _{\theta \phi }\) parameters agree with the data, the description of \(\lambda _{\theta }\) and \(\lambda ^*\) is still qualitative rather than quantitative, due to the huge experimental uncertainties.

Finally, we would like to note that there are significant theoretical uncertainties connected with the choice of the renormalization and/or factorization scales, the inclusion of NLO subprocesses, and the exact definition of NMEs. The detailed study of these uncertainties is beyond the scope of our present paper.

4 Conclusions

We have considered prompt \(\psi (2S)\) production and polarization in pp collisions at the LHC energy \(\sqrt{s} = 7\) TeV in the framework of \(k_T\)-factorization approach. We have used the LO non-relativistic QCD formalism including both color-singlet and color-octet contributions. Using the TMD gluon densities in a proton derived from the CCFM equation and from the Kimber–Martin–Ryskin prescription, we extracted the color-octet NMEs \(\left\langle {\mathcal O}^{\psi }\left[ ^1S_0^{(8)}\right] \right\rangle \), \(\left\langle {\mathcal O}^{\psi }\left[ ^3S_1^{(8)}\right] \right\rangle \) and \(\left\langle {\mathcal O}^{\psi }\left[ ^3P_0^{(8)}\right] \right\rangle \) for the \(\psi (2S)\) mesons from fits to transverse momentum distributions provided by the latest CMS measurements. Using the fitted NMEs, we have successfully described the data presented by the ATLAS, CMS, and LHCb Collaborations. We estimated the polarization parameters \(\lambda _\theta \), \(\lambda _\phi \), and \(\lambda _{\theta \phi }\), and we demonstrated that taking into account the off-shellness of the initial gluons in the color-octet contributions leads to unpolarized \(\psi (2S)\) production at high transverse momenta, in qualitative agreement with the LHC data.