Relativistic corrections to \eta_c-pair production in high energy proton-proton collisions

On the basis of perturbative QCD and the relativistic quark model we calculate relativistic corrections to the double $\eta_c$ meson production in proton-proton interactions at LHC energies. Relativistic terms in the production amplitude connected with the relative motion of heavy quarks and the transformation law of the bound state wave functions to the reference frame of moving charmonia are taken into account. For the gluon and quark propagators entering the amplitude we use a truncated expansion in relative quark momenta up to the second order. Relativistic corrections to the quark bound state wave functions are considered by means of the Breit-like potential. It turns out that the examined effects decrease total nonrelativistic cross section more than two times and on 20 percents in the rapidity region of LHCb detector.

The high energy and luminosity of the LHC makes it possible to observe the pair charmonium production processes and to measure the corresponding cross sections with sufficiently high accuracy. The result of the measurement of the pair J/ψ meson production by the LHCb Collaboration was published in [1] and discussed many times on different workshops [2]. This process together with the charmonium and associated open charm production can be considered as a probe of the quarkonium production mechanism. According to nonrelativistic QCD (NRQCD) [3] and collinear parton model [4][5][6] the predictions of observed cross section in the leading order in α s can be obtained by the use of parton distribution functions and a set of local non-perturbative charmonium production color singlet and color octet matrix elements [7][8][9][10][11]. In the proton-proton collisions, additional contributions from other mechanisms, such as the double parton scattering (DPS) or the intrinsic charm content of the proton to the total cross section are possible [12][13][14]. The processes of quarkonium production in proton-proton interaction are generally described using the scale-dependent parton density functions. They are calculated as functions of the Bjorken variable x at some factorization scale within the approach of the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations. However, the double charmonium production in pp collisions at high energies can be sensitive to the details of the parton kinematics. Therefore, it is more appropriate to use the parton distributions unintegrated over the transverse momentum k t in the framework of the k t -factorization [12]. There exists another source of theoretical uncertainty related with the pair charmonium production which gives essential modification of the cross sections. It is connected with the account of relative motion of heavy quarks forming the bound states. As was shown in [15] these relativistic corrections significantly change the cross section of the pair charmonium production in pp interaction obtained in non-relativistic approximation. The detailed investigation of this relativistic mechanism for exclusive double charmonium production in e + e − annihilation [16][17][18] evidently shows that it is impossible to obtain the reliable theoretical predictions for observed quantities without an account of relativistic corrections. Finally, the next to leading order QCD corrections to the production amplitudes also should be taken into account.
The strategy of experimental investigations can be directed on the study of such physical reactions in pp collisions in which one of the described mechanisms of quarkonium production is dominant. Unfortunately, as we know at present time all enumerated mechanisms have important effect in the pair charmonium production and their contributions to the total cross section should be taken into account to obtain high accuracy theoretical result.
In this work we continue the study of relativistic effects in the inclusive pair charmonium production by considering the process p+p → η c +η c +X. Our calculation of the production cross section is performed on the basis of relativistic quark model used previously for the investigation of relativistic corrections to the other reaction p+p → J/ψ+J/ψ+X in [15]. We work within the single-parton scattering (SPS) mechanism in which the basic contribution to the charmonium production is determined by the gluon-gluon fusion. The aim of the present study consists also in the analysis of some uncertainties regarding the choice of parton distribution functions (PDF). In spite of existing difficulties in the detecting of η c meson pairs it is thought that in new run of the LHC this process will be studied more successfully.
The differential cross section dσ for the inclusive double charmonium production in proton-proton interaction can be presented in the form of the convolution of partonic cross section dσ[gg → η c η c ] with the parton distribution functions of the initial protons [4][5][6][7]: where f g/p (x, µ) is the partonic distribution function for the gluon in the proton, x 1,2 are the longitudinal momentum fractions of gluons. The cross section formula (1) contains the factorization of the long distance PDFs and the short distance gluon fusion cross section dσ[gg → η c η c ] with the factorization scale µ. Neglecting the proton mass and taking the c.m. reference frame of initial protons with the beam along the z-axis we can present the gluon on mass-shell momenta k 1,2 = x 1,2 √ S 2 (1, 0, 0, ±1). √ S is the center-of-mass energy in proton-proton collision.
In the quasipotential approach the double charmonium production amplitude for the basic parton subprocess g + g → η c + η c can be expressed as a convolution of a perturbative production amplitude of two c-quark andc-antiquark pairs T (p 1 , p 2 ; q 1 , q 2 ) and the quasipotential wave functions of the final mesons Ψ ηc [17,18]: where p 1 and p 2 are four-momenta of c-quark andc-antiquark in the pair forming the first η c particle, and q 2 and q 1 are the appropriate four-momenta for quark and antiquark in the second meson η c . They are defined in subsequent transformations in terms of total momenta P (Q) and relative momenta p(q) as follows: This expression describes the symmetrical escape of the c-quark andc-antiquark from the mass shell. In Eq. (2) we integrate over the relative three-momenta of quarks and antiquarks in the final state. The systematic account of all terms depending on the relative quark momenta p and q in (1) is important for increasing the accuracy of the calculation. p = L P (0, p) and q = L Q (0, q) are the relative four-momenta obtained by the Lorentz transformation of four-vectors (0, p) and (0, q) to the reference frames moving with the four-momenta P and Q.
The relativistic wave functions of the bound quarks, accounting for the transformation from the rest frame to the moving one with four momenta P and Q, are the following [17,18]: where the hat symbol means a contraction of the four-vector with the Dirac gamma matrices; v 1 = P/M, v 2 = Q/M; ǫ(p) = m 2 + p 2 , m is c-quark mass, and M is η c charmonium mass, M = 2m. The amplitude (2) is projected onto a color singlet state by replacing v i (0)ū k (0) with a projection operator of the form v i (0)ū k (0) = γ 5 (1 + γ 0 )δ ik /2 √ 6. The relativistic wave functions in Eq. (4) are equal to the product of the wave functions in the rest frame Ψ ηc 0 and the spin projection operators that are accurate at all orders in |p|/m [17,18]. Our derivation of relation (4) accounts for the transformation law of the bound state wave functions from the rest frame to the moving one with four momenta P and Q, which was obtained in [19,20]. The physical interpretation of the double charmonium production amplitude is the following: we have a complicated transition of two heavy c-quark andc-antiquark produced in gg-fusion outside the mass shell and their subsequent evolution firstly on the mass shell (free Dirac bispinors) and then to the quark bound states. In the spin projectors (4) we have p 2 , q 2 = M 2 /4 − m 2 so that we can consider these structures as a transition form factors for the heavy quarks from the free states to the bound states.
In the leading order in the strong coupling constant α s , there are 39 Feynman diagrams contributing to the pair production of η c mesons. They can be divided into two sets shown in Figs. 1 and 2, respectively. Their total contribution to the production amplitude (2) can be presented in the following form: where we explicitly extracted the relativistic normalization factors √ 2M of quasipotential wave functions. The construction and transformation of the production amplitudes is performed by means of the package FeynArts [21] for the system Mathematica and Form [22].
The integrand term in (5) containing the trace of the amplitude M represents the contribution of 31 diagrams in Fig. 1 and equals up to the wave functions definitions (4) the analogous expression in the case of pair J/ψ production, which can be found in Ref. [15]. The second integrand term in (5), coming from additional 8 diagrams in Fig. 2 for η c production amplitude, has the form where the Mandelstam variables for the gluonic subprocess gg → η c η c are: φ is the angle between P and the z-axis. The transverse momentum P T of η c and its energy P 0 can be written as In order to calculate relativistic corrections contained in the production amplitude (6) we expand the inverse denominators of gluon and quark propagators as series in relative quark momenta p and q: There are 16 different propagators in the amplitude (5), which have to be expanded in the manner of Eqs. (11). Then, preserving in the expanded amplitude terms up to the second order in relative quark momenta p and q, we can perform angular integration using the following relations for S-wave charmonium: where R S (p) is the radial charmonium wave function.
As a result of the described transformations, we obtain the following general structure of pair η c production amplitude (5): where A i are the functions of variables s and t. Due to the bulkiness of corresponding expressions for coefficient functions A i we do not present them here in exact form.
In order to find the differential cross section for the gluonic subprocess we should calculate the squared modulus of the amplitude (13) summed over polarizations of the initial gluons by means of the following relation: Then we obtain the general form of the gg → η c η c cross section corresponding to the production amplitude (13): Making the substitutions for the functions A i , we find it useful to transform the result (15) as follows: The auxiliary functions F (i) entering the cross section (16) are written explicitly in Appendix. Note that the function F (0) describes non-relativistic result which coincides in the limit M ηc = 2m with the corresponding function obtained in Ref. [11] in the approach of NRQCD. Relativistic corrections in (16) are determined by a number of relativistic parameters ω i : In the non-relativistic limit, the parameterR(0) coincides with the definition of radial wave function at the origin, so it can be considered in some way as its relativistic generalization. All parameters, which contain the meson wave functions and describe the transition of the pairs (cc) to the bound state, are calculated in the framework of relativistic quark model [17,18]. This model is based on the Schrödinger equation with the Breit Hamiltonian in QCD and the nonperturbative confinement terms. Using the program of numerical solution of the Schrödinger equation [15,23], we obtain relativistic wave functions and bound state energies of S-wave charmonia. Numerical values of charmonium masses M th J/ψ = 3.072 GeV  [24]. The additional details on our relativistic quark model can be found in Refs. [15,18].
The numerical results of our calculation of the pair S-wave charmonium production cross sections in the case of non-relativistic approximation as well as with the account of relativistic corrections of order v 2 are presented in Table I. Along with total cross section values, we have also included there the cross section predictions corresponding to the rapidity interval 2 < y P,Q < 4.5 of the LHCb experiment [1] calculated with two different sets of linear PDFs: CTEQ5L [25] and CTEQ6L1 [26]. As shown in Table I, the cross section σ[pp → 2η c + X] at √ S = 7 (14) TeV is equal to 1.3 (2.4) nb for CTEQ5L and 1.0 (1.7) nb for CTEQ6L1. The most important production rates lie in the region of small P T (see Fig. (3)), where the color singlet contribution is dominant. Performing the numerical integration of differential cross section (16), we use the LO expression for the running coupling constant α s (µ) with the initial value α s (M Z ) = 0.118 and the renormalization scale µ = m T = M 2 + P 2 T , where M is the meson mass. In our numerical calculations of the cross sections we set M ηc = 2.980 GeV and m = 1.55 GeV. Therefore, we take into account non-zero bound state energy of η c charmonium state leading to the bound state corrections to the production cross section (16). Numerical results in Table I are determined by a number of parameters and functions: the c-quark mass, the factorization scale µ, parameters of the quark interaction operator, the bound state wave functions, the parton distribution functions and the strong coupling constant. Some of them (the c-quark mass, the quark-antiquark potential) are fixed in the relativistic quark model in the mass spectrum calculation. The factorization scale µ is taken in a commonly used form [8,10,11]. Other quantities lead to basic uncertainties of our numerical results.
It is evident from Table I, that relativistic corrections of order v 2 decrease the cross section values more than two times in both cases connected with the pair production of J/ψ or η c mesons. The only exception is the case of a pair η c production cross section in the rapidity region 2 < y P,Q < 4.5, where the relativistic effects decrease the cross section only by approximately 20 percents. The change of PDF from CTEQ5L to CTEQ6L1 brings the additional 20÷30 percent decreasing to the value of cross section. Along with the possibility of different PDF choices, there also exists the uncertainty dealt with the determination of every particular partonic distribution function. The sets CTEQ5L and CTEQ6L1 contain no means to estimate the uncertainties of such sort, however the set CTEQ6M [26] has all necessary functionality. Using 40 uncertainty eigenvectors from CTEQ6M we can roughly estimate the error of every cross section value in Table I dealt with the PDF uncertainty as 15 %. The only known calculation of the η c -pair production in proton-proton collision was performed in Ref. [11]. Their Table II contains the obtained numerical results for different PDFs CTEQ5L and CTEQ6L1 with P T > 3 GeV, which are of order of 4 nb. They used almost the same values of the c-quark mass, the factorization scale µ as in our calculation but a different numerical value for the parameter R(0). Our value of the radial wave function at the origin in nonrelativistic limit is equal to R(0) ≈ 0.8 GeV 3/2 , whereas in Ref. [11] the authors took the long distance matrix element O 1 S = 1.4 GeV 3 , which gave R(0) = 1.7 GeV 3/2 . In the region of large transverse momentum P T (P T > 3 GeV) the cross section falls considerably, so that the basic contribution to our result in Table I is determined by the region of small momenta P T . Therefore, our nonrelativistic results 1.5 nb and 1.2 nb differ significantly from the values of cross sections obtained in nonrelativistic SPS approximation in [11]. This difference is related with a choice of the parameter R(0) in [11] which exceeds our value more than two times.
Another possible source of uncertainties is connected with the determination of relativistic wave function in the momentum region p > ∼ m. The obtained charmonium wave function is strongly decreasing in this region. Its numerical value at p = m is more than 50 times smaller the maximum value. Nevertheless, relativistic factors p 2 and p 4 , entering in the integrals I 1,2 change this relation and increase the inaccuracy in the wave function determination at p > ∼ m. In spite of the fact that momentum integrals appear to be fully convergent, our relativistic model cannot provide a reliable calculation of the wave functions in the region of relativistic momenta p > ∼ m. Our definitions (17) of the parameters I 1,2 contain the cutoff at relativistic momentum of order m. Using indirect arguments related with the mass spectrum calculation accuracy we estimate in 10% the uncertainty of the wave function determination. Larger value of the error would lead to the essential discrepancy between the experiment and theory in the calculation of the charmonium mass spectrum. Then the corresponding error in the cross section (16) is not exceeding 20%. We do not consider a part of theoretical error related with radiative corrections of order α s because these corrections are omitted in our analysis. We also assume that relativistic corrections of order O(v 4 ) to the cross section (16) coming from the production amplitude should not exceed 30% of the obtained relativistic result. So, our total theoretical error is not exceeding 39%. To obtain this estimate we add the above mentioned uncertainties in quadrature.
Then, the coefficients F (i) entering the differential cross section (16) have the following form: a (a 2 1 + a 2 2 ), F (1) = k a k b (a 1 b 1 + a 2 b 2 ) + 8k 2 a (a 2 1 + a 2 2 ), F (2) = −4k 2 a (a 2 1 + a 2 2 ), F (3) = 6k a k b (a 1 b 1 + a 2 b 2 ) + k a k c (a 1 c 1 + a 2 c 2 ) + 24k 2 a (a 2 1 + a 2 2 ) + Note that we expand functions a 1,2 and b 1,2 in the mass difference (2m − M) up to the term linear in (κ − 1/2) and set the value κ = 1/2 in c 1,2 . Such simplifications allow us to significantly reduce the length of analytical expressions (A2)-(A7), while the numerical results of the cross sections change on 1 ÷5 percents. In Table I we present numerical results corresponding to exact functions a i , b i , and c i .