Central exclusive chi_{c,b} production at high energy colliders and gluon saturation approach

In this work we analyze the central exclusive production of $\chi_{c}$ and $\chi_{b}$ at the LHC and Tevatron energies using the recent unintegrated parton distribution (UGDs) functions disposable by package TMDlib. We provide predictions for LHC using recent UGDs based in CCFM formalism. We also study the underlying uncertainties on this production like the the choice for the unintegrated gluon distribution and factorization scale. Moreover, based on the parton saturation model for the gluon distribution we propose analytical expressions for the rapidity distributions including the prompt production of $J/\psi\gamma$ and $\Upsilon\gamma$.


I. INTRODUCTION
The central exclusive production (CEP) processes are considered as an useful way for testing many aspects of QCD [1]. CEP is a process where the incident hadrons remain intact after the interaction, and an additional simple central system is produced. In Regge language, CEP allows us to study the structure of Pomeron since we have a double Pomeron exchange leading to a specific final state like higgs boson, scalar and tensor mesons including charmonium states as χ c,b mesons. The Pomeron carries quantum numbers of vacuum, so it is a colorless entity in QCD and reactions initiated by Pomerons are experimentally characterized by the "rapidity gap" events. At the Large Hadron Collider (LHC), investigations on CEP are very active [2,3]. Our focus on this study is the production of χ c and χ b within the two-gluon exchange formalism and non-relativistic approach for evaluating the P -wave quarkonium decays [4]. At the LHC, the LHCb Collaboration has done analyses of CEP of χ c mesons, reconstructed in the decay χ c → J/ψγ → µ + µ − [5][6][7] at 7 and 8 TeV. The measured cross sections times branching ratios χ c states reach to dozens of picobarns, which it is in agreement with theoretical predictions. Search for the CEP of χ b mesons has been done, however the background contributions are not completely determined. In any case, most theoretical predictions for the cross section for χ b give values lower than 1 fb which implies very few events. The ALICE Collaboration has recorded zero bias and minimum bias data in proton-proton collisions at a center-of-mass energy of √ s = 7 TeV. Events containing double gap topology have been studied and they are associated to CEP [8,9]. In particular, central meson production was observed and it was verified that K 0 s and ρ 0 are highly suppressed while the f 0 (980) and f 2 (1270) with quantum numbers J P C = (0, 2) ++ are much enhanced (one of us calculated the CEP of these f 0 and f 2 mesons in Ref. [10]) . Such a measurement of those states is evidence that the double gap condition used by ALICE selects events dominated by CEP or related processes. ATLAS and CMS also measured inclusive χ c,b production and cross section ratios for different states are studied [11,12]. Further program on CEP is ongoing in both collaborations [13,14] and also at the Relativistic Heavy Ion Collider (RHIC) [15].
On the theoretical side, applications to χ c production is considered by the Durham group in [16] and references therein, including perturbative QCD and also a non-perturbative component. In Ref. [4] the author calculates the χc and χ b CEP cross sections for the Tevatron, in the forward approximation neglecting the χ c1 and χ c2 states. The Bialas-Landshoff model is applied to χ meson production in Ref. [17], consistent with the χ c0 cross section found by SuperCHIC MonteCarlo for the same parameters [18,19]. The Cracow/Lund group performed calculations [20,21], using a different approach than the Durham group for the unintegrated gluon distributions (UGDs), F g , and taking into account Quasi Multi-Regge Kinematics for the subprocess vertex. The cross sections are found to have a large dependence on the model parameters and the choice of gluon distributions. Interestingly, the cross sections vary by an order of magnitude when using distinct UGDs.
The focus of this work is the central exclusive production of heavy quarkonium (χ c , χ b ) at the Tevatron and at the LHC. One motivation for this study is that χproduction probes the gluon density down to fractional gluon momenta of x ∼ 10 6 , being potentially sensitive to saturation effects. Moreover, this kind of exclusive process is a standard candle in QCD calculations and brings information on the off-forward unintegrated gluon distribution. These objects are poorly constrained in the kinematics investigated here and it is a timely investigation. The formalism of Ref. [4] is considered taking into account the new fitted UGDs available in TMDlib (Transverse Momentum Dependent parton distributions) package [22]. The UGD's used here were those based on CCFM model with three different fitting sets. In addition we consider an analytical UGD based on parton saturation model, i.e. the celebrated GBW saturation model [23]. The purpose to use UDG's based on GBW model was to quantify the deviation of a simple model from a arXiv:2003.09552v1 [hep-ph] 21 Mar 2020 robust model like CCFM and investigate the role played by saturation physics in the UGDs at high energies. The novelty of the results is the updated computation of cross sections using the last CCFM-based UGDs and the predictions for the prompt production of J/ψ + γ and Υ + γ in the very same formalism. We provide analytical expressions for rapidy distributions for prompt χ and V +γ production, Eqs. (8) and (11), based on QCD parton saturation which are quite useful for further phenomenological studies. For the first time the estimation of nuclear effects are predicted for pA and AA collisions within the geometric scaling approach, shown in Eqs. (12). This is crucial for LHC analyses, where the nuclear saturation scale, Q 2 s,A (y) ≈ A 4/9 (10 −5 √ s/m) λ e λy (with λ 0.25) , should be close or larger than the meson mass for a given forward rapidity y.
This paper is organized as follows. In next section the theoretical formalism is presented, including the main building blocks and relevant parameters. In Sec. III, results of the calculations are presented and we compare them with the current literature. In last section we summarize our main conclusions and remarks.

