Prompt hadroproduction of $\eta_c(1S,2S)$ in the $k_T$-factorization approach

In this work, we present a thorough analysis of $\eta_c(1S,2S)$ quarkonia hadroproduction in $k_\perp$-factorisation in the framework of the light-front potential approach for the quarkonium wave function. The off-shell matrix elements for the $g^* g^* \eta_c(1S,2S)$ vertices are derived. We discuss the importance of taking into account the gluon virtualities. %We discuss the role of relativistic corrections, the Melosh spin-rotation and theoretical uncertainties by comparing our results with the standard NRQCD approach. We present the transverse momentum distributions of $\eta_c$ for several models of the unintegrated gluon distributions. Our calculations are performed for four distinct parameterisations for the $c\bar c$ interaction potential consistent with the meson spectra. We compare our results for $\eta_c(1S)$ to measurements by the LHCb collaboration and present predictions for $\eta_c(2S)$ production.


I. INTRODUCTION
present the most relevant numerical results for the differential η c (1S, 2S) production cross sections versus the available experimental data and discuss the related theoretical uncertainties. The basic concluding remarks and the main results are summarised in Section IV.

II. FORMALISM
A. Off-shell matrix element and cross section In Fig.1 we show a generic Feynman diagram for η c (1S) quarkonium production in protonproton collision via gluon-gluon fusion. This diagram illustrates the situation adequate for the k T -factorization calculations used in the present paper. The inclusive cross section for η c -production via the 2 → 1 gluon-gluon fusion mode is obtained from (2.1) The unintegrated gluon distributions are normalized such, that the collinear glue is obtained from where from now on we will no longer show the dependence on the factorization scale µ 2 where m is the mass of the η c -meson, and y is its rapidity in the pp cms-frame. The phasespace element is dΦ(2 → 1) = (2π) 4 δ (4) (q 1 + q 2 − p) d 4 p (2π) 3 δ(p 2 − m 2 ) . (2.5) In the k T -factorization approach, gluons are off-shell, q 2 i = −q 2 i , and their four momenta are written as ( √ s is the pp center-of-mass energy): with (2.7) We can then calculate the phase-space element as We therefore obtain for the inclusive cross section where the momentum fractions of gluons are fixed as x 1,2 = m T exp(±y)/ √ s. The off-shell matrix element is written in terms of the Feynman amplitude as (we restore the colorindices): In covariant form, the matrix element reads: To the lowest order, it is proportional to the matrix element for the γ * γ * η c vertex. In particular, the form factor and it can be represented in terms of the LF wave function as For details of the derivation and the normalization conventions and relation to the potential model wave function of the LF radial wave function ψ(z, k), see Ref. [9]. In this work, we will use the calculations of the form factor which were obtained in [9]. There, the representation of the γ * γ * η c transition form factor in terms of the LF wave function of the η c was derived. Several wave-functions obtained from potential models for the cc system which were previously obtained in Ref. [10] were used. Inserting the explicit form of the matrix element and averaging over colors, we obtain our final result: Parametrizing the transverse momenta as In our numerical calculations presented below, we set the factorization scale to µ 2 F = m 2 T , and the renormalization scale is taken in the form: The normalization of the inclusive cross section depends crucially on the value of the g * g * η c form factor for vanishing gluon virtualities q 2 1 = q 2 2 = 0. The latter in turn is directly related to the η c → gg decay width. From the proportionality of the g * g * η c and γ * γ * η c vertices to the leading order (LO), we obtain, that at LO, the γγ and gg widths are related by where the LO γγ width in turn is related to the transition form factor for vanishing virtualities through At NLO, the expressions for the widths read (see e.g. [21]) (2.20) In order to control the model uncertainty on the normalization, one may want to adjust its value F (0, 0) to the measured decay width. Here we face the ambiguity of fitting either to the hadronic or to the γγ width. As there are no other known radiative decays besides γγ, one may try to identify the gg-width with the total (hadronic) width.
In Tables 1,2, we show the values of |F (0, 0)| obtained in three different ways. In Table 1 we show the result extracted from the total decay width. Here we use the strong coupling α s = 0.26, which is appropriate to our choice of the renormalization scale in the production amplitudes. In Table 2 we extract the value of |F (0, 0)| from the radiative decay width in two different ways. The first result is obtained based on Eq.(2.19) using the experimental value for Γ(η c → γγ) on the left hand side, while the second one uses the NLO relation 2.20.
We observe a substantial difference between the two different extractions of |F (0, 0)|. While in the η c (2S) case, the error bars are too large to claim an inconsistency, the situation for the η c (1S) is not satisfactory. This is in fact an old problem and may hint at an insufficiency of the potential model treatment of the η c . Various possible solutions have been proposed, such as an admixture of light hadron states [22], a mixing with a pseudoscalar glueball [23], or nonperturbative instanton effects in the hadronic decay [24].

