Prompt inclusive production of $J/\psi$, $\psi'$ and $\chi_{c}$ mesons at the LHC in forward directions within the NRQCD $k_t$-factorization approach - search for the onset of gluon saturation

We discuss prompt production of $J/\psi$ mesons in proton-proton collisions at the LHC within NRQCD $k_t$-factorization approach using Kimber-Martin-Ryskin (KMR) unintegrated gluon distributions (UGDF). We include both direct color-singlet production ($g g \to J/\psi g$) as well as a feed-down from $\chi_c \to J/\psi \gamma$ and $\psi' \to J/\psi X$ decays. The production of the decaying mesons ($\chi_c$ or $\psi'$) is also calculated within NRQCD $k_t$-factorization approach. The corresponding matrix elements for $g g \to J/\psi$ g, $g g \to \psi'g$ and $g g \to \chi_c$ include parameters of the nonrelativistic space wave functions of the quarkonia at $r = 0$, which are taken from potential models in the literature. We get the ratio of the corresponding cross section ratio for $\chi_c(2)$-to-$\chi_c(1)$ at midrapidities much closer to experimental data than obtained in a recent analysis. Differential distributions in rapidity and transverse momentum of $J/\psi$ and $\psi'$ are calculated and compared with experimental data of the ALICE and LHCb collaborations. We discuss a possible onset of gluon saturation effects in the production of $J/\psi$ and $\chi_{c}$ mesons at forward/backward rapidities. We show that it is neccessary to modify the standard KMR UGDF to describe ALICE and LHCb data. A mixed UGDF scenario was proposed. Then we can describe the experimental data for $J/\psi$ production within model uncertainties with color-singlet component only. Therefore our theoretical results leave only a relatively small room for the color-octet contributions.


I. INTRODUCTION
There is a long-standing lack of convergence in understanding production of J/ψ quarkonia in proton-proton or proton-antiproton collisions. Some authors believe that the corresponding cross sections are dominated by the so-called color-octet contribution. On the other hand, some other authors expect that the color-singlet contribution dominates. The color-octet contribution cannot be calculated from first principle and is rather fitted to the experimental data. The fits lead to different sizes of the color-octet contribution, depending on the details of calculations of the color-singlet contribution(s). In many cases, successful fits were obtained, but in our opinion, there is no clear understanding of the problem. Different fits from the literature give different magnitudes of the color-octet contributions classified according to quantum numbers of the cc system.
In the present paper, we calculate the color-singlet contributions in the NRQCD k t -factorization approach and see how much room is left for the more difficult color-octet contribution.
It is known that a sizeable part of the J/ψ production comes from radiative decays of χ c mesons. Therefore, in the following, we have to include also this contribution very carefully trying to confront with experimental data for χ c production whenever possible.
In a very recent k t -factorization analysis of χ c production [1], the authors found very different values of the nonrelativistic wave function at the origin for χ c ð1Þ and χ c ð2Þ, jR 0 χ c ð1Þ ð0Þj 2 ≈ 5jR 0 χ c ð2Þ ð0Þj 2 ; ð1:1Þ from a fit based on the k t -factorization approach to LHC data. This large modification would put in doubt either the NRQCD approach and/or validity of the leading-order k t -factorization. In the standard potential model, one obtains the same radial wave function for different χ c species [2]. Here, we diccuss also this element of the whole construction. In the following, we use Kniehl-Vasin-Saleev matrix elements which are given explicitly in Ref. [3]. Finally, the ψ 0 quarkonium also has a sizeable branching fraction into J/ψX [4]. Fortunately, this contribution is much smaller than the direct one as is discussed in this paper. It was considered recently in almost identical approach in [5].
In the present approach, we concentrate rather on small transverse momenta of J/ψ or ψ 0 relevant for ALICE and LHCb data [6][7][8][9][10]. We expect that color-singlet contributions may dominate in this region of the phase space.

II. SOME THEORETICAL ASPECTS
In the following, we consider only color-singlet mechanisms and look at how much room is left for color-octet production.

