Central exclusive production of scalar and pseudoscalar charmonia in the light-front $k_T$-factorization approach

We study exclusive production of scalar $\chi_{c0}\equiv \chi_c(0^{++})$ and pseudoscalar $\eta_c$ charmonia states in proton-proton collisions at the LHC energies. The amplitudes for $gg \to \chi_{c0}$ as well as for $gg \to \eta_c$ mechanisms are derived in the $k_{T}$-factorization approach. The $p p \to p p \eta_c$ reaction is discussed for the first time. We have calculated rapidity, transverse momentum distributions as well as such correlation observables as the distribution in relative azimuthal angle and $(t_1,t_2)$ distributions. The latter two observables are very different for $\chi_{c0}$ and $\eta_c$ cases. In contrast to the inclusive production of these mesons considered very recently in the literature, in the exclusive case the cross section for $\eta_c$ is much lower than that for $\chi_{c0}$ which is due to a special interplay of the corresponding vertices and off-diagonal UGDFs used to calculate the cross sections. We present the numerical results for the key observables in the framework of potential models for the light-front quarkonia wave functions. We also discuss how different are the absorptive corrections for both considered cases.


I. INTRODUCTION
The central exclusive diffractive processes in proton-proton collisions at high energies have attracted recently a lot of attention. These processes lead to very unusual final states. For example, in the central exclusive production one produces one or a few particles at central rapidities which are fully measured. There are no other tracks in the detectors. The incoming protons remain intact (in the virtue of "elastic diffraction") or are excited into small mass hadronic systems, which disappear into the beam pipe. We consider here simultaneously two such reactions, pp → p χ c0 p and pp → p η c p, which are well suited to be analysed in the framework of the so-called Durham model formulated by Khoze, Martin and Ryskin (see Ref. [1] and references therein). From the experimental point of view, there is a rapidity gap, between each of the protons and the produced χ c0 or η c states. These processes hence provide a very clean environment for the study of the produced hadronic systems tightly connected to poorly known soft and semi-hard QCD dynamics. For a review of conceptual and experimental challenges with such central exclusive production (CEP) reactions, see for example Ref. [2].
The theory of the CEP of single χ cJ , J = 0, 1, 2 mesons, with a correct account for the spin of the mesons and precise kinematics of the production process has been worked out earlier by Pasechnik, Szczurek and Teryaev (PST) in a series of papers [3][4][5]. The numerical calculations were done for the Tevatron energies. In this analysis, the non-relativistic QCD (NRQCD) methods were applied. So far, only CEP of light pseudoscalar mesons was discussed in the literature [3,6]. There rather nonperturbative effects strongly dominate (see Ref. [6]). Very recently in Ref. [7] the production of χ c0 at the LHC was discussed in the k Tfactorisation and saturation dipole-model inspired approaches. The analysis was performed there in the NRQCD approach and using a single model for the unintegrated gluon distribution (UGDs) and a particular prescription for the off-diagonal UGD. Given a particular importance of the CEP of heavy quarkonia for ongoing and future experimental studies, we revisit and extend this analysis to account for additional effects and sources for theoretical uncertainties (such as the shapes of the charmonia wave functions and a treatment of the absorptive corrections, as well as an accurate treatment of the phase space and production kinematics) as well as incorporate the pseudoscalar η c final state for the first time.
Recently our group showed how to include relativistic corrections for the inclusive production of η c [8] and very recently for inclusive production of χ c0 [9] using the light-cone wave functions of the charmonia derived from the well-known cc interquark potential models. It is the aim of the present paper to do a similar study for the exclusive case. In addition, there is no such a study on the pp → ppη c CEP available in the literature. In contrast, η c (1S) was measured by the LHCb collaboration in the inclusive case for √ s = 7, 8, 13 TeV [10].
Is such a measurement possible for the exclusive production of η c ? This study is a first step to address this important question. An analysis of inclusive diffractive production of η c (1S) was done recently in [11].
In the present paper, we wish to discuss in parallel the exclusive production process of both the scalar χ c0 and pseudoscalar η c quarkonia For illustration of the corresponding production mechanism initially proposed by the Durham group [1], see Fig. 1. We start with a brief introduction into the formalism for pp → p χ c0 p and pp → p η c p reactions based upon the Durham model of CEP [1] setting up the necessary notation and conventions. In this model, a quarkonium state is produced via fusion of two virtual active gluons accompanied by an extra exchange with a screening gluon as is shown in Fig. 1. The additional exchange of the gluon provides colour conservation and hence the effective color singlet exchange in the t-channel. As a result, in the final state, in addition to the meson produced mainly at central rapidities, there are two forward protons that retain most of their initial energy. We intend to calculate the integrated cross sections for such processes as well as several differential distributions relevant for future measurements. We wish to discuss both the hard and soft processes involved in these reactions in the light-front QCD approach, to consider several prescriptions on how to calculate the off-diagonal UGDs (some of them have already been used previously in the literature) and to estimate the absorptive corrections in the differential distributions.
A. The light-cone amplitude for g * g * → χ c0 (1P ) process The gluon-gluon fusion vertex is proportional to the reduced amplitude T µν as follows: where α s is the strong coupling, N c = 3 and t a are the number of colors and SU(3) group generators in QCD, respectively, and while the relevant kinematical variables are displayed in Fig. 1. Here, we have the projector on transverse polarization states with X = (q 1 · q 2 ) 2 − q 2 1 q 2 2 . The longitudinal polarization vectors read as follows The convoluted form of reduced amplitude can be written as in terms of the light cone vectors n ± ν = (1, 0, 0, ±1). The form factors here G i (q 2 1 , q 2 2 ) have the integral representations in terms of the P -wave charmonia wave function ψ χ (z, k) (see Ref. [9] for more details) where z is a c-quark (orc-antiquark) momentum fraction, k is the relative cc transverse momentum, m c is the mass of c-quark, and the shorthand notations have been introduced.
B. The light cone amplitude for g * g * → η c process In Ref. [8], we introduced the covariant form of the vertex for two off-shell gluon fusion into η c meson: where The convoluted form reads: and (φ 1 − φ 2 ) is azimuthal angle between q 1 and q 2 . We then express I(q 2 1 , q 2 2 ) in terms of light-cone wave functions as follows [12] I(q 2 1 , q 2 2 ) = 4m c where ψ η (z, k) is the wave function of η c (1S) meson.