C. Unintegrated gluon distributions
We use a few different UGDs which are available from the literature, e.g. from the TMDLib package [26] or the CASCADE Monte Carlo code [27].
1. Firstly we use a glue constructed according to the prescription initiated in [28] and later updated in [29], which we label below as "KMR". It uses as an input the collinear gluon distribution from [30].
2. Secondly, we employ two UGDs obtained by Kutak in [31]. There are two versions of this UGD. Both introduce a hard scale dependence via a Sudakov form factor into solutions of a small-x evolution equation. The first version uses the solution of a linear, BFKL [5] evolution with a resummation of subleading terms and is denoted by "Kutak (linear)". The second UGD, denoted as "Kutak (nonlinear)" uses instead a nonlinear evolution equation of Balitsky-Kovchegov [32] type. Both of the Kutak's UGDs can be applied only in the small-x regime, x < 0.01.
3. The third type of UGD which we use has been obtained by Hautmann and Jung [33] from a description of precise HERA data on deep inelastic structure function by a solution of the CCFM evolution equations [34]. We use "Set 2" of Ref. [33].
For the case of the KMR UGD, it has recently been shown in [35], that it includes effectively also higher order corrections of the collinear factorization approach. In this sense should give within our approach a result similar to that found recently in the NLO approach [16] at not too small transverse momenta. In contrast to the collinear NLO approach in our approach we can go to very small transverse momenta close to p T = 0.

