Relativistic corrections to prompt double charmonium hadroproduction near threshold

We calculate the relativistic corrections to prompt $J/\psi$ pair and $J/\psi+\psi(2S)$ hadroproduction through the color-singlet channel within the framework of nonrelativistic QCD (NRQCD) factorization. The short-distance coefficients are obtained by matching full-QCD and NRQCD calculations at the partonic level, in which both squared amplitude and phase space are expanded in $v^2$. We find that such an expansion of the phase space spoils the convergence of NRQCD factorization near the production threshold. To fix this problem, we propose to modify the matching between full QCD and NRQCD by adopting the physical phase space. In this modified approach, the theoretical uncertainties due to the choice of charm quark mass $m_c$ are largely reduced and the overall agreement of our predictions with LHC data is significantly improved, both as for total and differential cross sections.


I. INTRODUCTION
The production mechanism of heavy quarkonium, the bound state of a heavy quarkantiquark pair, has been a long-standing puzzle in high-energy particle physics ever since the discovery of the first charmonium state J/ψ about half a century ago.Due to the hierarchy of energy scales in heavy-quarkonium production, the nonrelativistic QCD (NRQCD) factorization approach [1], which was established on top of the NRQCD effective field theory [2], has been generally viewed as the most suitable tool to study this process.In NRQCD factorization, the cross sections of heavy-quarkonium production are factorized into products of short-distance coefficients (SDCs) and long-distance matrix elements (LDMEs).The SDCs describe the creation of the heavy-quark pair and can be calculated perturbatively in powers of the strong-coupling constant α s .The supposedly universal LDMEs represent the non-perturbative evolution of heavy-quark pairs into the heavy-quarkonium states and scale as powers of v according to velocity scaling rules [3], where v is the relative quark velocity in the center-of-mass frame.In this way, the theoretical predictions are organized as double expansions in α s and v.A key feature of NRQCD factorization is the incorporation of all possible Fock states, including both color singlet (CS) and color octet (CO) ones.The possible dominance of CO states, dubbed CO mechanism, greatly improves our understanding of the heavy-quarkonium production.However, the verification of LDME universality in global fits at next-to-leading order (NLO) to experimental data still poses a serious challenge, especially in the most important case of the J/ψ meson (see Ref. [4] for a recent review).
In recent years, J/ψ associated hadroproduction with another quarkonium state, a W , or a Z boson [5] have attracted considerable attention.Among them, prompt J/ψ pair production is of particular interest not only because of doubled LDME sensitivity, but also because this provides a laboratory for the study of double parton scattering (DPS), which is important at large invariant mass m ψψ and large rapidity separation |∆y|, and the extraction of the key DPS parameter σ eff [6].Moreover, this is the background for the newly discovered tetra-quark state X(6900) [7][8][9], which mainly decays into J/ψ pairs.On the experimental side, prompt double J/ψ hadroproduction has been measured by the D0 Collaboration [10] at the Fermilab Tevatron and by the LHCb [11,12], CMS [13], and ATLAS Collaborations [14] at the CERN LHC.Under the assumption that the single parton scattering (SPS) contribution can be neglected in the large-|∆y| region, both the D0 and LHCb Collaborations were able to extract σ eff , yet with different results.Recently, the LHCb Collaboration meticulously measured total and differential cross sections of inclusive 2J/ψ [15] and J/ψ+ψ(2S) [16] hadroproduction at center-of-mass energy √ s = 13 TeV, separately for SPS and DPS in the 2J/ψ case, enabling meaningful comparisons with NRQCD predictions.
On the theoretical side, the idea to probe LDMEs in double J/ψ hadroproduction was first raised by Barger et al. in Ref. [17], where the fragmentation contribution to the 2cc 1 ) production channel was considered.Later on, it was found that the leading-order (LO) 1 ) mainly contributes at small and moderate values of transverse momentum p T [18][19][20], where it is largely enhanced at NLO [21], while in the large-m ψψ and large-|∆y| regions, contributions involving t-channel gluon exchange dominate by orders of magnitude [22].Upon including all possible LO contributions, the SPS predictions in the large-m ψψ and large-|∆y| regions were found to underestimate experimental measurements by more than one order of magnitude [22].To reduce the large gap, higher-order QCD corrections to some channels [23][24][25] were calculated, and large logarithm of the form (α s ln |s/t|) n were resummed using high-energy factorization [23].Notice that, in the case of double P -wave Fock state production, a new type of infrared divergences causes NRQCD factorization as we know it to break down [26] pushing complete NLO NRQCD predictions out of reach for the time being.Besides experimental determinations [10][11][12], also theoretical groups, armed with advanced knowledge of SPS contributions, fitted σ eff to experimental data finding different results [27,28].For a recent review of theoretical progress on prompt double J/ψ hadroproduction, we refer to Ref. [29].
While NRQCD predictions for prompt double J/ψ hadroproduction were generally found to undershoot the experimental data in the large-m ψψ region, the situation is very different close to the production threshold, m ψψ 2m J/ψ [22].In fact, the LHCb measurement at center-of-mass energy √ s = 7 TeV undershoots the NRQCD prediction by almost one order of magnitude in the first bin, 6 GeV < m ψψ < 7 GeV [22].Close to threshold, the CO contribution is negligibly small and the QCD-corrected CS channel gg → 2cc( 3 S [1] 1 ) significantly falls short of the LHCb data.
Besides NLO QCD corrections, higher-order relativistic corrections may also have an important impact given that v 2 ≈ α s (2m c ) for charmonium.This is supported, e.g., by studies of J/ψ associated production with η c mesons, charmed [30] and light hadrons [31,32] in e + e − annihilation and by J/ψ yield [33,34] and polarization [35] in single inclusive photoand hadroproduction.
Higher-order relativistic correction to inclusive double J/ψ hadroproduction were first studied in the context of the relativistic quark model and were found to reduce the LO total cross section by a factor of 6, down to a value undershooting the LHCb measurement [11] by more than a factor of 3 [36].In the analogous study for differential cross sections in the NRQCD factorization framework [37], the relativistic corrections were found to be just a few percent, and it was argued that their effect on the phase space is suppressed at the hadron level by an additional factor of 1/ ln( √ s/m T ), where m T = m 2 J/ψ + p 2 T is the J/ψ transverse mass.This ambiguous situation is unsatisfactory and together with the new high-precision data from LHCb [15,16] motivates us to revisit Ref. [37], which has not been confirmed yet.In the following, we thus perform an independent calculation of the NLO relativisitic corrections to inclusive double charmonium hadroproduction via SPS in the NRQCD factorization framework, placing special emphasis on the threshold region, which actually requires a careful treatment of the phase space.
The remainder of this paper is organized as follows.In Sec.II, we describe the calculation of SCDs in detail.In Sec.III, we present our numerical results for the NRQCD predictions including relativistic corrections, both for 2J/ψ and J/ψ + ψ(2S) production, and compare them with respective LHCb data [15,16].Our conclusions are summarized in Sec.IV.