II. THEORETICAL FORMALISM
Here we analyse the central exclusive χ production, p+p(p) → p+χ J +p(p) in the two gluon exchange model [4], where the hard sub-process gg → χ J is initiated by gluon-gluon fusion and the second t-channel gluon (with transverse momentum k ⊥ ) is needed to screen the color flow across the rapidity gap intervals. For the hadronization, a non-relativistic approach is used to compute the P -wave quarkonium decays. Given the forward scattering amplitude, M, the rapidity distribution of χ production will be where y is the rapidity of χ state. Moreover, t i is the momentum transfer squared at the proton (anti-proton) vertices, and B is slope for the proton form factor, which will be taken as B = 4.0 GeV −2 . After integrating on momentum transfer t 1 and t 2 , one obtains [4], where F g are the unintegrated off-forward (skewed) gluon distribution functions, computed at a perturbative scale µ 2 . For the masses and first derivative of radial Pwave functions [24], we use m(χ c0 ) = 3.414 GeV with |R P (0)| 2 c = 0.075 GeV 5 and m(χ b0 ) = 9.859 GeV with |R P (0)| 2 b = 1.42 GeV 5 (notice that the recent values for wave functions 0.1296 and 1.6057 will increase cross section by a factor 1.73 and 1.1.13, respectively) . The rapidity gap survival factor S 2 for central exclusive χ J production can be calculated using the formalism of [16], which gives: S 2 (Tevatron) = 0.058 , S 2 (LHC) = 0.029.
Regarding the UGDs, they can be obtained from the conventional gluon density as [16] where the factor R g takes into account the skewed effects of the off-forward gluon density compared with the conventional gluon density in the region of x x. The factor T 2 [16] will reduce to the conventional Sudakov form factors in the double logarithmic limit. The R g factor used in the literature are R g (Tevatron) = 1.4 , R g (LHC) = 1.2. This factor produces a factor equal to 3.84 (Tevatron) and 2.07 (LHC), since it appears as R 4 g in Eq.
(2). However, we used R g =1 in order to compare predictions to others works in literature.
We will consider here two implementations of UGDs. The first one is the new fitted UGDs available in TMDlib (Transverse Momentum Dependent parton distributions) library [22], based on CCFM model with three different fitting sets. In this case, in the numerical calculations we used the α s (m 2 χ ) to a one-loop order (LO) and 4-flafours (n f = 4). For each CCFM set it was used a specific λ QCD , using prescription given in Ref. [22]. The second UGD we have considered is taken from the saturation model [23], which is analytical and with parameters fitted from DIS data at small-x. It reads as, where Q s (x) = (x 0 /x) λ/2 is the saturation scale giving the transverse momenta transition between the dilute and saturated gluon system. It presents the geometric scaling property, i.e. the UGD depends on the scaling variable k 2 ⊥ /Q 2 s (x) and not separately on x and k ⊥ . In numerical calculation, we use the updated values for the model parameters (fit result including charm): σ 0 = 27.32 mb, λ = 0.248 and x 0 = 4.2 × 10 −5 [25]. Also, at large rapidities we multiply the GBW UGD by the large x threshold, (1 − x) 5 .
In next section, we investigate the uncertainties on theoretical predictions and look closer in the parton saturation model applied to the CEP of quarkonium.