A. Main contributions
The main color-singlet mechanism of the J/ψ meson production is illustrated in Fig. 1. In this case, J/ψ is produced in association with an extra "hard" gluon due to C-parity conservation.
We calculate the dominant color-singlet gg → J/ψg contribution taking into account transverse momenta of initial gluons. In the k t -factorization approach, the differential cross section can be written as dσðpp → J/ψgXÞ dy J/ψ dy g d 2 p J/ψ;t d 2 p g;t π jM off-shell g Ã g Ã →J/ψg j 2 × δ 2 ð⃗ q 1t þ ⃗ q 2t − ⃗p H;t − ⃗p g;t Þ × F g ðx 1 ; q 2 1t ; μ 2 ÞF g ðx 2 ; q 2 2t ; μ 2 Þ; ð2:1Þ where F g are unintegrated (or transverse-momentumdependent) gluon distributions. The matrix elements were calculated as was explained, for example, in [11,12]. The corresponding matrix element squared for the gg → J/ψg is jM gg→J/ψg j 2 ∝ α 3 s jRð0Þj 2 : ð2:2Þ Running coupling constants are used in the present calculation. Different combinations of renormalization scales were tried. We decided to use The factorization scale in the calculation was taken as μ 2 F ¼ ðm 2 t þ p 2 t;g Þ/2. Similarly, we do calculations for the P-wave χ c meson production. Here, the lowest-order subprocess gg → χ c is allowed by positive C-parity of χ c mesons. In the k t -factorization approach, the leading-order cross section for the χ c meson production can be written somewhat formally as where F g are unintegrated (or transverse-momentumdependent) gluon distributions and σ gg→χ c are gg → χ c (off-shell) cross sections. The situation is illustrated diagramatically in Fig. 2. The matrix element squared for the gg → χ c subprocess is jM gg→χ c j 2 ∝ α 2 s jR 0 ð0Þj 2 : ð2:5Þ FIG. 1. Leading-order diagram for direct J/ψ (ψ 0 ) meson production in the k t -factorization approach.
FIG. 2. Leading-order diagram for χ c meson production in the k t -factorization approach.
After some manipulation, that can be also used to calculate rapidity and transverse momentum distribution of the χ c mesons. In the last equation, ⃗p t ¼ ⃗ q 1t þ ⃗ q 2t is the transverse momentum of the χ c meson and ⃗ q t ¼ ⃗ q 1t − ⃗ q 2t is the auxiliary variable which is used for the integration of the cross section. Furthermore, m t;χ c is the so-called χ c transverse mass and x 1 ¼ 1 4 is the Jacobian of transformation from ð⃗ q 1t ; ⃗ q 2t Þ to ð ⃗p t ; ⃗ q t Þ variables. As for the J/ψ production, running coupling constants are used. Different combination of scales were tried. The best choices are The factorization scale(s) for the χ c meson production are fixed traditionally as μ 2 F ¼ m 2 t . The J/ψ mesons are produced then by the χ c → J/ψγ decays which are dominated by E1 transitions [13,14]. This channel cannot be easily eliminated experimentally as the produced photons are usually rather soft. Due to the same reasons, χ c mesons can be measured at large transverse momenta or very forward/backward directions.

B. Unintegrated gluon distributions
In the present analysis, the Kimber-Martin-Ryskin KMR UGDFs [15] are used, which are generated from conventional collinear MMTH2014LO PDFs [16]. In actual calculations of distributions, we interpolate them on a three-dimensional grid in log 10 ðxÞ, log 10 ðk 2 t Þ and log 10 ðμ 2 Þ prepared before the calculation of the production cross section or differential distributions.
The KMR UGDF was succesfully used, for example, for production of charm and charmed mesons [17,18] as well as for production of two pairs of cc [19,20].
In a standard approach, the KMR UGDFs are calculated for larger values of gluon transverse momenta and are usually frozen at small gluon transverse momenta. The value at which the freezing is applied is independent of all other variables, longitudinal momentum fraction, in particular. The UGDFs used in calculations neglect possible effects of saturation. For small initial gluon momenta, k 2 1t < Q 2 s or k 2 2t < Q 2 s , and for forward/backward production, some effects of gluon saturation may be expected. The saturation scale as is often parametrized as One could correct the original KMR distributions by assuming saturation of UGDFs for k 2 it < Q 2 s , We call this model of UGDF "saturation A" for brevity. For comparison, we consider also faster damping of the smallk t region by multiplying the F A by an extra damping factor, We call this model of UGDF "saturation B" for brevity. Another, called "mixed UGDF", scenario is discussed in the text. Some consequences of the small-k t corrections are discussed in the following section.