III. MATRIX ELEMENT FOR pp → ppM REACTION
The amplitude for the CEP process for a given meson V ≡ χ c0 , η c reads 1 : in terms of the "active" (fusing into V ) x 1,2 and "screening" x ′ (connecting both proton lines) gluon momentum fractions. The screening gluon carries a transverse momentum Q, while the transverse momenta of active gluons are denoted by q 1 , q 2 . The generalized unintegrated gluon distributions (UGDs) also depend on the hard scale of the process µ (see below). The 2 → 3 total cross section can be calculated generically as follows: or, following a simplification done in Ref. [13], as σ = 1 2s is the relative azimuthal angle between the outgoing protons, s is the pp center-of-mass energy squared, y is rapidity of the outgoing meson V .

IV. DIFFERENT APPROACHES TO OFF-DIAGONAL GLUON DENSITIES
In the forward limit of small t 1,2 → 0 corresponding to Q 2 ≃ q 2 1,2 ≡ Q 2 ⊥ , the generalized UGDs in Eq. (3.1) are simplified and are considered as functions of only one transverse momentum, i.e.
The Khoze-Martin-Ryskin (KMR) prescription for the off-diagonal ("skewed") UGD includes a Sudakov form factor T g (q 2 ⊥ , µ 2 ) and is typically written as [1] with gluon virtualities q 2 ⊥ ≡ q 2 etc playing a role of the momentum scale squared in the collinear gluon density xg(x, q 2 ⊥ ), and with the nucleon form factor F (t) often parameterised in the following two ways with m p being the proton mass, corresponding to the isoscalar nucleon form factor [14] or the QCD elastic profile factor, respectively. The Sudakov form factor is taken as: with the hard scale µ 2 = M 2 V and ∆ = k ⊥ /(k ⊥ + M V ). Regarding the longitudinal momentum fractions, central diffractive production is dominated by the region x ′ ≪ x 1,2 ≪ 1. We thus compute the skewedness correction R g in Eq. (4.2) using a method proposed and derived for the collinear off-diagonal gluon distributions [15]: In a slightly off-forward case t 1,2 = 0, the choice of Q ⊥ in the off-diagonal KMR gluon in Eq. (4.2) becomes somewhat arbitrary. In practical calculations, we use the so-called "minimum prescription" proposed by the Durham group, by substituting Q 2 ⊥ → min(Q 2 ⊥ , q 2 ⊥ ) in Eq. (4.2), with transverse momentum of an active gluon q ⊥ and transverse momentum of the screening gluon Q ⊥ . In addition, we suggest a geometrical average of active and screening gluon momenta as Q 2 ⊥ → Q 2 ⊥ q 2 ⊥ -an option, called BPSS in the following, for brevity. We vary our results by also using the modified off-diagonal CDHI gluon defined as [16] whereQ 2 = (Q 2 ⊥ + q 2 ⊥ )/2. In order to take into account the saturation effects, we use of the simplest saturation-based UGD inspired by the Golec-Biernat-Wüsthoff (GBW) model [17]. In order to extrapolate it into the off-diagonal domain, we use the prescription proposed in Ref. [3] (further referred to as the PST prescription): where f GBW g is the diagonal GBW UGD, and In practical calculations, we have used the following fitted values of the GBW parameters: σ 0 = 29.12 mb, λ = 0.277, x 0 /10 −4 = 0.41, with fixed α s = 0.2.