III. RESULTS AND DISCUSSIONS
Here, we will focus on the exclusive production of mesons χ c,b in proton-proton collisions at LHC energies. The present investigation is relevant for the AT-LAS, CMS and ALICE experiments. We will consider the theoretical formalism presented in previous section and investigate its theoretical uncertainties. In particular, we   consider the uncertainty coming from the choice for the unintegrated gluon distribution taking into account different prescription for the renormalization/regularization scale µ 2 . As a cross check we perform predictions also for the lower energy at the Tevatron. The distribution for the meson rapidity is presented and for completeness we compute the corresponding integrated cross sections.

A. Unintegrated gluon distribution
In this section we compare the different set of UGD's for distinct choices for the hard scale. Namely, we investigate the play roled by µ 2 , using the prescriptions m 2 χ /4 ≤ µ 2 ≤ m 2 χ . Let us start by the χ c0 production. In Fig. 1 is shown the behavior on transverse momentum, k 2 ⊥ , for diferent sets of gluon distribution at central rapidity, y = 0, at 14 TeV (LHC energy). At midrapidities the typical gluon momentum fraction is x 1 = x 2 ∼ 10 −4 with a not so hard scale 3 < µ 2 < 11 GeV 2 . In this kinematic range, parton saturation physics (taming the gluon distribution at small-x) could be important. Three sets for CCFM UGD are presented (JH-2013-set1, set A0+ and set B0), as well as the gluon saturation UGD from GBW model and the UGD from GRV94-LO. We can see that the peak occurs for larger k 2 ⊥ in CCFM compared to GRV94 and GBW UGD's. This is directly related to the starting scale Q 2 0 in hard scale evolution and the extrapolation for small gluon transverse momenta. For UGD's extracted from parton saturation physics, the peak occurs around the saturation scale, Q 2 . Therefore, at central rapidity at the LHC the saturation scale is of order Q 2 s 1 GeV 2 , which is confirmed by the numerical results. All results shown in Fig. 1 are computed with µ 2 = m 2 χc . We turn the corresponding analyses for χ b0 production. In this case, at midrapidity gluons have x 1,2 ∼ 10 −3 probed at scale 25 < µ 2 < 100 GeV 2 . In order to single out the uncertainty related to choice of hard scale we compare the set A0+ taking µ 2 = m 2 χ b and µ 2 = m 2 χ b /4.  Fig. 2, using same notation as previous figure. Basically, it is found that for a larger scale the contribution from gluon with large transverse momenta is increasingly important. Once again, the GBW UGD is peaked near the saturation scale and large transverse momenta contributions are exponentially suppressed. In what follows, we investigate the numerical results for the sets we have discussed above.