II. ANALYTIC RESULTS
According to the factorization theorem of the collinear parton model, the hadroproduction rate of inclusive double prompt J/ψ is given by where f i/A (x) is the parton distribution function (PDF) of parton i in hadron A. In NRQCD factorization through relative order O(v 2 ), the partonic cross section can be further factor-ized as where O H (m) and P H (m) are the LO and O(v 2 ) four-quark operators pertaining to the are the corresponding SDCs of the partonic subprocess i + j → cc(m) + cc(n) + X, and Br(H → J/ψ + X) is the branching fraction, which equals unity for H = J/ψ.Near threshold, the J/ψ pair is mainly produced through the CS channel of the gluon fusion 1 ), on which we concentrate in the following.This also includes the feed-down contribution from ψ(2S).
The relevant four-quark operators for the production of H = J/ψ, ψ(2S) are defined as [1] O where ψ (χ) is the Pauli spinor that annihilates (creates) a heavy quark (antiquark), σ i In NRQCD factorization, the SDCs are the same for the J/ψ and ψ(2S) mesons, since the Fock states carry the same quantum numbers.Therefore, we only need to calculate 1 ) perturbatively, through the matching condition between full QCD and the NRQCD, The left-hand side of Eq. ( 4) can be computed easily by implementing the covariantprojection method [38].In this method, the product of Dirac spinors v(p c)ū(p c ) is written in a Lorentz-invariant form to all orders of v 2 after projection onto spin singlet, S = 0, or spin triplet, S = 1, configurations.The four-momenta of the c and c quarks of the (cc) pair in an arbitrary reference frame can be related to those in the (cc) rest frame via an appropriate Lorentz boost matrix L µ ν as where P µ r = (2E q , 0), E q = q 2 + m 2 c , and 2q is the relative three-momentum of the c and c quarks in the (cc) rest frame.Then, in an arbitrary frame, the projection on spin triplet reads [38] where ǫ µ is the polarization vector of the spin triplet state and the Dirac spinors are normalized as ūu = −vv = 2m c .In this way, the full QCD scattering amplitude can be expressed as where the factor m c /E q i compensates the relativistic normalization of the (cc) i pair and with being the standard Feynman amplitude and 3, space for the cc pair to be projected onto the CS state.
To expand Eq. ( 8) as a double series in q 1 and q 2 , it is convenient to define so that For S-wave production, only even powers of q i contribute.To extract the higher-order corrections in v 2 , we need to further decompose the higher-rank tensor products into the irreducible representations and retain the S-wave terms.At O(v 2 ), we have Through O(v 2 ), we thus have In Eq. ( 12), the second and third terms within the square brackets are already of O(v 2 ).
However, the first term A(0, 0) still depends on q 2 1 and q 2 2 implicitly via where k 1 and k 2 are the four-momenta of the incoming gluons.Notice that s is independent of q 2 1 and q 2 2 .To expand t and u in v 2 , we find it is more convenient to introduce two new variables, The two-body phase space measure thus takes a very simple form, In the nonrelativistic limit, q 2 1,2 → 0, we have t → t 0 and u → u 0 with s + t 0 + u 0 = 8m 2 c .Analogously to Eq. ( 14), we introduce It is straightforward to see that where Thus, the two-body phase space integral in Eq. ( 15) can also be re-expressed in t0 and û0 as With the help of Eqs. ( 14), (16), and ( 17), both A(0, 0) and the two-body phase space measure in Eq. ( 19) can be easily expanded in v 2 , too.
We are now in a position to present the master formulas for the SDCs in the matching condition in Eq. ( 4), with , where the symbol implies summing over all polarizations in the initial and final states and then averaging the spins and colors of the incoming partons and the outgoing cc Fock states.The factor K = −4/(s − 16m 2 c ) originates from the expansion of the k factor in Eq. (19) and is divergent at threshold, although the complete phase space integral is still convergent.
From the explicit expression of the k factor, we observe that the actual expansion parameter for the phase space is not v 2 = q 2 /m 2 c , but rather 8m 2 c v 2 /(s − 16m 2 c ), which is not small near threshold.We hence conclude that, for J/ψ pair production near threshold, the usual fixed-order NRQCD predictions will not be reliable unless the numerical value of 8m 2 c /(s − 16m 2 c ) is of order unity or smaller.To overcome this problem, we alternatively choose the physical phase space in our numerical evaluation, without the expansion in v 2 , and replace t 0 and u 0 by their physical counterparts, t p and u p , using a similar relation as Eq. ( 17).The numerical consequence of this treatment is addressed in Sec.III.
We use the program packages Feyn Arts [39] to generate Feynman diagrams and Feyn Calc [40] to handle Dirac algebra, color indices, and expansion of the squared amplitude.When we drop the relativistic corrections to the phase space, our result for 1 ) agrees with that in Ref. [37] except for a factor of 1/2 due to the identical-boson symmetry.We employ the program package CUBA [41] to perform the numerical phase space integrations.

