The exclusive $J/\psi$ process at the LHC tamed to probe the low $x$ gluon

The perturbative QCD expansion for $J/\psi$ photoproduction appears to be unstable: the NLO correction is large (and of opposite sign) to the LO contribution. Moreover, the predictions are very sensitive to the choice of factorization and renormalization scales. Here we show that perturbative stability is greatly improved by imposing a $`Q_0$ cut' on the NLO coefficient functions; a cut which is required to avoid double counting. $Q_0$ is the input scale used in the parton DGLAP evolution. This result opens the possibility of high precision exclusive $J/\psi$ data in the forward direction at the LHC being able to determine the low $x$ gluon distribution at low scales.


Introduction
It would be valuable to be able to constrain the gluon parton distribution function (PDF) at low x using J/ψ photoproduction data measured at HERA and at the LHC, via exclusive pp → p + J/ψ + p events, especially events in the forward region measured by the LHCb collaboration. Indeed, for LHCb kinematics at 13 TeV we can reach down to x 3 × 10 −6 . Exclusive J/ψ production is driven by the subprocess γ * p → J/ψ +p, see Fig. 1. Unfortunately, it turns out that the NLO corrections calculated in the conventional MS collinear approach are found to be very large and to depend strongly on the choice of factorization and renormalization scales [1,2,3]. Indeed, for an 'optimum' choice of scales it is found that the NLO correction has the opposite sign to the LO contribution and even changes the sign of the whole amplitude, see the continuous curves in Fig. 2. Thus one may doubt the convergence of the whole perturbation series.

Optimum scale
What do we mean by the 'optimum' scale? It was shown in Ref. [3] that it is possible to find a scale (namely µ F = m c ) which resums all the double logarithmic corrections enhanced by large values of ln(1/ξ) into the gluon and quark PDFs, where ξ is the skewedness parameter of the Generalised Parton Distributions (GPDs) describing the proton-gluon (and proton-quark) vertices. That is, it is possible to take the (α S ln(1/ξ)ln(µ 2 F )) term from the NLO gluon (and quark) coefficient functions and to move it to the LO GPDs. This allows a resummation of all the double log (α S ln(1/ξ)ln(µ 2 F )) n terms in the LO contribution by choosing the factorization scale to be µ F = m c . The details are given in Ref. [3], see also Ref. [5].
The result is that the γp → J/ψ + p amplitudes are schematically of the form where the GPD can be related to the conventional PDF via the Shuvaev transform for ξ < |x| 1 [6]. With the choice µ F = m c there is a smaller remaining term in the NLO coefficent funcions, and so the residual dependence on the scale µ f is reduced.
Unfortunately, even after this, the NLO corrections, and their variations with scale, although reduced, are still unacceptably large, as shown in Fig. 2. The dashed and dot-dashed curves correspond to NLO predictions for two different values of the residual scale µ f : namely µ 2 f = 4.8 and 1.7 GeV 2 respectively, while the continuous curves correspond to the 'optimum' scale choice µ 2 The choice µ R = µ F is justified in subsection 3.1.

Double counting
So for the QCD prediction to be useful we should search for some other sizeable physical contribution to the NLO correction. Here we consider a power correction which may further reduce the NLO correction and, moreover, may reduce the sensitivity to the choice of scale. The correction is O(Q 2 0 /M 2 ψ ) where Q 0 denotes the input scale in the parton evolution. It turns out to be important for the relatively light charm quark, m c M ψ /2. Let us explain the origin of this 'Q 0 correction'. We begin with the collinear factorization approach at LO. Figure 3: (a) LO contribution to γp → V + p. (b) NLO quark contribution. For these graphs all permutations of the parton lines and couplings of the gluon lines to the heavy-quark pair are to be understood. Here P ≡ (p + p )/2 and l is the loop momentum.
Here, we never consider parton distributions at low virtualities, that is for Q 2 < Q 2 0 . We start the PDF evolution from some phenomenological PDF input at Q 2 = Q 2 0 . In other words, the contribution from |l 2 | < Q 2 0 of Fig. 3(b) (which can be considered as the LO diagram, Fig. 3(a), supplemented by one step of DGLAP evolution from quark to gluon, P gq ) is already included in the input gluon GPD at Q 0 . That is, to avoid double counting, we must exclude from the NLO diagram the contribution coming from virtualities less than Q 2 0 . At large scales, Q 2 Q 2 0 this double-counting correction will give small power suppressed terms of O(Q 2 0 /Q 2 ), since there is no infrared divergence in the corresponding integrals. On the other hand, with Q 0 ∼ 1 GeV and µ F = m c (∼ M ψ /2) a correction of O(Q 2 0 /m 2 c ) may be crucial. In the present paper we re-calculate the NLO contribution for J/ψ photoproduction excluding the contribution coming from the low virtuality domain (< Q 2 0 ). We find that for J/ψ this procedure substantially reduces the resulting NLO contribution and, moreover, reduces the scale dependence of the predictions. It indicates the convergence of the perturbative series.
An outline of the procedure is given in [9], where also the NLO description of exclusive J/ψ production in the k T factorization and collinear factorization schemes are compared.
2 Avoiding double counting of the low Q 2 contribution

The NLO quark contribution
We start with the NLO quark contribution to the γp → J/ψ + p process. The corresponding Feynman diagrams are that of Fig. 3(b) together with the diagram where both gluons couple to the same heavy quark line. Here we will use the non-relativistic approximation for the J/ψ wave function. Since the momentum fractions (x + ξ) and (x − ξ) carried by the left and right quarks are different we have to use the skewed (generalized) parton distribution (GPD), F q (x, ξ, Q 2 ). The skewedness parameter ξ = M 2 ψ /(2W 2 − M 2 ψ ), where W is the γp energy. We see that the upper part of diagram Fig. 3(b) is the same as the diagram for the LO gluon Fig. 3(a) contribution. For the LO contribution the integral over the gluon virtuality |l 2 | starts from the input scale Q 2 0 , while all the contributions from low virtualities |l 2 | < Q 2 0 are collected in the input gluon GPD, F g (x, ξ, Q 2 0 ). Note that this input distribution already includes that part of the quark contribution of Fig. 3(b) coming from |l 2 | < Q 2 0 . Thus to avoid double counting when computing the NLO quark coefficient function, C NLO q , of Fig. 3(b) we have to include the theta function Θ(|l 2 | > Q 2 0 ) in the integration over l 2 . Depending on the ratio this can be a significant correction. The corresponding integral has no infrared or ultraviolet divergence and can be calculated in D = 4 dimensions.
Actually, the calculation is performed in the physical scheme (with D = 4). On the other hand, parton distributions are usually presented in the MS factorization scheme where dimensional regularization is used. The problem is that when we calculate the coefficient function in D = 4 + 2 we have finite contributions of / origin. Formally these / terms come from unphysically large distances ∝ O(1/ ). In fact, these / terms are compensated by a corresponding re-definition of the PDFs. In order to retain the / terms and to use the MS scheme we do not calculate diagram 3(b) in D = 4 dimensions for |l 2 | > Q 2 0 , but instead calculate the part corresponding to small |l 2 | < Q 2 0 . We consider this part as the correction which should be subtracted from the known NLO MS coefficient function [1,10]. Recall that after the subtraction of the contribution generated by the last step of the LO evolution, P LO ⊗ C LO , there is no infrared divergence and the subtracted part of C NLO coming from |l 2 | < Q 2 0 does not contain / terms.

The NLO gluon contribution
The NLO 'Q 0 corrections' for the gluon coefficient function are more complicated. Besides the ladder-type diagrams analogous to Fig. 3(b), but with the light quark line replaced by a gluon line, there are other diagrams which have a structure similar to vertex corrections, see [1,10]. However the 'dangerous' contribution is again from the ladder-type diagrams, where to avoid double counting we have to exclude the |l 2 | < Q 2 0 domain whose contribution is already included in the LO term using the input gluon GPD, F g (x, ξ, Q 2 0 ). Qualitatively this is exactly the same calculation as that for the NLO quark. The only difference is that the lower line in the diagrams of Fig. 5 is now replaced by a gluon line and the lower part of the diagram is now given by the product of two three-gluon vertices averaged over the incoming gluon transverse polarizations. The notation is identical to that for the quark contribution. Both the quark-and the gluon-induced contributions are determined as described in the Appendix. They involve the calculation of the diagrams of Fig. 5 (given in the Appendix), and the analogous diagrams for the gluon-induced contribution.  3 Results Fig. 4 shows the LO and NLO contributions to the imaginary part of the J/ψ photoproduction amplitude when the Q 0 cut in the NLO contribution is taken into account. It should be compared to Fig. 2 which had exactly the same scale choices, but without the Q 0 cut imposed. The improvement in going from Fig. 2 to Fig. 4 is dramatic. First, the NLO contribution is now much smaller than the LO contribution. Second, the scale variation is much smaller. The continuous curves in Figs. 2 and 4 show the LO and NLO comparison for the choice of scales µ F = µ R = m c ≡ M ψ /2, which we had previously argued was optimal [3]. The stability achieved by imposing the Q 0 cut means that J/ψ photoproduction (γp → J/ψ p) data and LHC exclusive J/ψ (pp → p + J/ψ + p) data can now be included in the global parton analyses.

The choice of scales
Let us discuss the above scale choices in more detail. By choosing the 'optimal' factorization scale µ F = m c we resum all the higher-order double-logarithmic corrections (α s ln(1/ξ) ln µ 2 F ) n (enhanced at high energies by the large value of ln(1/ξ)) into the gluon generalized parton distribution (gluon GPD) [3].
The renormalization scale is taken to be µ R = µ F . The arguments are as follows. First, this corresponds to the BLM prescription [11]; such a choice eliminates from the NLO terms the contribution proportional to β 0 (i.e. the term β 0 ln(µ 2 R /µ 2 F ) in eq. (3.95) of [1]). Second, following the discussion in [12] for the analogous QED case, we note that the new quark loop insertion into the gluon propagator appears twice in the calculation. The part with scales µ < µ F is generated by the virtual component (∝ δ (1 − z)) of the LO splitting during DGLAP evolution, while the part with scales µ > µ R accounts for the running α s behaviour obtained after the regularization of the ultraviolet divergence. In order not to miss some contribution and/or to avoid double counting we take the renormalization scale equal to the factorization scale, µ R = µ F .

Discussion of the results
Note that in the present paper we have calculated the imaginary part of the γp → J/ψ p amplitude. The real part of the amplitude can be restored via dispersion relations assuming positive signature, as in eq. (5) of Ref. [13]. Recall that we obtain the necessary GPDs from the CTEQ6.6 parton set [4] using the Shuvaev transform [6]. We use a relatively old parton set [4] in which the low x gluons are forced to be positive so as to make a meaningful comparison with our earlier work. The goal of this paper is not to make a quantitative description of the data, but to demonstrate that we can achieve stability of the perturbative QCD description of relatively low scale J/ψ production by imposing the Q 0 cut. We have shown this is a power correction -a correction which is needed to avoid double counting. This will allow future high precision exclusive J/ψ production data obtained at the LHC to be incorporated in global parton analyses.
The general procedure to include the HERA γp → J/ψ p data and, in particular, the LHCb data for exclusive J/ψ production, pp → p + J/ψ + p, in a global analysis follows that used to produce Fig. 4 of Ref. [13]. These processes are driven by the gluon PDF and the LHCb data probe the gluon at very low values of x. However, in Ref. [13] we approximated the NLO corrections to the coefficient functions by accounting for the explicit l ⊥ integration in the last step of the interaction. Moreover, we just fitted the J/ψ data and used a parametric form for the gluon which approximated its x and Q 2 dependence. So the analysis of Ref. [13] was quite simplified, although very informative; see, for example, Fig. 5 of [13] which compared the resulting gluon PDF with those of different global analyses 2 .
The present paper, on the other hand, retains collinear factorization and calculates the complete NLO contribution. We may expect the high γp energy, W , data points in the updated version of Fig. 4 of Ref. [13] to require a larger gluon distribution in the region from x < ∼ 10 −3 down to 10 −5 , at low scales, than coming from extrapolations of the NLO gluon PDFs from global fits to data not including the J/ψ data. An indication in favour of a larger gluon PDF in this domain comes also from the recent LHCb data on open charm (and beauty) [14].
Finally, it is useful to compare our approach with that of [15], where it was demonstrated that the re-summation of the BFKL-induced (α S ln(1/ξ)) n terms in the coefficient functions additionally reduces the factorization scale dependence. Recall that our choice of Figure 5: Two diagrams (a,b) computed for the NLO quark coefficient function. Note that p and p refer to the incoming and outgoing quark lines. In the corresponding diagrams computed for the NLO gluon coefficient function the light quark line is replaced by a gluon. The other two diagrams of the different coupling of the two t-channel gluons to the heavy quarks are implicitly included.
resums only the double logarithmic, (α S ln(1/ξ) ln µ F ) n contributions 3 . The remaining part, which does not contain ln µ F , should be considered, in the collinear factorization approach, as higher-order, NNLO, N 3 LO, ... corrections. Of course, it would be good to account for these corrections as well. However, to properly calculate these corrections one has to exclude the low (< Q 2 0 ) virtuality contribution. Otherwise we will face the problem of double counting again. The present paper shows these (power) corrections (necessary to avoid double counting) are crucial to achieve perturbative stability.

Appendix
Here we describe the calculation of the piece that we subtract from the full result. Only the imaginary part of the ladder-type cut diagrams shown in Fig. 5 and the corresponding diagrams where the light-quark line is replaced by gluons is computed.
All momenta appearing in the calculation may be decomposed in terms of light-like momenta p, n and a transverse four-momentum l ⊥ , where l is the loop momentum and h 1 , h 2 are the momenta of the outgoing heavy quark and heavy anti-quark, respectively. Here p can be chosen as the momentum of the incoming light parton and n the momentum of the incoming on-shell photon. With this convention we have p · p = n · n = 0, p · n =ŝ/2, p · l ⊥ = n · l ⊥ = 0, whereŝ is the photon-parton centre-of-mass energy squared. The four momenta of the incoming and the outgoing light partons are proportional. We may write p µ and p µ = Xp µ with To leading order in the heavy quark relative velocity, the S-wave spin-triplet component of J/ψ can be computed using the projection [16,17,18] Hereū, v are the spinors of the outgoing heavy quark and anti-quark which form the J/ψ. The indices α and β label their spin. N J/ψ is an overall factor which contains the non-perturbative NRQCD matrix element describing the J/ψ formation. The vector J/ψ describes the polarisation of the J/ψ with momentum K = h 1 + h 2 and mass M ψ = 2m c .
In our calculation we split each diagram of Fig. 5 into two parts. An "upper" part which contains a trace over the heavy quark fermion line and a "lower" part which in the quark channel contains a trace over the light quark line and in the gluon case consists of two triple gluon vertices contracted with g µν ⊥ .
First we discuss the "upper" part which is different for the diagrams (a) and (b) of Fig. 5 but identical for the quark and gluon channels. Where it appears, we replace the contraction of l ⊥ with the polarisation vectors using which follows from tensor decomposing the l ⊥ integral after the integration over the l azimuthal angle. We can simplify the calculation by noting that the sum of the "upper" parts of diagrams (a) and (b) obey the gauge condition where Here T(h.loop) is the upper part of the amplitude, which besides the trace over the quark loop, includes the heavy quark propagator 1/(q 2 − m 2 c ). Using the gauge condition the only contractions of the "upper" part that appear in the sum of diagrams are Tr(h.loop) µν a p µ l ⊥ν = N J/ψ 4m c ( γ · * J/ψ )l 2 ⊥ŝ /2, Tr(h.loop) µν a l ⊥µ p ν = N J/ψ 4m c ( γ · * J/ψ )l 2 ⊥ŝ /2, The contractions involving p µ n ν , n µ p ν , n µ n ν , n µ l ν ⊥ , l µ ⊥ n ν , l µ ⊥ l ν ⊥ appear in the computation of individual diagrams but cancel for the sum of diagrams.

Quark-induced NLO correction
For an unpolarized light quark the trace over the "lower" light quark line gives where the normalization factor includes the colour factor C F and the quark GPD, F q .
This light quark part should be contracted with the trace, Tr(h.loop) µν , given by the heavy quark (upper) loop. Due to the gauge condition (11) we have that (p − l) µ acts as p µ , while (p − l) ν acts as p ν = Xp ν giving for diagram (a) and for diagram (b). The term M q accounts for terms which cancel between the two diagrams. The denominators come from the uncut propagators: 1/l 2 for the left, and 1/l 2 for the right gluon and 1/(q 2 − m 2 c ) for the uncut heavy quark propagator. The result is to be integrated over the gluon transverse momentum (dl 2 ⊥ ) while the longitudinal components are fixed by the quark on-mass-shell conditions. It is easy to perform this integral numerically accounting for the condition which was introduced in Section 2 in order to avoid double counting. Recall, however, that we are not going to calculate the whole NLO contribution, but just the correction to the known MS coefficient function. So, in order to compute the correction, in the integration over the l ⊥ we only consider the region of |l 2 | < Q 2 0 . Actually, we integrate over dl 2 directly; the factor (1 − β) coming from the relation l 2 = l 2 ⊥ /(1 − β) is exactly cancelled by the residue from the light quark on-mass-shell pole. So we obtain the correction to the quark-induced part of the γp → J/ψ + p amplitude where the 'hard matrix elements' M q a,b are given by (24) and (26). The factor 1/ŝ 2 comes from the delta functions needed to put the lower light quark and the heavy quark coupled to the right gluon in Fig. 5 on-mass-shell. The factor m 4 c accounts for the normalization N J/ψ , defined to be consistent with the normalization of eqs. (3.93) and (3.95) of [1] for which the correction was calculated; actually the last factor (...) is the correction to f q of (3.93) of [1]. 4 For the gluon correction ∆M g there is an additional factorŝ/2m 2 c = 1/ξ due to the definition of the gluon GPD, F g ; see the extra factor of ξ in eq. (3.94) of [1], see also [10].
Note that we have explicitly calculated the NLO diagrams (a) and (b) of Fig. 5 which contain both LO 5 and NLO contributions. To identify the NLO part we therefore have to subtract the contribution generated by the LO evolution equation, which is of the form of the convolution P LO ⊗ C LO , before we integrate over l 2 ⊥ . This subtraction completely cancels the logarithmic infrared divergence dl 2 /l 2 . Note that the subtraction must be done only in the region of |l 2 | < µ 2 F since at the factorization scale µ F the DGLAP evolution stops. 6 Also note that in the LO approximation the convolution P LO ⊗ C LO is larger than the value of the matrix element given by explicit calculation of the diagrams shown in Fig. 5. Thus the final result has the sign opposite to that for the LO amplitude.
In this way we obtain the quark NLO coefficient function. Since we are looking for the power correction needed to avoid double counting of the low |l 2 | < Q 2 0 contribution 7 , we actually have to integrate the matrix element M q over |l 2 | < Q 2 0 only (as explained above) and to subtract the result from the known NLO coefficient function given in the MS scheme.
In the notation of Ref. [1] this should be considered as the new form of Imf q (y) of their eq. (3.93), after allowing for the changes made by our introduction of the 'Q 0 cut'.

Gluon NLO correction
In the gluon case the tensor A g µν corresponding to the lower part of Fig. 5 diagrams (with the lower quark line replaced by a gluon line) was calculated explicitly. It can be written in the form A g µν = N g (ag µν + b 11 p µ p ν + b 22 h µ h ν + b 12 p µ h ν + b 21 h µ p ν + c 1 p µ l ⊥ν + d 1 l ⊥µ p ν + c 2 h µ l ⊥ν + d 2 l ⊥µ h ν ) , Here the normalization factor is 8 .
(30) 5 The integration of the pure logarithmic form dl 2 /l 2 up to µ F actually reproduces the LO contribution already included in Fig. 3(a). On the other hand some non-logarithmic corrections originating from higher powers of l 2 , together with the whole contribution above µ F , are NLO α s corrections which are not enhanced by the large collinear (l 2 ) logarithms. 6 This is the origin of the ln(4m 2 /µ 2 F ) factor in the first term of f q (y) of eq. (3.93) of [1]. Since now we integrate over the |l 2 | < Q 2 0 < µ 2 F the correction does not depend on µ F . 7 This contribution is already included in the input value GPD(Q 0 ). 8 Here the denominator (x + ξ − i )(x − ξ + i ) arises from the particular definition of the gluon GPD.
Note that X is defined in (4) and β is given by (7). Recall that we are looking for the imaginary part of the amplitude (i.e. s-channel discontinuity).
This expression should be convoluted with the "upper" part of the diagram. The result for the sum of diagrams can again be simplified using the gauge conditions (11). That is vector h µ = (p − l) µ acts as p µ , while h ν acts as p ν = Xp ν .
As before, the result is multiplied by the terms 1/l 2 and 1/l 2 from the t-channel gluon propagators and by the term 1/(q 2 − m 2 c ) from the corresponding heavy quark propagator. Then we have to subtract the part generated by the LO evolution equation which is given by the convolution P LO ⊗ C LO . Finally we integrate over l 2 ⊥ , accounting for the condition |l 2 | < Q 2 0 (the longitudinal components are fixed by the heavy quark and gluon (p − l) 2 = 0 on mass-shell conditions). In this way we obtain the power correction which should be subtracted from the known NLO gluon coefficient function Imf g (y) given by eq. (3.95) of [1] (see also [10]), which we then use to obtain the Q 0 subtracted NLO gluon contribution.