B. Differential cross section
We calculate the predictions for the rapidity distribution, y, for the exclusive χ c,b production. For sake of completeness we have done a cross check for Tevatron energies, which is shown in Fig. 3 (the curve label is the same as previous figures). Here, µ 2 = m 2 χc0 . The behavior is similar for different sets except for the GBW UGD. The suppression at large rapidities compared to CCFM and GRV94-LO is evident and this trend is more dramatic for LHC energies. The predictions for LHC are presented in Fig. 4, where the choice for distinct sets for UGD's leads to one order of magnitude difference at midrapidities. This can be traced out to the uncertainty on the determination of gluon distribution at very small-x. One has dσ dy (y = 0) ∼ 100 nb for Tevatron and dσ dy (y = 0) ∼ 300 nb with R g = 1 and a sizable spread for LHC case. The evaluations for χ b production is presented in Figs. 5 (LHC) and 6 (Tevatron). In both cases the cross section normalization is strongly dependent on the chosen UGD. Moreover, it is verified that the sensitivity to the hard scale µ 2 is not so strong in the rapidity distribution. This is shown in Fig. 5 for the CCFM set A0+ at LHC energy. One has dσ dy (y = 0) ∼ 100 pb for Tevatron and dσ dy (y = 0) ∼ 500 pb with R g = 1 and a sizable spread for LHC case. It is clear from the present investigation that the main source of uncertainty in the calculations comes from Interestingly, the GBW UGD allows us to obtain an analytical expression for the rapidity distribution. Defining an effective saturation scale, , and using the analytical expression in Eq. (4) one has for their product the following: where the effective saturation scale tends toQ 2 s ≈ Q 2 s (x 2 ) at large backward rapidities whereasQ 2 s ≈ Q 2 s (x 1 ) at large forward rapidities. Moreover, at central rapidity one hasQ 2 s = Q 2 s (x)/2 where x = x 1 = x 2 . In Eq. momentum integration of Eq. (2) in the following form: where, the remaining integral is given by, where A = (3σ 0 /4π 2 α s ) 2 and ξ = m 2 χ /Q 2 s . The integration over τ can be explicitly done, which reads as I s = e ξ (2ξ 2 + ξ 3 )Ei 2 (ξ) − ξ 2 with ξ >> 1 for the values of m c,b . By using the leading terms in the asymptotic expansion of the exponential integral function, Ei 2 (ξ) ≈ e −ξ ξ [1 − (2/ξ) + (6/ξ 2 ) + . . .], an approximate expression for rapidity distribution can be obtained. In the complete case the rapidity distribution is given as, By writing down the expression above in terms of energy and rapidity, one obtains, with an overall normalization given by Here, we consider α s = 0.335 and α s = 0.25 for χ c and χ b , respectively.
As a cross check of evaluation of Eq. (8) (with R g = 1), one obtains dσ theo χ c0 dy (y = 0) = 77 nb for Tevatron, which is consistent with measured value (76 ± 14) nb [26]. Also, LHCb have reported preliminary results on exclusive χ c meson production in the χ c → J/ψ + γ channel [7], in the rapidity kinematic region 2.0 < η < 4.5. The cross section times branching ratios (taken from PDG [27]) for production in the LHCb acceptance (ε s = 0.76) given by saturation model for χ c0 is 29 pb with large uncertainty compared to the measured value 9.3 ± 4.5 pb. Notice that the χ(J = 1, 2) production amplitudes vanish in the perturbative two-gluon exchange model we are using. However, by considering the normalization of gg → χ J and the mass difference, we estimate that the cross sections for those states could be a factor ∼ 0.7 and 0.06 times the cross section fo J = 0 state. This gives 20.3 pb and 1.74 pb, compared to experimental values 16.4 ± 7.1 pb and 28 ± 12.3 pb, respectively. For comparison, the corresponding prediction from SuperCHIC [18] is 14 pb, 9.8 pb and 3.3 pb, respectively.
Finally, we consider the predictions for quarkonium CEP cross sections at different collider energies. In Table I we show the differential cross section for the central exclusive χ c (and χ b0 ) process at RHIC, Tevatron and LHC energies. It was verified that the predictions are a factor 2 higher than those from Durham model for χ c0 [16].
The perturbative two-gluon exchange model can also be used to compute the prompt production of V = J/ψ, Υ in the process p + p(p) → p + V γ + p(p). The CEP cross section for this channel is given by [4], (10) where y γ and y V are the photon and meson rapidities. The meson transverse momentum is denoted by p ⊥ with a transverse mass m ⊥ = m 2 V + p 2 ⊥ . Now, For the masses and radial S-wave functions at origin [24], we use m(ψ) = 3.096 GeV with |R S (0)| 2 c = 0.81 GeV 5 and m(Υ) = 9.46 GeV with |R S (0)| 2 b = 6.48 GeV 5 . Once again, the saturation model gives an analytical solution for the integral I V . Therefore, the differential cross section is written as, The numerical calculation for the differential cross sections for production of J/ψ + γ and Υ + γ at central rapidity are presented in Table I (integrated over photon rapidity and meson transverse momentum). We see that the cross sections for S-wave quarkonia are comparable or larger than those for P-wave states times Br(χ → V γ) ∼ 10 −2 at least at y = 0. This is in disagreement with conclusions presented in Ref. [4], which predict that the leading contributions to CEP of S-wave quarkonium are the feeddown contributions from P-wave decays. In Figs. 7 and 8 we present the differential cross sections, Eq. (11) integrated over photon rapidities, in terms of meson transverse momentum at y V = 0 for the differents UGDs discussed before. We found that the main contribution for the meson p ⊥ -spectra comes from the region p ⊥ < ∼ m V . Before discussing the integrated cross sections for different models of UGDs, let us discuss the extrapolation of the saturation model to nuclear collisions. It is found in Eqs. (8) and (11) that the rapidity distributions take the form dσ/dy ∝ (σ 0 ) 2 Let us consider the label 1 for projectile and 2 for the target and take into account the geometric scaling property in nuclear reactions demonstrated in Ref. [28]. Namely, for the unintegrated gluon distribution in a nucleus we could replace in Eq.
. Here, δ = 0.79 and Q s (x) is the saturation scale for the proton case. Explicitly for pA collisions we get: and for AA collisions similarly we obtain, The crude approximation above based on geometric scaling property can be compared to the sophisticated calculations using SuperCHIC 3 Monte Carlo [19], which implements CEP in nuclear collisions.