III. NUMERICAL RESULTS
In the numerical analysis, we adopt LO set CTEQ6L1 [42] of proton PDFs, available through the LHAPDF6 program library [43], and the LO formula for α and Λ (4) QCD = 215 MeV [42].We set the renormalization and factorization scales to be for the direct production of charmonia H 1 and H 2 with common transverse momentum p T .We take ξ = 1 as default and vary ξ between 1/2 and 2 to estimate the theoretical uncertainties from unknown higher orders in α s .For definiteness, we choose m c = 1.5 GeV as is frequently done in similar analyses, so as to facilitate comparisons with the literature.The values m J/ψ = 3.097 GeV, m ψ(2S) = 3.686 GeV, Br(ψ(2S) → J/ψ + X) = 61.4% for masses and branching fraction are taken from the latest Review of Particle Physics [44].As for the J/ψ and ψ(2S) CS LDMEs, we use the LO ones evaluated from the wave functions at the origin for the Buchmüller-Tye potential [45], and estimate the O(v 2 ) ones by the ratio obtained for the J/ψ case in Ref. [46], In the evaluation of feed-down contributions to differential cross sections, we approximate the momentum of the J/ψ meson from ψ(2S) decay as We are now in a position to confront our improved theoretical predictions for the cross section of prompt double J/ψ hadroproduction in SPS with experimental data to investigate how the inclusion of O(v 2 ) relativistic corrections affects the agreement close to the production threshold.Our prime observable is the m ψψ distribution in the first few bins.Owing to the LO relation m ψψ = 2 4m 2 c + p 2 T cosh(∆y) [22], this corresponds to small values of rapidity separation |∆y|, so that we also focus on the first few bins of the |∆y| distribution.
Because the m ψψ distribution rapidly falls off with increasing value of m ψψ [22], the total cross section is dominated by the threshold region.As mentioned above, the bulk of the cross section close to threshold is due to the CS channel gg → 2cc 1 ), while the CO contributions are suppressed there [22].In fact the CO contributions to the total cross sections of direct and prompt double J/ψ hadroproduction under LHCb experimental conditions at √ s = 7 TeV [11] only amount to 8.2% and 5.9%, respectively [22].
Of all available experimental data [10][11][12][13][14][15][16], we select for our dedicated comparisons the recent LHCb data of Ref. [15] because they are most precise and the SPS contribution is separated from the DPS one, which is not included in our predictions.Reference [15] is based on a data sample corresponding to an integrated luminosity of 4.2 fb −1 , refers to prompt double J/ψ hadroproduction at √ s = 13 TeV, with both J/ψ mesons in the ranges 0 < p T < 14 GeV and 2.0 < y < 4.5, and specifies total cross sections and various distributions.
Besides the m ψψ and |∆y| distributions in their lower ranges, also the distribution in the rapidity y of either J/ψ meson is of interest for our purposes because it is also dominated by the small-m ψψ region.In all these cases, LO NRQCD, at O(α 4 s ), is expected to yield useful approximations, to be improved by the inclusion of relativistic corrections of relative order O(v 2 ).In the following, we only consider SPS data from Ref. [15].We start by considering the total cross section and compare the LHCb result [15], where the first error is statistical and the second one systematic, with our LO and O(v 2 )corrected NRQCD predictions, respectively.As is familiar from similar comparisons in the literature [29], the LO NRQCD prediction greatly overshoots the experimental result, by more than a factor of 2. This overshoot is significantly reduced by the inclusion of the O(v 2 ) corrections, to about 50%. 1 ), and the residual contribution, from the phase space.Detailed analysis reveals that the latter is dominant in the first m ψψ bin, while G 1,2 ( 3 S We note in passing that inclusion of our relativistic corrections also leads to improved descriptions of prompt double J/ψ hadroproduction at √ s = 7 TeV [11] and prompt J/ψ + ψ(2S) hadroproduction at √ s = 13 TeV [16] as measured by the LHCb Collaboration, albeit without SPS-DPS separation.To illustrate this, we merely consider the total cross sections.
Specifically, we have for Ref. [11] for Ref. [16].It is reasonable to expect that the SPS fractions of the LHCb results in Eqs. ( 28) and ( 29) are in the same ball park as for the LHCb results in Ref. [15], where it is about 50%, as may be inferred by comparing Eq. ( 26) with [15] σ SPS+DPS LHCb = 16.36 ± 0.28 ± 0.88 nb .
As explained in Sec.II, the conventional expansion of the phase space in v 2 diverges in the limit when the two-gluon invariant mass √ s approaches the four-quark threshold 4m c , which renders the naive evaluation of the relativistic corrections ill-defined for m ψψ values close to 2m J/ψ .To overcome this problem, we proposed to evaluate the phase space using the meson mass m J/ψ instead of the quark mass m c and not to expand it in v 2 .
We illustrate the numerical effect of this modification for the m ψψ distribution of direct double J/ψ hadroproduction in SPS in Fig. 2. We observe from Fig. 2, that the naive evaluation leads to a negative result of similar magnitude in the first m ψψ bin, violating NRQCD factorization.This reflects the fact that, in the treatment of the phase space, the formal expansion parameter v 2 is endowed with the factor 8m 2 c /(s − 16m 2 c ), which is large close to threshold.On the other hand, the relative difference between the naive and modified evaluations rapidly decreases with increasing value of m ψψ and becomes insignificant beyond 10 GeV.A rigorous solution of this problem is likely to require a resummation to all orders of v 2 , which is beyond the scope of this work.We refer the interested reader to Ref. [47], where such a resummation was discussed for the exclusive process e + e − → J/ψ + η c .
We found that the conventional power counting for the expansion in v 2 is offset for the phase space close to threshold due to an enhancement factor of the form 8m 2 c /(s−16m 2 c ), where √ s is the partonic center-of-mass energy.As a cure to this problem, we proposed an alternative treatment of the phase space implemented with physical charmonium masses.
We confronted our improved theoretical predictions with recent high-precision measurements by the LHCb Collaboration, in which the SPS contributions were extracted by appropriate acceptance cuts [15].In fact, the inclusion of O(v 2 ) corrections turned out to significantly improve the NRQCD description of the measured m ψψ distribution in the first few bins.As expected, this improvement was found to carry over to the |∆y| distribution in the first few bins, to the y distribution, and to the total cross section, which all receive dominant contributions from the near-threshold region of the m ψψ distribution.
In conclusion, relativistic corrections play an important role in our understanding of the mechanism underlying double charmonium hadroproduction at the LHC, especially in kinematic configurations close to the production threshold, and enhance the predictive power of NRQCD factorization.Further improvements of the theoretical description near threshold are expected from the resummations of relativistic corrections [47] and logarithmic corrections due to multiple soft-gluon emission [48], which lie beyond of the scope of this work.
space component of the covariant derivative D µ .

FIG. 1 .
FIG. 1.The (a) m ψψ , (b) |∆y|, and (c) y distributions of prompt double J/ψ hadroproduction in SPS measured by LHCb [15] at √ s = 13 TeV are compared with our CS NRQCD predictions at LO (dashed lines) and NLO with respect to relativistic corrections (solid lines).The theoretical uncertainties are indicated by the shaded blue (LO) and yellow (NLO) bands.

1 )
wins out in higher bins.In the first few m ψψ bins, the inclusion of the O(v 2 ) corrections significantly improves the agreement between NRQCD predictions and LHCb data.The same is true for the |∆y| and y distributions in Figs.1(b) and (c), respectively.In fact, as for the first three bins of the |∆y| distribution and most of the bins of the y distribution, the errors of the O(v 2 )-corrected NRQCD prediction and the LHCb data overlap.

FIG. 2 .
FIG. 2. The m ψψ distribution of direct double J/ψ hadroproduction in SPS under LHCb experimental conditions evaluated in NRQCD with O(v 2 ) relativistic corrections evaluated in our modified approach (solid lines) is compared with its counterpart evaluated by rigidly expanding the phase space in v 2 (dashed lines).The theoretical uncertainties are indicated by the shaded blue (rigi-NRQCD) and yellow (our NRQCD) bands.Negative values are shown in red.