III. NUMERICAL RESULTS
Before presenting results for the cross sections let us understand first the kinematical situation relevant for the LHCb experiment. The gluons entering the g * g * → χ c vertex (see Fig.1) are characterized by their longitudinal fractions (x 1 or x 2 ) or gluon virtualities at high-energies directly related to the their transverse momenta q 1,2 . In Fig.2 we show two dimensional distributions for the first gluon (x 1 , q 1T = |q 1 |) (left panel) and for the second gluon (x 2 , q 2T = |q 2 |) (right panel). We observe a large asymmetry of the two distributions related to asymmetric LHCb configuration: y ∈ (2, 4.5). The relatively large lower cut on η c transverse momentum p T > 6.5 GeV causes that q 1T and q 2T are themselves not small and are essentially in the perturbative regime where UGDs should be rather reliable. Moreover, large p T of the η c also entails the large factorization scale. In this calculation we used the KMR UGD described above.
The projections on the x i and q iT axes are shown in Fig.3. The asymmetric LHCb kinematics causes that x 1 is rather large and x 2 is very small -even smaller than 10 −5 , much smaller than for other perturbative partonic processes. We observe also a clear asymmetry in q 1T and q 2T . The q 2T transverse momentum corresponding to small x 2 is substantially larger than q 1T of the large-x 1 gluon. The low-x gluon therefore transfers the bulk of the  transverse momentum of the η c at large p T . We can exploit the good separation in x 1,2 to investigate the small-x behaviour of the unintegrated glue. From our choice of UGDs, the parametrizations of Kutak are available only for x < 0.01, so we will use these UGDs only for the small-x gluon. To avoid a proliferation of plots, we will use the Kutak UGDs always together with the KMR UGD for the large-x gluon. A similar strategy was taken previously in Ref. [6].
The distributions for different UGDs in Fig. 3 are rather similar which makes are conclusions more universal. Now we wish to show the behaviour of the different unintegrated gluon distribution on gluon transverse momentum k 2 T for small x 2 = 10 −5 and typical scale parameter µ 2 = 100 GeV 2 relevant for the LHCb experiment. In the left panel of Fig. 4 we plot F (x, k 2 , µ 2 F )/k 2 , which corresponds to the distribution of the gluon transverse momentum squared k 2 . In the right panel the dimensionless UGD F (x, k 2 , µ 2 F ) is plotted. We show all the UGDs used in the present work. The left panel of Fig. 4 better shows the behaviour at smaller k 2 = k 2 T , while the right panel emphasizes the large-k T tails. We first observe, that the "linear" Kutak UGD looks quite similar to the KMR UGD, although both are constructed by different procedures. They have in common, that by their construction, both procedures lead to integrated gluon distributions which well describe jet cross sections at the LHC. The nonlinear Kutak UGD is considerably smaller than the linear one, especially at low transverse momenta. At very large transverse momenta the difference between linear and nonlinear UGDs becomes much smaller. The Jung-Hautmann distribution does not have an extended tail in k 2 T as the other distributions, it is however much larger at low gluon transverse momenta. Can these different UGDs be tested in η c production?
In Fig. 5 we compare the results of the k T -factorization approach with the LHCb experimental data for √ s = 7 and 8 TeV from Ref. [7] and with the data for √ s = 13 TeV from the recent PhD thesis [8]. The theoretical calculations use an off-shell form factor normalized to the γγ-decay width at NLO. The off-shell g * g * η c form factor was calculated using the LF wave function obtained for one of the potentials (the so-called power-law potential) in [10]. It is up to the color factor proportional to the γ * γ * η c form factor obtained in [9]. The dependence on the choice of the potential will be discussed below. The description of data for √ s = 7, 8 TeV is reasonable for all UGDs. The theoretical results tend to be somewhat lower than the experimental data, especially at large p T . The best description is obtained for the KMR UGD and the linear UGD by Kutak.
Our calculations fare a bit worse in the comparison to the data at √ s = 13 TeV. Here all the UGDs give results substantially below the data. Please note, that we include only the direct production mechanism in the color-singlet channel. A possible feed down from higher resonances, e.g. from the h c → η c γ radiative decay, is not included. See the recent Ref. [19] for an estimate which finds a few percent contribution from the feed down. We also do not include a possible color-octet contribution. In general the difference between different UGDs is the largest at low p T . Here one may e.g. expect effects related to nonlinear evolution and gluon saturation. While it is true that the cross section peaks precisely in this most interesting region around p T ∼ 2 GeV it seems to be exceedingly difficult to measure at p T < 6 GeV. At least this is true in the pp decay channel used by the LHCb collaboration. Perhaps other decay channels would be better in this respect. The γγ channel seems interesting as another option. A simulation of the signal and background would be valuable in this context.
In Fig.6 we show our predictions for the so far unmeasured (at the LHC) η c (2S). In this calculation we also adjusted the value of |F (0, 0)| to the γγ width at NLO. The cross section is in a similar ballpark as for η c (1S) and the results for different UGDs show a similar variation as for the η c (1S) production.
We now wish to turn to the dependence of our results on the cc potential used for the calculation of the LF wave function. In Ref. [9] the γ * γ * → η c (1S, 2S) transition form factors, which is closely related to the g * g * η c (1S, 2S) form factors were obtained for different cc potentials. In Fig.7 [7] for √ s = 7, 8 TeV and preliminary experimental data [8] for √ s = 13 TeV. Different UGDs were used. Here we used the g * g * → η c (1S) form factor calculated from the power-law potential. dictated by the NLO expression for the experimental η c (1S, 2S) → γγ decay width.
In Fig.8 we relax the normalization to the values predicted by the different potentials. Here the spread of results is bigger than in the previous case. We need to caution the reader, that generally the results from the phenomenological potentials undershoot the experimental widths. However, we wish to notice that experimental decay widths are known only with some precision [25]. In [36] a different values was measured. It appears, that the behaviour of the off-shell form factors at large q 2 i is a more reliable result than their value at the on-shell point, though.
In Fig.9 we show the integrated cross section within the LHCb cuts for the three different energies measured up to now. In this plot we used the KMR UGD. We observe that the data appear to indicate a faster than linear rise of the cross section, while the theoretical calculations predict a slightly slower rise of the cross section with √ s. We checked that similar results are obtained with other UGDs. Once more we remind of the missing feeddown component. A possible color-octet component is not expected to change the energy dependence, as it is driven by the similar g * g * → cc process. It is interesting to investigate what is the role of the off-shell form factor. For example in the approach of Ref. [17] gluon virtualities are neglected in the hard matrix element. The curves in Fig. 10 clearly show that the effect of the inclusion of gluon virtualities in the transition form factor is essential and cannot be neglected.
Up to now we concentrated on the kinematics of the LHCb experiment. Could η c quarkonia be measured by other experiments at the LHC? In Fig.11 we show the ranges of x i and q T i carried by gluons for a rapidity interval −2.5 < y < 2.5 typical of the central detectors of ATLAS or CMS. As in the case of LHCb, we assumed a lower cut on transverse momentum, p T > 6.5 GeV. In the center-of-mass rapidity interval symmetric around zero, of course both UGDs enter symmetrically and therefore the distributions of x 1 and x 2 coincide, as do the ones for q 1T and q 2T We find that in the central rapidity region one would test x 1 , x 2 ∼ 10 −4 -10 −2 i.e. the region where the gluon is already known reasonably well from the HERA experiments. The q iT -distribution has a large plateau at perturbatively large values, so that we suppose that the predictions of the KMR UGD should be reliable in this case.
The corresponding transverse momentum distributions are shown in Fig.12 for two differ-

