Scale dependence of open $c\bar{c}$ and $b\bar{b}$ production in the low $x$ region

The `optimal' factorization scale $\mu_0$ is calculated for open heavy quark production. We find that the optimal value is $\mu_F=\mu_0\simeq 0.85\sqrt{p^2_T+m_Q^2} $; a choice which allows us to resum the double-logarithmic, $(\alpha_s\ln\mu^2_F\ln(1/x))^n$ corrections (enhanced at LHC energies by large values of $\ln(1/x)$) and to move them into the incoming parton distributions, PDF$(x,\mu_0^2)$. Besides this result for the single inclusive cross section (corresponding to an observed heavy quark of transverse momentum $p_T$), we also determined the scale for processes where the acoplanarity can be measured; that is, events where the azimuthal angle between the quark and the antiquark may be determined experimentally. Moreover, we discuss the important role played by the $2\to 2$ subprocesses, $gg\to Q\bar{Q}$ at NLO and higher orders. In summary, we achieve a better stability of the QCD calculations, so that the data on $c\bar{c}$ and $b\bar{b}$ production can be used to further constrain the gluons in the small $x$, relatively low scale, domain, where the uncertainties of the global analyses are large at present.


Introduction
The present global PDF analyses (e.g. NNPDF3.0 [1], MMHT2014 [2], CT14 [3]) find that there is a large uncertainty in the low x behaviour of the gluon distribution. There is a lack of appropriate very low x data, particularly at low scales. However, recently measurements on open charm and open beauty in the forward direction have been presented by the LHCb collaboration [4,5,6,7]; moreover, the ATLAS collaboration has measured open charm production in the central rapidity region [8]. These data sample the gluon distribution at rather low x: namely in the domain 10 −5 < ∼ x < ∼ 10 −4 . A discussion of the data in terms of existing global PDFs has been presented in [9,10], and they have been incorporated in a fit with the HERA deep inelastic data in [11].
In the ideal case it would be good to have such data where both the heavy quark and the heavy antiquark were measured, since when we observe only one quark (one heavy hadron) the value of x that is probed is smeared out over an order of magnitude by the unknown momentum of the unobserved quark in the QQ-pair [12,9], where Q ≡ c, b, see e.g. Fig.1 in [9]. Nevertheless, even measurements of the inclusive cross section of one heavy quark can be used to check and further constrain the existing PDFs.
Another problem, which was emphasized in [10], is that the QCD prediction at NLO level strongly depends on the factorization scale, µ F , assumed in the calculation. We might expect that the major source of the strong µ F dependence arises because in the DGLAP evolution of low x PDFs the probability of emitting a new gluon is strongly enhanced by the large value of ln(1/x). Indeed, the mean number of gluons in the interval ∆ ln µ 2 F is [13] n α s N C π ln(1/x) ∆ ln µ 2 F , leading to a value of n up to about 8, for the case ln(1/x) ∼ 8 with the usual µ F scale variation interval from µ F /2 to 2µ F . In contrast, the NLO coefficient function allows for the emission of only one gluon. Therefore we cannot expect compensation between the contributions coming from the PDF and the coefficient function as we vary the scale µ F . It was shown in [14,15] that this strong double-logarithmic part of the scale dependence can be successfully resummed by choosing an appropriate scale, µ 0 , in the PDF convoluted with the LO hard matrix element, which in our case is M(gg → QQ).
The outline of the paper is as follows. In Section 2 we recall the method of performing the resummation to determine the optimal scale µ 0 . In Section 3 we justify choosing the renormalization scale equal to the factorization scale. Then in Section 4.1 we use the procedure discussed in Section 2, to resum the ln(1/x) terms so as to determine the optimum factorization scale, µ 0 . Unfortunately for heavy QQ production (unlike the Drell-Yan process) a large sensitivity to the choice of scale remains. In Section 4.2 we identify the source of the problem to be the important 2 → 2 (that is gg → QQ) diagrams at NLO and higher orders. We argue that it is possible to also resum these diagrams. We then find the scale sensitivity is reduced. It would be advantageous if both heavy mesons (arising from Q andQ) could be measured experimentally, but, at present, the statistics are limited. However, a possibility to circumvent this problem is discussed in Section 5. In Section 6 we return to open single inclusive cc and bb production and compare the QCD predictions with the optimal scale with LHC data; and are able to make an observation about the gluon PDF at low x. In Section 7 we present our conclusion.

Way to choose the optimum factorization scale
Here we recall the procedure proposed in [14,15], which provides a reduction in the sensitivity to the choice of factorization scale by resumming the enhanced double-logarithmic contributions from a knowledge of the NLO contribution. The cross section for open heavy quark production at LO + NLO at factorization scale µ f may be expressed in the form 1 where the coefficient function C (0) does not depend on the factorisation scale, while the µ f dependence of the NLO coefficient function arises since we have to subtract from the NLO diagrams the part already generated by LO evolution.
We are free to evaluate the LO contribution at a different scale µ F , since the resulting effect can be compensated by changes in the NLO coefficient function, which then also becomes dependent on µ F . In this way eq. (2) becomes is calculated now at the scale µ F used for the LO term, and not at the scale µ f corresponding to the cross section on the left hand side of the formula. Since it is the correction which remains after the factorization scale in the LO part is fixed, we denote it C (1) rem (µ F ). Note that although the first and second terms on the right hand side depend on µ F , their sum, however, does not (to O(α 4 s )), and is equal to the full LO+NLO cross section calculated at the factorization scale µ f .
Originally the NLO coefficient functions C (1) are calculated from Feynman diagrams which are independent of the factorization scale. How does the µ F dependence of C (1) rem in (3) actually arise? It occurs because we must subtract from C (1) the α s term which was already included in the LO contribution. Since the LO contribution was calculated up to some scale µ F the value of C (1) after the subtraction depends on the value µ F chosen for the LO component. The change of scale of the LO contribution from µ f to µ F also means we have had to change the factorisation scale which enters the coefficient function C (1) from µ f to µ F . The effect of this scale change is driven by the LO DGLAP evolution, which is given by where P left and P right denote DGLAP splitting functions acting on the PDFs to the left and right respectively. That is, by choosing to evaluate σ (0) at scale µ F we have moved the part of the NLO (i.e. α s ) corrections given by the last term of (4) from the NLO to the LO part of the cross section. In this way C (1) becomes the remaining µ F -dependent coefficient function C does not contain the double-logarithmic (α s ln(µ F )ln(1/x)) n contributions. It is impossible to nullify the whole NLO contribution since the function C (1) (µ) depends also on other variables; in particular, it depends on the mass,ŝ, of the system produced by the hard matrix element. On the other hand we can choose such a value of µ which makes C (1) (µ,ŝ) = 0 in the limit of largê s m 2 Q . Recall that the ln(1/x) factor arises in the NLO after the convolution of the largê s asymptotics of hard subprocess cross section with the incoming parton low-x distributions satisfying At NLO level the change µ f to µ F is irrelevant: eq.(3) is an identity (it just changes the higher order terms). However, in this way we simultaneously resum all the higher order doublelogarithmic contributions in the PDFs(µ F ) of the LO part. As a result we are able to suppress the scale dependence caused by large values of log(1/x).
Thus the choice of µ = µ 0 , which nullifies C (1) at largeŝ m 2 Q , excludes the Double Log (DL), α s ln µ 2 F ln(1/x), contribution from the NLO correction by resumming the series of double-logarithmic terms in the PDFs, which are then convoluted with the LO coefficient functions. To find the appropriate value of µ 0 we must choose the NLO subprocess driven by the same ladder-type diagrams (in the axial gauge) as the ladder diagrams that describe LO DGLAP evolution. The appropriate subprocess is gluon-light quark fusion, gq → QQq. In the high energy limit, where the subprocess energy satisfiesŝ(gq) m 2 Q , the cross section described by this subprocess contains double-logarithmic terms log(µ 2 F /µ 2 0 ) log(ŝ/m 2 Q ). The subprocess gq → QQq is contained in the sketch of Fig. 1(b), where it is shown pictorially how the enhanced double-logarithmic terms are transferred to the PDFs in the LO term.

Extension to higher orders
We note that, in general, this decomposition can be continued to higher order. For example, if the NNLO contribution is known, then we will have three scales: µ f , µ F = µ 0 and µ 1 , (6) where the scale µ 1 is chosen to nullify the final term in the small x limit.
In fact in Section 4.2 we will use this equation to include the important 2 → 2 (that is, gg → QQ) subprocess at NLO and higher orders. We will show reasons why the scale choice µ 1 = µ 0 will give a good approximation for the resummation of these higher-order 2 → 2 contributions.

Comparison with k t factorization
The approach we have introduced is based on collinear factorization. However, actually it is close in spirit to the k t -factorization method. Indeed, there, the value of the factorization scale is driven by the structure of the integral over k t , see Fig. 2. In the k t -factorization approach this k t integral is written explicitly, while the parton distribution unintegrated over k t is generated by the last step of the DGLAP evolution, similar to the prescription proposed in Refs. [16,17]. Then, using the known NLO result, we account for the exact k t integration in the last cell adjacent to the LO hard matrix element. This hard matrix element M, provides the convergence of the integral at large k t . In this way it puts an effective upper limit of the k t integral, which plays the role of an appropriate factorization scale.

The renormalization scale µ R
Besides the factorization scale, the QCD prediction, truncated at NLO, strongly depends on the renormalization scale µ R , since the LO term is already proportional to α 2 s (µ R ). Let us discuss the possible choice of µ R . First, it is reasonable to have µ R > ∼ µ F , since we expect all the contributions with virtualities less than µ F to be included in the PDFs, while those larger Q q k t Figure 2: The diagram for AA * , where A is the amplitude for the subprocess gg → QQq shown in Fig. 1(b). However, in the k t factorization approach k t is integrated over and the effective upper limit of the convergent integral essentially plays the role of the appropriate factorization scale.
than µ F to be assigned to the hard matrix element. This is in line with the fact that the current scale of the QCD coupling increases monotonically during the DGLAP evolution. So the coupling responsible for heavy quark production should have a scale µ R equal to, or larger than, that in the evolution.
Another argument is based on the BLM prescription [18], which says that all the contributions proportional to β 0 = 11 − 2 3 n f should be assigned to α s by choosing an appropriate scale µ R . A good way to trace the β 0 contribution is to calculate the LO term generated by a new type of light quark, so n f → n f + 1. Note that the new quark-loop insertion appears twice in the calculation. The part with scales µ < µ F is generated by the virtual (∝ δ(1 − z)) component of the LO splitting during DGLAP evolution, while the part with µ > µ R accounts for the running α s behaviour obtained after the regularization of the ultraviolet divergence. In order not to miss some contribution and to avoid double counting we take the renormalization scale µ R = µ F . The argument for this choice was made in more detail in [19] for the QED case.
Of course, all these are only the arguments why we expect µ R = µ F , and are not a proof. Formally we can only say that we expect µ R to be of the order of µ F . Thus there could be further uncertainty in the scale dependence of the predictions due to the possibility that µ R = µ F . However, based on these arguments, below we study the factorization scale dependence using the renormalization scale µ R = µ 0 .
We emphasize (see also [10]) that the renormalisation scale dependence affects just the normalization of cross section, but not its energy behaviour. It is cancelled in the ratio of the cross sections measured at the LHC energy of 7 (or 8 or 5) TeV to that at 13 TeV or in the ratio of the cross sections obtained at different rapidities. Thus these ratios will probe the low x dependence of the gluons at scale µ F = µ 0 essentially without any uncertainties due to possible variations of the µ R scale.

Sensitivity of predictions to the factorization scale
Here we implement the proposals of eqs. (3) and (6) in an attempt to reduce the factorization scale dependence of the QCD predictions for QQ production in high-energy pp collisions. Note, however, that calculating the NLO contribution of the diagram in Fig. 1, we have integrated over the momenta of other particles; in particular, over the transverse momentum, −k t , of the light quark.

The optimum scale to resum (α s ln(1/x)lnµ 2 ) n terms
We use the formulae from appendix B of the [20] paper to calculate the gg → QQq matrix element in the high energy limit in order to find a scale, that nullifies the double-logarithmic NLO contribution: that is, to find a scale µ 0 at which the DGLAP-induced contribution (P left ⊗ C (0) + C (0) ⊗ P right ) replaces the NLO correction calculated explicitly. Note, however, that calculating the NLO contribution of the diagram in Fig. 1, we have integrated over the momenta of other particles; in particular, over the transverse momentum, −k t , of the light quark. Since we are going to consider the upper (heavy quark) box in Fig. 1 as the 'hard' subprocess, and would like to keep the DGLAP k t ordering, we put an additional cut −|k t | < min{m T Q , m T Q }; otherwise the lower part of diagram (which may be either qQ or qQ scattering) may have k t > m T , and would then be treated as the hard subprocess. Here m T = m 2 Q + p 2 T . The values that we find for the 'optimal' scale µ 0 are presented in Fig. 3 as the function of p T /m Q ratio, where p T is the transverse momentum of the observed heavy quark. It turns out that the values of the optimal scale are close to the value µ 2 F = m 2 T ≡ p 2 T + m 2 Q that is used conventionally; that is F = 1. However, we now have a physics justification for the scale choice shown in Fig. 3, which to a good approximation is µ 0 0.85m T , that is F 0.85. Now that we have the value of µ 0 , we can study the factorization scale, µ f , dependence of the QCD predictions for cc and bb production. The results are shown in Fig. 4 for a Q quark of pseudorapidity η = 3 -typical of the LHCb experiment. The curves for the first two procedures 2 , mentioned in the caption of Fig. 4, are obtained from (3) setting µ F (and µ R ) equal to µ 0 , and then varying µ f in the range (m T /2, 2m T ). We use the CT14 [3] PDFs as an example of a recent set of partons which have no negative gluon distributions and take the corresponding heavy quark masses: m c = 1.3 GeV and m b = 4.75 GeV. Subroutines from the MCFM [21] and FONLL [22] programmes were used for the computations. shown in Fig. 4 by the dashed red curves. We repeat the cross section prediction, but now use (2) with the conventional choice µ f = (0.5, 1, 2)m T , which gives the blue curves. Not surprisingly, since the optimum scale is close to the conventional choice µ 0 = m T , the scale uncertainties are comparable.
Unfortunately we still have rather strong dependence of the predicted cross section on the choice of the value of µ f . It is caused by the relatively low mass contribution coming mainly from the 2 → 2 (gg → QQ) component of C (1) . This component does not contain a ln(1/x) dependence, but, at our low scales, it is numerically large; it gives up to twice as large a contribution as the LO one. Moreover, being convoluted with low-x PDFs, which strongly depend on µ f , it produces a large scale uncertainty.

The optimum scale to resum the higher-order 2 → 2 diagrams
In order to reduce the scale dependence it turns out to be important to fix the scale of the 2 → 2 NLO contribution, that is to know the value of µ 1 in the 2 → 2 part of the second term on the right hand side of eq.(6). Strictly speaking to do this we have to know the NNLO expression. At the moment there exist only numerical NNLO result for t-quark pair production, see [23] and references therein. Nevertheless we can extract some use from these calculations. As was demonstrated in [24] (see Fig. 6 for example) the corrections to the 2 → 2 NLO contributions are mainly of Sudakov origin 3 ; -these are 'soft' corrections corresponding to a relatively small momentum transferred along an additional gluon. Such corrections does not change essentially the original kinematics of the 2 → 2 subprocesses or the dependence of the corresponding 'hard' matrix element on the virtuality of incoming parton. Thus it looks reasonable to convolute these terms with the same PDFs as those used for the LO evaluation. This will provide the correct resummation of the higher-order DL terms, (α s ln µ F ln(1/x)) n (with n = 2, 3, ...) inside the incoming parton distributions. Referring to (6), it means that we may argue that the 2 → 2 part of the NLO coefficient function C (1) must be convoluted with partons taken at the scale µ 1 = µ 0 . In other words we write the cross section as with µ F = µ 0 and α s (µ R ) = α s (µ 0 ), where we have divided the C (1) correction into two terms (2→3) , with only the second term evaluated at the residual factorization scale µ f . The corresponding results, calculated from (8), are shown in Fig. 4 by the dotted red curves. We see that the remaining µ f dependence is much reduced. We consider this observation as a strong argument in favour of the possibility of using open charm or beauty data to constrain the low x partons at the scale µ f = 0.85m T .
Since the major contribution to cc and bb production comes from gluon-gluon fusion (see Fig. 5), including these data in global parton analyses will allow a better study of the gluon low-x behaviour, and hence to strongly diminish the present uncertainty observed in this region.

Azimuthal cut to reduce optimal scale for bb events
In the case of open bb production the optimal scale is rather large; typically µ 2 0 > 30 GeV 2 . On the other hand, the main uncertainties in the gluon PDF are observed at much lower scales ∼ 2 -4 GeV 2 . One possibility to reduce the scale at which the process probes the partons is to observe both heavy quarks (i.e. both the quark and the antiquark), and then to select the events where the transverse momentum of the pair is small. This proposal was discussed in [12] (and in [14] for Drell-Yan pair production). Unfortunately, the transverse momenta of B mesons can only be measured for a few particular decay modes, and the product of the branching ratios for the two B mesons is small. It means that we do not have sufficient statistics.
Another idea was proposed by Alexey Dzyuba 4 . As a rule the vertex of B meson decay can be observed experimentally, and it is possible to measure the azimuthal angle, φ, between the two heavy mesons. That is, we may select BB events with good coplanarity. In such a case the transverse momenta of the incoming partons must be small, otherwise the coplanarity will  Table 1: The optimal factorization scale, µ 0 , corresponding to a coplanarity cut ∆φ < φ 0 for events with both heavy quarks in the rapidity interval 2 < y < 4.5. The cross section is integrated over the heavy quark transverse momenta p T .
be destroyed. In other words, for events with a small ∆φ = π − φ we deal with lower scale partons. For example, in Table 1 we show the optimal scale µ 0 (∆φ) calculated for events with ∆φ < φ 0 corresponding to the LHCb rapidity interval 2 < y < 4.5. As expected, for low φ 0 we have µ 0 ∝ φ 0 . For instance, for bb production with ∆φ < 10 o one can probe gluons at a rather low scale, namely µ 1.5 GeV.
6 Comparison with cc and bb data Now that we have the optimal factorization scale, µ 0 0.85m T , we can make an exploratory comparison with the existing LHC data for open single inclusive heavy-flavour production. To compare with the data we use the subroutines from MCFM and FONLL programmes [21,22]. The QCD description of the present data in the low-x, low-µ domain is shown in Fig. 6. As an example, we consider just the D + (B + ) meson cross sections using the probabilities of the quark to meson transition P (c → D + ) = 0.25 (see, for example, [25] p.208) and P (b → B + ) = 0.4 (see [26] and [27] p.63). We account for the fact that the D/B meson momentum is less than that of the parent quark by making the assumption that p D ∼ 0.75p c and p B ∼ 0.9p b (see [27,28]). That is, we calculate the meson cross sections as where the last factor (1/0.75 or 1/0.9) accounts for the ratio of the dp T,D(B) and dp T,c(b) intervals.
It is seen that the QCD predictions obtained using the 'optimal' factorization scale and the central values of CT14 NLO partons underestimate the LHCb cc data. Note, however, there are large uncertainties in the behaviour of the low-x gluon distributions obtained from the global parton analyses. This uncertainty may be reduced for the NLO partons 5 by including the open charm/beauty data in the global analysis and using the 'optimal' scale to calculate the corresponding cross sections. Of course, there are also the uncertainties due to higher α s order  Figure 6: The QCD predictions for the cross section for heavy meson (D + , B + ) production compared with LHCb data [5,7], as a function of the p T of the heavy meson. In the upper plot the blue (red) curves and data points, taken at √ s=13 TeV, correspond to the D + rapidity bins 2 < y < 2.5 (4 < y < 4.5) respectively; whereas the lower plot corresponds to B + rapidity in the interval 4 < y < 4.5 for a collider energy of √ s=7 TeV. CT14 NLO PDFs [3] are used. The optimal factorization scale is taken µ F = µ 0 = µ 1 = 0.85m T ; and the renormalization scale is taken to be µ R = µ F , see Section 3.
contributions not included into the calculations. When the NNLO formulae become available it will be possible to extend our procedure and to include open charm data into the NNLO global parton analyses.

Conclusion
We have calculated the 'optimal' factorization scale, µ 0 , which allows a resummation of the higher-order α s corrections, enhanced at high energies by the large ln(1/x) factor; that is to resum the double logarithmic, (α s ln µ 2 F ln(1/x)) n , terms and move them into the incoming parton distributions. The result is given in Fig. 3. It is essentially for single open inclusive heavy quark production, where p T is the transverse momentum of the observed heavy quark.
We also considered the case when the azimuthal angle, φ, between the heavy quark and the antiquark can be measured. We showed that by selecting events with small ∆φ = π − φ we are able to probe smaller factorization scales µ 0 . This is an advantage for bb production: compare the results of Table 1 with eq, (11). The disadvantage is that the rate is smaller for such events, even though we do not require that the transverse momentum of both the heavy quarks are measured.
The choice µ F = µ 0 reduces the uncertainty of the perturbative QCD calculations. It will allow LHC data on cc and bb production to be included in global parton analyses to constrain the behaviour of the gluon distribution in the region of very small x and low scale, equal to µ 0 , where the uncertainties of the present global parton analyses are especially large.