C. Integrated cross section
Based on the rapidity distribution obtained above, the integrated cross section can be computed. Results for Tevatron and LHC energies are shown in Tab. II. The output for the different UGS's sets are presented for χ c0 (χ b0 ) in units of nanobarns and disregarding the skewedness effects. At the LHC, the larger cross section corresponds to the GRV94-LO UGD, whereas GBW UGD gives lowest values. Based on the theoretical uncertainty from UDG's one obtains σ(χ c0 ) = 3619 ± 3241 nb and σ(χ b0 ) = 4.7 ± 3.5 nb at LHC with R g = 1.
In case of χ c,b to be measured by detecting their radioactive decays to quarkonium plus photon the final cross section for quarkonium production from χ 0 feeddown decays would be: for χ c0 exclusive production. On the other hand, for the χ b0 production one obtains, where it has been assumed for simplicity that the production cross section for χ b in the 2P and 1P states are of same order of magnitude. We now compare our results to other models available in literature. The SuperCHIC MC generator implements the Durham model and the χ c cross sections according to it at √ s = 7 TeV, over the full kinematic range and including the branching ratios of χ c → J/ψγ → µ + µ − are 194 pb, 133 pb and 44 pb for χ c0 , χ c1 and χ c2 respectively. The predictions in this work are considerably larger than SuperCHIC [18], with the saturation model being the closest one (≈ 300 pb). The measured value by LHCb is 134 pb. Interestingly, high cross sections were also obtained in Ref. [29], using Bialas-Landshoff (BL) formalism implemented in DPEMC Monte carlo. The BL model was also applied to χ production in Ref. [17], with a cross section of 350 nb for χ c0 production at the LHC. Predictions are also consistent in order of magnitude with results presented by Cracow/Lund group in Ref. [20], including the large uncertainty from the choice for the UGDs. The same occurs for results from 3-Pomeron model [30], which predicts σ(χ c0 ) = 212 ± 53 nb at 7 TeV (future version of ExDiff Monte Carlo [31], based on theoretical framework of Ref. [30] will include quarkonium production).
The predicted χ b0 cross section is much higher than Durham group, probably due to a different coupling of two gluons to the χ b . The non-perturbative two-gluon model (BL) from Ref. [17] predicted a total χ b0 cross section of 0.3 nb at √ s = 14 TeV, which it is consistent with present calculations using CCFM-setB0 and CCFM-setB0+ or CCFM-setA0 . Moreover, in Ref. [4] was predicted a total cross section of 0.88 nb at the Tevatron, in agreement in order of magnitude with present work. A Regge-eikonal approach for CEP is investigated in Ref. [32], which predicts σ(χ b0 ) 0.16 nb and 1.3 nb at Tevatron and LHC, respectively. Once again, results presented in Table II are consistent with those calculation.
As final comment, besides being considered the theoretical uncertainties on UGDs and hard scales, other quantities are source of additional uncertainty as the slope of the proton form factor, B, the gap survival factor and value of the wave-function at the origin.

IV. SUMMARY
We have investigated the central exclusive production of χ c,b0 in hadron-hadron collisions. In the theoretical calculations, we take into account the perturbative twogluon model and non-relativistic approximation for meson wave functions. The numerical results are obtained for different models for the unintegrated gluon distribution, including a analytical parametrization from parton saturation approach. It was found that the main uncertainty in the prediction comes from the choice for the UGD. We verified that the different prescription for the hard scale µ 2 has a small effect for χ b production. We have also shown the saturation model for the UGD allow us to obtain an analytical expression for the rapidity distribution both for χ J and prompt production V γ. It depends explicitly on the effective saturation scale, Q s (x 1 , x 2 ), and can be easily extended to pA or AA collisions using arguments of geometric scaling. That is, the nuclear saturation scale is rescaled compared to the nucleon one, Q 2 s,A ∝ A 1/3 Q 2 s,p . We found that the corresponding scaling is σ pA coh ∼ A α pA σ pp for proton-nucleus (with α pA = (4δ − 2)/3δ ≈ 4/9) and σ AA coh ∼ A α AA σ pp (with α AA = (8δ + 2)/3δ) ≈ 10/3) for nucleus-nucleus collisions, respectively.
In summary, our study demonstrated that the CEP of mesons is a powerful tool to investigate the perturbative QCD dynamics and in proton-proton collisions at the LHC. This shall stimulate further experimental and theoretical studies.