V. NUMERICAL RESULTS
In the CEP processes at high energies, it is mandatory to consider gluons carrying very small longitudinal momentum fractions x. For this purpose, in practical calculations we use a few parton distribution functions (PDFs) introduced by the Dortmund group: JR14NLO [18] (Q 2 0T = 0.8 GeV 2 ), GJR08NLO [19] (Q 2 0T = 0.5 GeV 2 ) and GRV94NLO [20] (Q 2 0T = 0.4 GeV 2 ). In Fig. 2, we illustrate the shape of the corresponding gluon PDFs at scales and longitudinal momenta fractions typical for the considered pp → ppη c and pp → ppχ c,0 CEP processes. In the range of scales under discussion, the gluon PDFs from the literature differ considerably. We do not employ the Durham or CTEQ PDFs for which the initial scales for evolution are rather high making them difficult to be applied in the context of the exclusive reactions discussed here. The total cross sections computed over the full phase space for each PDF mentioned above are listed in Tables 1 and 2, for χ c0 and η c , respectively. The integrated cross section for pp → ppχ c0 at √ s = 13 TeV is shown in Table 1 for different off-diagonal UGD prescriptions for the effective Q 2 iT summarized as follows: ) .
These prescriptions lead to similar cross sections of the order of 1 µb before including absorption effects. The corresponding gap survival factor is of the order of 0.1 as will be discussed at the end of this section. In Table 2 we present similar results for η c production. The total cross section for the η c production is 3-4 orders of magnitude smaller than that for χ c0 , i.e. surprisingly small. The cross section for the PST prescription for off-diagonal gluon is quite similar as for the In order to obtain the cross section, several gluon distributions were used with Q 2 0T ≥ 0.4 GeV 2 for GRV94NLO, Q 2 0T ≥ 0.5 GeV 2 for GJR08NLO, and Q 2 0T ≥ 0.8 GeV 2 for JR14NLO. The light-cone form factor for the gg → χ c0 coupling was calculated using the Buchmüller-Tye potential (for more details, see Ref. [9]) No gap survival factor is included here.
Durham and CDHI prescriptions in the case of χ c0 , while the spread in the total cross section for η c is much higher. In Fig. 3 we show rapidity distribution of χ c,0 (left) and η c (right) quarkonia CEP. We show results for the Durham (min) and CDHI prescription for the off-diagonal UGDFs. We present results for R g = 1 as well as with R g calculated according to the Shuvaev prescription  Table 1 but for η c meson. The light-cone form factor for the gg → η c (1S) coupling was calculated using the power-law potential (for more details, see Ref. [8]) (see Eq. (4.5)). Inclusion of R g increases the cross section by a factor of 3-4. While for χ c0 the difference of the results for the Durham prescription and the CDHI prescription is small, for η c the difference is of the order of magnitude size.
The distribution in transverse momentum are shown in Fig. 4. The distribution for η c and χ c0 CEP are somewhat different. The maximum of the cross-section for η c is at p T ∼ 1 GeV and the dip at vanishing p T is more pronounced.
In Fig. 5 we show two-dimensional distributions in (t 1 , t 2 ) (t 1 , t 2 are four-momenta squared transferred in the proton lines), for pp → ppη c (1S). In the left panel we show the result  for the Q 2 it = q 2 it Q 2 t prescription proposed above for effective transverse momenta used in the KMR method as well as the PST prescription for off-diagonal UGD using the diagonal GBW UGDF. The off-diagonal KMR UGD with the newly proposed prescription leads to some numerical irregularities in the (t 1 , t 2 ) space and extremely small cross sections. In contrast, the distribution for the PST prescription with the GBW UGDF is regular and the corresponding cross section is much larger (see also Table 2).
In Fig. 6 we show similar results for two other prescriptions for the off-diagonal KMR UGDFs: with the Durham prescription -left panel, and the CDHI prescription -right panel. The two prescriptions do not lead to numerical issues but generate extremely small cross sections.
For completeness in Fig. 7 we show similar results for the χ c0 production. In contrast to the η c CEP process, in this case both prescriptions for effective transverse momenta (Durham and CDHI prescriptions) lead to fairly similar results. Here the cross sections are peaked at t 1 = 0, t 2 = 0.
In Fig. 8 we show relative azimuthal angle (between outgoing protons) distributions. The distribution for χ c0 (left) is very different than that for η c (right). While for χ c0 there is one maximum for the back-to-back configurations, there are two maxima for η c . The cross section vanishes in the back-to-back kinematics in the case of η c CEP. The exact position of the maxima depends on the details of the treatment of the off-diagonal UGDs so their experimental identification could pin down the correct theoretical modelling of these objects.
Finally, we wish to compare our results for the exclusive reactions pp → ppη c and pp → ppχ c0 with their inclusive production counterparts as calculated recently in Refs. [8] and [9]. In Figs. 9 and 10 we show the numerical results (rapidity and transverse momentum distributions) for η c and χ c0 , respectively. While for η c production the cross section for the exclusive process is a few orders of magnitude lower than that for the inclusive case, this is quite different for χ c0 meson. Both for rapidity and transverse momentum distributions the results for the exclusive case are very different compared to the inclusive case.