III. RESULTS
A. ψ 0 production We start with ψ 0 production. In Fig. 3 For very small x, one may expect saturation effects. The KMR UGDF does not include such effects. To illustrate the potential effect instead of using both UGDFs of the same type (KMR) for large x 1 /x 2 , we take the KMR UGDF, whereas for small x 2 /x 1 , we take, for example, the Kutak-Staśto nonlinear UGDF [21]. It is marked by the KS acronym in Fig. 4. A slightly smaller cross section has been obtained than with the KMR UGDF. The effect is, however, not significant.
Since the ψ 0 meson decays ψ 0 → J/ψX with BF ¼ 0.61 [4], it constitutes also a contribution to the J/ψ channel and is taken into account in the rest of the paper.

B. J/Ψ direct production
There are three components of prompt J/ψ production: direct production (see Figs. 5-8) and feed down from ψ 0 and χ c decays.
In this subsection, we present results for the direct component for J/ψ production. In Fig. 5, we show exclusively this contribution for three different collision energies: and ffiffi ffi s p ¼ 13 TeV (right panel). As for ψ 0 production, we show our results for two different prescriptions for α s and for the KMR UGDF. As for the ψ 0 production, there is large uncertainty related to the choice of running coupling constant (see the yellow band). The direct contribution is large, but there is a room for other contributions, which are discussed in the following.   The results are compared with the ALICE [6][7][8] and LHCb [9,10] experimental data.
For completeness in Fig. 6, we show also our results with the mixed UGDFs scenario described above. Here, the effect of UGDF modification is similar as for ψ 0 .
We look also for transverse momentum distributions. For example, the LHCb Collaboration measured such distributions for different intervals of rapidity [9,10]. In Fig. 7, we show such distributions for ffiffi ffi s p ¼ 7 TeV for  two different prescriptions of α s . Our direct component exhausts a large fraction of the cross section for small p t . At larger p t , clearly some other contributions are missing.
In Fig. 8, we show similar distributions for ffiffi ffi The situation is very much the same as for ffiffi ffi Now we proceed to the important contribution of J/ψ originating from the feed down from the χ c production and decay. The χ c ð0Þ meson has a very small branching fraction for decay χ c ð0Þ → J/ψγ (BR ¼ 0.0127 [4]). Therefore, in the following, we take into account only production and decays of χ c ð1Þ (BF ¼ 0.339 [4]) and χ c ð2Þ (BF ¼ 0.192 [4]). In Fig. 9, we show rapidity distributions of resulting J/ψ mesons for two different prescriptions of α s (compare top and bottom panels) for three different energies ffiffi ffi s p ¼ 2.76, 7, 13 TeV. This calculations were performed with the KMR UGDF. For comparison, we present also existing data of the ALICE and LHCb Collaborations. We show both contributions of each of the mesons and a sum of them. The first prescription leads to a clearly too large cross section, having in mind the other missing contributions. This is especially clearly seen for ffiffi ffi s p ¼ 13 TeV, at large rapidities. The second prescription is not in conflict with the data, but one should remember other, not yet included, contributions (direct one and ψ 0 feed down). The situation with χ c production seems more problematic than for the direct contribution and not yet discussed ψ 0 feed down. What is specific for χ c production? In Fig. 10 of searching for saturation/nonlinear effects. We show results when the averaging is performed in different regions of χ c transverse momenta.
The very small values of longitudinal momentum fractions relevant for χ c production in the forward directions fully justify the use of the "mixed" UGDFs, discussed already in the context of direct production. In Fig. 11, we show corresponding rapidity distributions. The results obtained for the "mixed" distributions are quite different than those obtained solely with the KMR UGDFs, especially for ffiffi ffi s p ¼ 13 TeV. Is it a sign of the onset of saturation?
This should be clarified in the future by dedicated measurements of χ c mesons for different rapidities. This process seems to be very promising in this context.  D. χ c production So far χ c mesons were measured only at ffiffi ffi s p ¼ 7 TeV, at midrapidities and rather large transverse momenta. Then, the corresponding longitudinal momentum fractions are not so small. In Fig. 12, we show our results for both ffiffi ffi s p ¼ 7 TeV, together with ATLAS experimental data [22], and our predictions for ffiffi ffi s p ¼ 13 TeV. We get a reasonable, but not ideal, description of the experimental transverse momentum distributions for χ c ð1Þ (left panel) and χ c ð2Þ (middle panel). We slightly overestimate the data for χ c ð2Þ especially for smaller values of transverse momenta. For completeness, we show the ratio χ c ð2Þ/ χ c ð1Þ. In principle, we could try to treat parameters of χ c ð1Þ and χ c ð2Þ independently and better fit them to the ATLAS data, but we leave it for the future when next-to-leading order corrections will be included. Summarizing this short subsection, we have shown that our parameters for χ c are reasonable. They are to some extent effective as only leading-order k t -factorization is done here. How it changes at next-to-leading order clearly goes beyond the scope of the present analysis.