IV. CONCLUSIONS
In the present paper we have discussed in detail the production of pseudoscalar η c quarkonia in proton-proton collisions at the LHC energies. The calculations have been performed in the k T -factorization approach using g * g * → η c vertex calculated from the cc wave functions. The latter vertices are closely related to the γ * γ * → η c transition form factors which were obtained recently in [9] and which we used in the present work.
We have used different unintegrated gluon distributions available in the literature: firstly the KMR UGD, which effectively includes higher-order effects of the collinear approach, secondly two UGDs by Kutak which were obtained by using linear and nonlinear small-x evolution equations, and thirdly the Jung-Hautmann UGD from 2013, obtained from a fit to HERA data in a CCFM approach.
We have calculated transverse momentum distributions of both η c (1S) and η c (2S) charmonia. The results of the η c (1S) have been compared with the experimental data obtained by the LHCb collaboration for √ s = 7, 8 GeV [7] and to a measurement at √ s = 13 TeV published in a PhD thesis [8].
A quite good agreement with the data was obtained with the KMR UGD. We have shown the range of x 1 , x 2 and gluon transverse momenta q 1T , q 2T probed in the kinematics of the LHCb experiment [7,8]. For the LHCb experiment one of the x-values, x 1 ∈ (10 −2 , 10 −1 ) is large, while the second one, x 2 ∈ (10 −5 , 10 −4 ) takes very small values. For the LHCb experiment we have shown also a large asymmetry in q 1T (larger) and q 2T (smaller). It turns out that at large p T of the meson the bulk of transverse momenta is transferred by the small-x gluon. The Kutak UGD cannot be used for the range of x 1 relevant for the LHCb experiment. Therefore in this case we have used two different UGDs: Kutak UGD for small x 2 and KMR UGD for large x 1 . We have used both the linear and nonlinear UGDs of Kutak. The mixed UGD scenario with the linear Kutak UGD leads to very similar transverse momentum distribution of η c as that for using the KMR UGD on both sides. Indeed in the considered range of x and k T the linear Kutak UGD is very similar as the KMR UGD. The nonlinear version of the Kutak UGD leads to smaller cross sections, especially for small η c transverse momenta.
A measurement of η c at low transverse momenta in the LHCb kinematics would therefore be very valuable in the context of searching for nonlinear effects and onset of gluon saturation. One could, for example, consider another measurement for the γγ decay channel.
We have shown that it is crucial to include the dependence on gluon virtualities in the g * g * → η c (1S, 2S) vertex. We have also discussed uncertainties related to the g * g * → η c (1S,2S) form factor. We have shown results of calculations with the form factor obtained from different cc potentials from the literature. The associated uncertainty is somewhat smaller than that related to the choice of UGD. A better measurement of η c (1S,2S) → γγ would be important in this context.