VI. ABSORPTIVE CORRECTIONS
It is understood that the Born-level cross sections receive absorptive corrections through hadronic rescatterings at large distances. These are related to the interactions of spectator partons [21]. They give rise to the so-called gap survival probability in exclusive reactions. The calculation of the latter poses a difficult problem, which has not been solved yet in a way fully consistent with the perturbative QCD approach to the production amplitude.
Numerous approaches exist in the literature, some of them are based on soft multi-Pomeron exchanges [22][23][24], while other approaches avoid the decomposition into Born term and absorptive correction altogether treating the absorptive effects dynamically [25,26] and at the amplitude level in the dipole picture [27][28][29], and some of them relate the gap survival probability to the absence of multiparton interactions [30,31].
It is also understood that the gap survival must depend on the kinematics of the process. Here, we wish to discuss the absorptive corrections at the amplitude level, in a simple quantum-mechanical treatment. To this end, we adopt a simple effective Reggeon Field Theory motivated approach.
In the simplified case where only "elastic rescattering" is taken into account, the amplitude looks as follows: Here Y = log(s/m 2 p ) is the rapidity difference between the colliding beams at center-ofmass energy √ s, y is the cm-rapidity of the produced meson V , and p 1,2 are the transverse momenta of outgoing protons. The absorptive correction is then computed as follows At √ s = 13 TeV we take σ pp tot = (110.6 ± 3.4) mb and the nuclear slope B el = (20.36 ± 0.19) GeV −2 [32]. In a double-Regge approach, the Born-level amplitude has the form Here, R IP (y, p 2 ) are the Pomeron Regge-propagators, and V (p 1 , p 2 ) is the IPIP → Meson vertex.
Let us now briefly discuss the vertices V (p 1 , p 2 ). The most general form of the Pomeron-Pomeron-particle vertex for a spinless particle can be written as a Fourier expansion: For a scalar particle, all V − n = 0, while for the pseudoscalar V 0 = 0, V + n = 0. For definiteness, let us concentrate on only the first order, n = 1. We thus adopt (V 0+ for scalar and V 0− for pseudoscalar state): We further neglect a possible dependence of vertices V ± i on p 2 1 and p 2 2 .      Our amplitude is normalized in such a way that the expression dσ = 1 512π 4 s 2 |M(Y, y, p 1 , p 2 )| 2 dyd 2 p 1 d 2 p 2 d 2 p δ (2) (p + p 1 + p 2 ) (6.7) holds. We will now concentrate on central diffractive production, i.e. we fix the meson rapidity to be y = 0. Below we adopt √ s = 13 TeV. We can therefore forget about the Regge-propagators in Eq. (6.4), and without loss of generality we write Then, using the vertices of Eq. (6.6), the transverse momentum distributions of the mesons at the Born level are obtained as Now, the absorptive corrections require the evaluation of the loop integral δM(Y, 0, p 1 , It is useful to introduce the dimensionless quantities Then, the absorptive corrections are obtained as for the scalar and pseudoscalar meson, respectively. We now adjust the constants V 0 , V ± 1 , as well as B D , to our numerical results obtained for the Born-level amplitude.
In Tables 3 and 4, we show the parameters obtained for different prescriptions for the generalized unintegrated gluon distribution, GBW as well as the CDHI and Durham prescriptions for the GJR08NLO gluon distribution. We also show the gap survival factors (6.14) We observe that depending on the gluon distribution used, we obtain for the χ c the gap survival values of S 2 = 0.2 ÷ 0.32, while for the η c production rather similar values S 2 = 0.26 ÷ 0.3 emerge. This similarity may be rather surprising, as the η c amplitude, due to the vanishing in forward direction, is more peripheral than the one for χ c production. However notice that the both reactions have significantly different values of the effective diffraction slopes B D . Our simplified double-Regge approach works reasonably well. In Figs. 11, 12 and 13 we show the cross section dσ/dydp T at y = 0 for the χ c for the three different generalized UGD prescriptions. Shown is the exact numerical result of the Born amplitude (solid line) as well as the result of our effective Regge amplitude fit (long-dashed line). By the short-dashed line we show the differential cross section including absorptive corrections on top of the Regge amplitude Born term. We see from these figures that the effective Regge amplitude form is reasonably accurate for p T ∼ < 1.5 GeV, with a slight ambiguity in the slope B D . In the case of the η c shown in Fig. 14, the effective Regge fit works almost perfectly for the PST-GBW and CDHI prescriptions, while for the Durham case the description is rather poor. Note however that for the η c case the PST-GBW gives by far the largest cross section.

VII. CONCLUSION
In the present paper we have calculated the key observables of central exclusive χ c0 and η c quarkonia production in proton-proton collisions at the LHC within a formalism proposed earlier ago by the Durham group for central exclusive Higgs boson production.
The χ c0 meson CEP was already computed in the literature previously, while η c production has been analysed here for the first time. Compared to the previous calculations we have used modern versions of collinear gluon distributions to generate off-diagonal unintegrated gluon distributions.
In the present analysis we have also used the gg → η c and gg → χ c0 transition amplitudes calculated using the light-cone cc wave functions obtained in the framework of potential models. We have performed similar calculations for inclusive production of η c and χ c0 very recently and showed that one can very well describe the experimental data for η c (1S) meson measured in last few years by the LHCb collaboration. Our previous results showed that in the inclusive case the cross section for η c is significantly larger than that for χ c0 .
It was the main aim of the present paper to make a similar analysis for the exclusive production case, not done so far in the literature. In contrast to the inclusive case, we have found that for the CEP the situation reverses i.e. the corresponding cross section for exclusive η c production is considerably smaller than its counterpart for exclusive χ c0 production, at least, for the hard part obtained using the Durham or Cudell et al prescriptions for calculation of the scale in the off-diagonal unintegrated gluon distribution. The reason is a specific interplay of the off-diagonal UGDs and virtual gluon -virtual gluon -quarkonium vertex.
We also proposed a way to calculate the soft effects (in the region of small gluon transverse momenta) using the GBW UGDF and a simple (PST) prescription for its off-diagonal extrapolation. In this case, the cross section is only slightly smaller for η c than for χ c0 production. We have also discussed to which extent the absorption effects for pp → ppη c are different than those for pp → ppχ c0 . We find that the absorptive corrections are rather similar in spite of very different (t 1 , t 2 ) dependence for the corresponding Born amplitudes. However, the very different slopes of the Born-amplitudes compensate this effect to some extent.
It would be desirable to measure the cross section for pp → ppη c by identifying η c e.g. in the pp decay channel as was done in the inclusive case. It could be interesting to estimate the signal-to-background ratio before the real experiment. The pp → pppp continuum was calculated previously by Lebiedowicz, Nachtmann and Szczurek [33] and a first experimental evidence was obtained very recently by the STAR collaboration at RHIC [34]. Also the pp → ppγγ reaction could be considered as an alternative to measure the pp → ppη c reaction.