E. All contributions for J/ψ production
Having reviewed all components separately, we are ready to include all of them together. In the following, we adopt always prescription 2 (lower limits above) for α s as an example.
In Fig. 13 The disagreement is larger for larger rapidities (smaller longitudinal momentum fractions). This may be related to onset of saturation in this region of phase space and is worthy of further study. In Fig. 14, we show similar results for the "mixed" scenario. We get too much damping of the cross section, especially for largest ffiffi ffi s p . This may signal also the presence of other, nonincluded, mechanisms or may signal that the KS saturation effects are too strong. They may also appear too early in x.
Since, as discussed above, the longitudinal momentum fractions for J/ψ and ψ 0 are about an order of magnitude larger than those for χ c production, we consider also a new scenario. Here, we take standard KMR UGDFs for the Swave quarkonia and "mixed" UGDFs for the χ c mesons. The resulting distributions are shown in Fig. 15. The agreement with the experimental data is very good, but we cannot draw too strong conclusions. More systematic      studies of low-p t distributions of J/ψ and χ c mesons would certainly be very useful in this context.

F. General comments and relations to other
approaches in the literature In very forward directions, often a so-called hybrid approach is applied for different reactions, like forward jet or dijet production. In this approach, one uses one collinear parton distribution and one unintegrated gluon distribution. The hybrid model is claimed to be a sensible approximation in very forward directions [large x F of the produced object(s)], when one of the longitudinal fractions (x 1 or x 2 ) is very small and the second longitudinal fraction is large. At the LHCb or forward ALICE measurements for x's are rather small x 1 , x 2 < 10 −1 (see Fig. 10). For the LHCb rapidity coverage, x F < 0.05. Therefore, for the high energy collisions and so-called "forward" LHCb rapidity region, this is still a rather "central" rapidity region in the sense of the longitudinal momentum fractions. The corresponding larger gluon longitudinal fractions (maxfx 1 ; x 2 g) are similar as those for central rapidities at RHIC.   15. Rapidity distribution of J/ψ mesons for all considered mechanisms for three different energies. In this calculation, only for χ c production, the mixed UGDFs scenario was used.  We showed in our previous works that the inclusive D meson production at the LHCb can be described nicely within the original k t -factorization, see, e.g., [19]. We checked that the hybrid approach is not the best there and leads to a deficit of the cross section for ffiffi ffi s p ¼ 7 and 8 TeV. Now, we show some results of the hybrid model and a comparison to the results obtained within the original k t -factorization approach. In Fig. 16, for example, we show rapidity distribution of χ c ð1Þ mesons (2 → 1 partonic process). The hybrid model gives almost the same result as the original k t -factorization approach. In Fig. 17, we show corresponding transverse momentum distribution obtained by integrating over the LHCb rapidity region 2.0 < y < 4.5. Also, here, the agreement is quite good. Finally, in Fig. 18, we show rapidity distribution for a directly produced J/ψ meson (the gg → J/ψg 2 → 2 partonic process). Here, the agreement between the two approaches is not so good. This may be partially understood by the fact that forward J/ψ does not automatically mean also forward gluon. The situation here is similar as for open charm production.
In the present study, we have focused on application of the k t -factorization approach. There are other models in the literature. One of them is the next-to-leading order collinear approach (see, e.g., [23] and references therein). This approach is similar to our approach.
Another alternative approach is the color-evaporation model [24,25]. The old color-evaporation model has been corrected recently [26] and is also able to describe experimental data on inclusive quarkonium production, approximately even some polarization variables [27]. This approach is based, however, on a very different philosophy compared to our approach. We include only color-singlet terms with essentially no free parameters and include all known feed down contributions. As written already in our abstract, we do not see significant room for color-octet terms. So far, we have not seen studies of the improved color-evaporation model for rapidity distributions in forward region. Certainly, such an analysis would be interesting and worthy to do, but it clearly goes beyond the scope of the present approach where we focus on quite a different approach. Compared to the color-evaporation model, there are essentially no free parameters in our approach.

IV. CONCLUSIONS
In the present paper, we have focused on the calculation of cross sections for inclusive prompt production of J/ψ and ψ 0 in forward directions within the k t -factorization approach. In this calculation, NR QCD matrix elements were used with parameters of quarkonia cc wave functions at the origin taken from potential model(s).
In the present calculation, we have used two different sets of unintegrated gluon distribution functions: the Kimber-Martin-Ryskin UGDF based on DGLAP collinear gluon distribution function and the Kutak-Staśto UGDF which includes nonlinear effects at small x values and describes exclusive production of J/ψ [28].
We have included both the direct component and the component related to radiative decays of χ c mesons. In general, they give similar contributions for the integrated cross section.
We have compared our results with the recent results of the ALICE and LHCb Collaborations (small transverse momenta and forward directions) at ffiffi ffi s p ¼ 7 TeV. We have found that using standard KMR UGDF we overestimate the forward production of J/ψ in this case. The biggest contribution is given by radiative decays of χ c mesons. We have proposed how to modify UGDFs to include possible onset of saturation effects. In this mixed UGDF scenario, a reasonable description of the data is possible. We have found that within model uncertainties (UGDFs, renormalization scale, parameters of the nonrelativistic wave function) we can almost describe the production of J/ψ or ψ 0 at low transverse momenta and forward direction including only color-singlet contribution.
We have discussed theoretical uncertainties related to the choice of renormalization scales. In addition, we have discussed some open issues related to the KMR UGDFs. We have shown how to modify the KMR UGDFs to include possible saturation effects. A possible onset of saturation or, in general, nonlinear effects for UGDFs was discussed, especially in the context of the LHCb data. We have found that production of χ c mesons in forward directions is a very good way to search for the onset of saturation, because, as discussed in our paper, it probes smaller values of longitudinal momentum fraction than the J/ψ production and is therefore better suited for that purpose. Such an analysis could be done by the LHCb Collaboration.
We have also made a comparison of our k t -factorization results to the results of the hybrid model where one collinear and one unintegrated gluon are used. This latter approach is used, for example, for forward jet production. We have found that the two approaches (hybrid and k t -factorization) have given very similar results for χ c production (gg → χ c ) and a bit different results for direct gg → J/ψg production. A short explanation is given in the main text.
Other approaches have been briefly discussed, but a comparison of results goes beyond the scope of the present studies.