Open charm production and low x gluons

We compare the rapidity, $y$, and the beam energy, $\sqrt{s}$, behaviours of the cross section of the data for $D$ meson production in the forward direction that were measured by the LHCb collaboration. We describe the observed cross sections using NLO perturbative QCD, and choose the optimal factorization scale for the LO contribution which provides the resummation of the large double logarithms. We emphasize the inconsistency observed in the $y$ and $\sqrt{s}$ behaviours of the $D$ meson cross sections. The $y$ behaviour indicates a very {\it flat} $x$ dependence of the gluon PDF in the unexplored low $x$ region around $x\sim 10^{-5}$. However, to describe the $\sqrt{s}$ dependence of the data we need a steeper gluon PDF with decreasing $x$. Moreover, an even steeper behaviour is needed to provide an extrapolation which matches on to the well known gluons found in the global PDF analyses for $x\sim 10^{-3}$. The possible role of non-perturbative effects is briefly discussed.


Introduction
The LHCb collaboration have published measurements of the cross section for D meson [1,2,3] and B meson [4,5] production in the forward direction with rapidities in the region 2 < y < 4.5. D meson production is measured via the decays D → Kπ, Kππ, and the data are available, differential in both rapidity and transverse momentum, at three beam energies √ s = 5, 7 and 13 TeV. The B meson data, obtained via b → J/ψX decays, are available at 7 TeV differential in both y and p t [4], but the more recent measurements via semileptonic decays at 7 and 13 TeV are presented only inclusively in p t [5].
Here, we concentrate on D meson production, since we wish to probe low scale distributions which are more sensitive to the input gluon PDF. The data can be described by the production of a cc-pair followed by the fragmentation of the c quark into the D meson. Moreover, gluon fusion, gg, is the major contributor to forward D meson production. Therefore, since the mass of the c quark, m c , is not too high and that there is a large rapidity, this process allows a probe of the gluon distribution g(x, µ 2 F ) at very small x [6]− [10] x ∼ (m T / √ s)e −y ∼ 10 −5 (1) and at a small factorization scale µ F ∼ m T . Here m T = m 2 c + p 2 t,c , where p t,c is the transverse momentum of the c quark and √ s is the centre-of-mass proton-proton energy. Due to the absence of data probing the gluon in the low x domain the gluon PDF is practically undetermined in this region by the global PDF analyses. The uncertainty at a scale µ F = 2 GeV is illustrated in Fig. 1. The low x gluon density as given by recent NLO global parton analyses NNPDF3.0 [11], MMHT2014 [12], CT14 [13], calculated using the PDF interpolator LHAPDF [14].
Thus the open cc data are potentially very valuable. However a major problem is the sensitivity to the choice of factorization and renormalization scales. Note that in this low x region the gluon density strongly depends on the factorization scale µ F due to the presence of Double Log (DL) terms [(α s N C /π)ln(1/x)lnµ 2 F ] n with a large ln(1/x) ∼ 10. The conventional choice of scale is µ F = m T . Indeed, it was shown that the DL terms can be resummed into the incoming parton PDFs by choosing µ F 0.85m T [15]. Though the dependence on the residual scale, µ f , becomes weaker, it is nevertheless still appreciable.
It is known that the ratio of measured cross sections at different energies or at different rapidities is less sensitive to the value chosen for the factorization or renormalization scales [6]− [10]. However, ratios do not determine the overall normalization of the gluon density. On the other hand it would be valuable to fix the absolute value of the gluon density from these low x data. Even though the cc data have much lower statistical weight in comparison to the data used in the global PDF analyses, a reliable estimate of the gluon density in the unexplored region x ∼ 10 −5 , Q 2 < ∼ 10 GeV 2 would therefore be of great value. Among the papers that have studied the impact of cc and/or bb LHCb data on the behaviour of the gluon in the low x region, two papers have direct bearing on our study. First, [7] adds the cc and bb LHCb data to an ensemble of other data in a global PDF fit. The paper describes two types of fit. The most relevant is the fit where they attempt to determine the normalization of the low x gluon PDF. Since the result must be consistent with the gluons obtained from the whole set of data included in their global analysis, they find this can only be achieved by allowing the factorization and renormalization scales for both cc and bb production to be free parameters. They find that the data require the scales Am T to be We note that µ c R = 0.44m T was also found in [9]. Such unnaturally low scales indicate that actually the heavy quark data needs low x gluons larger than those extrapolated from the global PDF analyses [11]− [13]. In other words the value of the cross section, σ ∼ α 2 s (µ 2 R ) g(x 1 , µ 2 F ).., which was described by a low scale µ R , may be reproduced using a natural scale and a larger gluon, g(x 1 ).
The second paper, Gauld [10], concentrates on attempting to determine the low x gluon from the bb LHCb data [4,5], and compares with results he obtains from the cc data. He finds a large discrepancy with these forward production data at 13 and 7 TeV, both through fits to the normalized cross sections and to various cross section ratios. He concludes that the LHCb data are not consistent with evolution via perturbative QCD, assuming that there is no dramatic change in the input gluon behaviour in the region x ∼ 10 −4 − 10 −3 .
Here we proceed differently, although we, too, work at NLO. For the partonic subprocesses a(x 1 )b(x 2 ) → ccX, which drive the forward open cc production, we take the partons from the global PDF analyses with the exception of the gluon PDF at very low x ≡ x 1 say. Recall that the beam energies of the LHCb experiments sample, for low p t and large y, the x intervals Although we include the contributions from production by NLO gq and qq fusion, we note the dominance of production by gg fusion, so the assumption that only the low x gluon PDF is unknown is reasonable. We therefore perform fits to the LHCb open charm data using various parametrizations of the low x gluon density. We proceed in three stages, which may be summarised as follows.
(1) We first explore how well a simple two-parameter form for the low x gluon describes the √ s and y behaviour of the data in each p t bin independently, see Section 3. Since the gluon distribution at or near the input of DGLAP evolution is the most interesting 1 , here we consider the LHCb data in the four p t intervals from 1 to 5 GeV.
(2) Then, more ambitiously, we fit to all the LHCb data in the interval 1 < p t < 5 GeV simultaneously using a simple two-parameter Double Log (DL) parametrization (which is known to well approximate DGLAP evolution at low x throughout this p t range). We obtain a satisfactory description of all the open charm LHCb data, namely χ 2 = 141 for 120 D meson data points, see Section 4. So far, so good. Thus, up to the sizeable uncertainty due to the choice of scales we know the behaviour of the gluon PDF around x = 10 −5 . The problem is that when our fit is extrapolated from x = 10 −5 to x = 10 −3 , into the domain where the gluon is well known from the global PDF analyses we see a mismatch of the behaviour 2 .
(3) Indeed, the extrapolated gluon is more than a factor of two above the well known gluon PDF of the global analyses at x = 10 −3 . To attempt to cure this problem we re-fit the data using a lower scale than the optimal choice µ F (= µ R ) = 0.85m T . To be precise we take µ f (= µ R ) = 0.5m T and also we adjust the DL parametrization, see Section 5. As a result we achieve a good matching at x = 10 −3 for all Q 2 of interest. However, this fit is not satisfactory, as we shall explain below.

Description of the data
The LHCb collaboration have measured open cc production at three different beam energies 3 7, 13, 5 TeV [1,2,3]. The data, d 2 σ/dp t,D dy, are presented for five D meson rapidity intervals in the range 2 < y < 4.5 and we use four transverse momentum bins covering the range 1 < p t,D < 5 GeV. We fit to the data for D ± , D 0 ,D 0 meson production.
The c → D fragmentation functions D(z) was taken from [16], where they were determined from e + e − annihilation data in the Υ(4S) region. Actually the relative normalization of c quark fragmentation to the different channels is only known from previous data to about 10% accuracy. Therefore we allow, via a parameter N D , for an additional renormalization of the D 0 ,D 0 relative to the D ± data.
The FONLL programme [17] was used to calculate the open charm cross section at NLO. The running coupling, the charm mass (1.4 GeV), the quark and high x gluon distributions are given by those found by the MMHT2014 NLO global fit [12]. We now describe, in turn, the three methods mentioned above, to determine the low x gluon from the cc LHCb data. 1 At large scales the gluon PDF is formed mainly by evolution depending on parton PDFs with larger x at low scale. 2 Here the optimal scale µ F = µ R = 0.85m T , which resums the large DL terms, is used. 3 Note that the measurements at 13 and 5 TeV have recently been corrected. We use these updated data.

Gluons at fixed scales with a simple parametrization
The data in each of the four p t intervals (1 − 2, 2 − 3, 3 − 4, 4 − 5 GeV) were fitted separately assuming a simple two-parameter form for the low x behaviour of the gluon with x = x 1 and x 0 = 10 −5 . In this way we obtain gluons at four different scales, µ F . Note that the value of µ F is a bit larger than p 2 t,D + m 2 D since after fragmentation the transverse momentum of the c quark is given by In fact µ F = 2.0, 2.9, 3.9, 4.9 GeV for p t,D = 1.5, 2.5, 3.5, 4.5 GeV respectively.
As mentioned above, we allowed an extra normalization parameter, N D , between the D 0 ,D 0 and the D ± data. It was found to be N D 1.1 in every case. In detail, this means we take the fractions 0.246 and 1.1(0.565) for the D ± and (D 0 +D 0 ) charm quark decay channels respectively, leaving 0.133 for the D s + Λ c channels. The results are essentially unchanged if we set N D = 1, but then χ 2 is a bit larger.
The normalization N and the power λ in (4) and their uncertainties are obtained by fitting to the data in each p t,D interval using the MINUIT numerical minimization code [18,19]. The results of the four fits are presented in Table 1 and by the dashed curves in Fig. 2.  Table 1: The parameters N and λ giving the low x behaviour of the gluon distribution, xg(x, µ F ) = N (x/x 0 ) −λ from individual fits to the LHCb open charm data [1,2,3] in the four different p t,D intervals. The scale µ F and p t,D are given in GeV. The total χ 2 , χ 2 all , in each interval is shown, together with the contributions from the 5, 7 and 13 TeV data sets.
From the individual χ 2 contributions in Table 1 we see hints of tension between the 5+13 TeV data on the one hand and the 7 TeV on the other hand. Surprisingly if the 5+13 TeV data are fitted without the 7 TeV data (and vice-versa) we find similar values of N and λ, although the values of χ 2 are reduced.

Combined fit with a Double Log (DL) parametrization
As seen from the results of the simple fits shown in Table 1, the gluon density increases strongly with scale, while the power of the x behaviour has a weaker scale dependence. It is not evident whether such a behaviour is consistent with DGLAP evolution. Since in the low x region the major effect comes from the DL contribution, it is reasonable to attempt to fit all four groups of data simultaneously with the formula with With three light quarks (N f = 3) and N c = 3 we have β 0 = 9. The resummation of the leading double logarithmic terms (α s ln(1/x)ln(µ 2 )) n is written explicitly in the exponential, while the remaining single log terms are now parametrized by the powers a and b. We take, in (7), Λ QCD = 200 MeV and Q 0 = 1 GeV. As before, we set x 0 = 10 −5 .
Such an ansatz was used successfully to describe J/ψ and Υ photoproduction data [20]. Moreover, it was checked that in the x, µ 2 F region of interest this formula was consistent with NLO DGLAP evolution. So we fix the power b, which is responsible for the µ 2 F behaviour, to be the same as that found in the fit to J/ψ photoproduction. Now we are left to describe 120 LHCb data points with only two free parameters: N DL and a, plus the parameter N D introduced in Section 3.
The parameters of the combined fit are presented in the first row of Table 2, and in Fig. 2 the results are compared with those of the simple fits. The low x gluons obtained in the two fits are very similar. The DL description, which embodies NLO DGLAP evolution, should be more reliable than the results of the simple fits. Indeed, the dashed curves in Fig. 3 show that the 'combined' fit gives an acceptable description of all the D meson data in the interval 1 < p t < 5 GeV.
However, a problem is clearly evident. If the gluon PDF, determined from the LHCb data in the region of x ∼ 10 −5 , is extrapolated up to x ∼ 10 −3 then it greatly exceeds the well known gluon densities determined by the global PDF analyses. Can the gluon density determined from the cc data be reconciled with the global gluon PDF?

A fit which matches to the 'global' gluon at x = 10 −3
Here we attempt to find a parametrization of the gluon PDF that fits the LHCb cc data and which is consistent with the large x gluon which is well known from the global PDF analyses [11]− [13]. That is, we need to find a parametrization which reduces the extrapolated 'cc ' gluon by more than a factor of two at x ∼ 10 −3 for all Q 2 of interest. We therefore choose lower y dσ/dydp t D + +D -5 TeV  factorization and renormalization scales 4 µ f = µ R = 0.5m T . The main effect comes from 4 In the LO part of the contribution we still choose the 'optimal' scale µ F = 0.85m T which provides the resummation of the DL [(α s N C /π)ln(1/x)lnµ 2 F ] n terms, see [15]. There still remains a dependence on the choice of the residual factorization scale µ f .   [1,2,3] (120 data points in total) using the DL parametrization given by eqs. (6) and (7) of Section 4. We also show the contributions of the three data sets to the total χ 2 = 141. For comparison we show the parameters of the gluon obtained in a similar fit [20] to J/ψ data. Note, however, that the gluons obtained in [20] from J/ψ data are not the MS gluons but correspond to those of the physical factorization scheme. The last row gives the values of the parameters of the fit described in Section 5; the much smaller value of N DL is compensated by a larger argument of the exponential factor in (6) (due to Q 0 reducing from 1 to 0.5 GeV.) the change of renormalization scale, since the LO cross section is proportional to α 2 s (µ 2 R ). The gluon density does decrease, but it still strongly exceeds the global density at x = 10 −3 . To overcome this problem we reduce the value of the parameter Q 0 in the DL form (6). We take Q 0 = 0.5 GeV. In this way we reach a satisfactory matching 5 at x = 10 −3 , at the price of considerably worsening the fit to the data, as can be seen by the solid lines in Fig. 3 and the χ 2 values in the last row of Table 2. Is this a problem of the fit to the data or it is a problem of the data themselves?
It appears that the rapidity dependence of the cc cross section disagrees with the energy, √ s, dependence of the data. Indeed, Fig. 3 shows that in order to describe the rapidity behaviour of the data, we need a curve that decrease faster with rapidity -that is, a low x gluon density which should increase slower with decreasing x in comparison with that in the fit. On the contrary, to obtain a larger σ(13 TeV)/σ(5 TeV) cross section ratio we need a low x gluon which increases faster with decreasing x.
Another view of the same problem can be observed if we compare the data at larger beam energy √ s 1 and large rapidity with data at a smaller beam energy √ s 2 and a smaller rapidity shifted by ln(s 1 /s 2 ). In this case we deal with the same value of large x 2 , so the ratio of cross sections should be equal to the ratio of the gluons at x 1 and x 1 (s 1 /s 2 ). The ratio of these cross sections is close to 1 for D 0 , and is about 1.1 − 1.2 for D + , indicating that λ 0 − 0.1 in terms of (4), in agreement with Table 1. This implies, that in this interval of x, the gluons are almost constant. However, in order to match these 'flat' small x gluons with the 'global' gluons at x = 10 −3 , the 'cc ' gluons have to decrease rapidly (by a factor of about 4) in the interval from x = 10 −5 to x = 10 −3 ; which corresponds to λ 0.3. It is very unnatural to have such a rapid qualitative change in the x behaviour of the gluon density. No reasonable physical effect can generate this behaviour.
Note also, that contrary to cc production, the same analysis of the ratio of B meson cross sections indicates the behaviour of the low x gluon is xg ∝ x −λ with an unexpectedly large λ > 0.6 for x > 10 −4 [10].

Discussion
Note that it looks quite unnatural to describe data with such a low renormalization scale as that in Section 5. First, the idea of the collinear factorization approach was that the infrared behaviour of the amplitude is described by the incoming parton PDF, while the short range interaction is collected in the hard matrix element. A low µ R means that we include the running of the QCD coupling up to a large distance. Increasing the distance (as, for example, in the scales given in (2) for charm) means that we enter the region already described via the PDFs if µ R < µ F . Technically if µ R < µ F we have double counting of the loop insertion in the gluon propagator. It is included in the DGLAP evolution of the PDF up to µ F and the same loop is included in the renormalisation 6 of α s . That is why we take µ R = µ F , and, in Section 5, Thus we come back to the problem of matching at larger x with the gluons from the global analyses. To provide such a matching we need, in the interval of x = 0.5 · (10 −4 − 10 −3 ), a larger value of λ > 0.3 − 0.4. In general, the decrease of λ, with decreasing x, is expected. This is caused both, by stronger absorptive corrections at smaller x and by the DL effects. The problem is that, already at x = 10 −3 , the 'global gluons' (gluons from global analyses) are already rather flat (see Fig. 2, especially at a lowest scale (p t = 1.5 GeV). For example, for µ F = 2.0 (2.9) GeV the 'global gluons' correspond to λ < ∼ 0.17 (0.23) at x = 10 −3 and there is no reason for λ to increase at x < 10 −3 .

Alternative parametrization
In the present paper we have used the DGLAP form of the gluon distribution in order to compare the low x gluon PDF obtained from the analysis of the LHCb open charm data, with that given by the parton global analyses. On the other hand, at such small values of x the QCD dynamics may be better described by the BFKL equation. The solution of BFKL equation leads to a power-like x −λ behaviour (analogous to used in eq. (4) above) with a prefactor which weakly depends on x. That is, the analysis in Sect. 3 can be considered as a BFKL-based description of the data. Since the analysis was done for each p t bin individually, we do not have to worry about the form of the scale (or p t ) dependence of the gluon parametrization.
If we were to include the prefactor, and write the asymptotic BFKL amplitude as then the matching with the high x global gluons would be even worse. Indeed, the effective power would decrease with increasing x and the extrapolation of our solution to larger x would give even larger gluons than those shown in Fig.2.
However, the low value of λ ∼ 0.01−0.1 indicates that in our relatively low-p t domain we do not deal with a pure BFKL amplitude (that is a single QCD Pomeron). The expected BFKL intercept should be larger. At NLL level the value of λ = 0.2 − 0.3 and weakly depends on the scale, see, for example [22,23].
It might be thought that the low values of λ, that we obtain in Table 1, may be explained by a possible saturation effect; since, as x decreases, the growth of absorptive corrections damps down the rise of the low-x gluons. As is seen from Fig.2, the extrapolation of the conventional 'global gluons' leads, at x ∼ 10 −5 , to gluon densities which are consistent with our gluons (needed to describe the LHCb open charm data), and if then the value λ will be strongly diminished due to the saturation effect, this could provide reasonable matching with the larger x 'global gluons'. However, this explanation does not work in practice.
The problem is the p t dependence. For a smaller value of p t the cross section of an additional interaction is larger, and consequentially the absorptive corrections are stronger. So, for example, in the p t = 4.5 GeV bin (µ F = 4.9 GeV) we expect a much weaker effect than that for p t = 1.5 GeV (µ F = 2 GeV). Moreover, for larger p t , we sample larger x. On the other hand, we see from Fig. 2 that at, say, x = 10 −4 and µ F = 3.9 or 4.9 GeV, the density of 'global gluons' is already about twice smaller than that from our analysis (needed to describe the data) and including the saturation effect we will only diminish it further, and make the discrepancy larger. That is, while there is a chance to get more or less a satisfactory description of the lowest p t bin, by tuning the parameters of the absorptive (saturation) effect, we will surely fail to describe the data at larger values of p t .
Besides this, we have to recall, that no saturation effects were observed in the analysis of data for exclusive J/ψ production at the LHC. These data, which probe more or less the same kinematic region 7 , are well described by the DL parametrization of eq. 6 (see e.g. [20]).

Consistency of heavy quark data
As discussed in [10], the disagreement between the experimentally observed charm/bottom cross sections and the QCD predictions may be caused by a deficiency in calculating the detector efficiencies in the different y intervals. Recall that the experimental rapidity behaviour for the D-meson cross section leads to a very flat low x gluon PDF, xg, see the low value of λ in Table 1. When this gluon is extrapolated to larger x we have a large mismatch with the well known gluons of the global PDF analyses. On the other hand, to describe the B-meson rapidity distribution [4] (which is rather flat in y) we need much steeper gluons. Fig. 4 shows an example of comparing D and B meson production. The figure compares the QCD predictions calculated with different sets of gluon PDFs for the 7 TeV data taken at p t = 4.5 GeV for D-meson production (the left plot) with a corresponding plot of B-meson data at p t = 1.25 GeV. These two plots correspond to practically the same scale (the same value of m T ) and the same x interval. We see that the gluons found in the present paper, which describe the D-meson cross sections, strongly overestimate B-meson production. Moreover, the B-meson data ask for even steeper gluons than those given by the central values of the gluon PDFs of the CT14 and MMHT2014 global parton analyses. Indeed, to reproduce such a flat rapidity dependence of the B-meson cross section shown in Fig. 4, the gluons should increase faster with decreasing x.
Note also that even using the low renormalization scale µ R = 0.5m T , and the correspondingly smaller gluons which at x=0.001 match those from the global analyses, we overestimate the D-meson cross section measured at √ s = 1.96 TeV by the CDF collaboration [24] in just the x ∼ 0.001 region. For p t = 2 GeV we find 53.8 µb/GeV as compared to the CDF result of 32.7 ± 6.5 ± 4.2 µb/GeV. The predicted cross section is too large -not due to the new (larger) gluon PDFs, but due to the large value of α s (µ R ) taken at the low scale µ R = 0.5m T .
All these facts indicate an inconsistency in the data.
Another possibility is that at these relatively low scales (that is, large distances) relevant for charm production we feel non-perturbative effects. Note that a large discrepancy is observed in Fig. 3 for the smallest p t bin. There could be the non-perturbative production of a cc pair or non-perturbative effects in the c → D fragmentation. Unlike the e + e − case, for inclusive pp processes the c quark is surrounded by a large number of light quarks from the underlying event. Then, besides fragmentation, the D meson may be formed by the fusion of a charm and a light quark. To investigate the latter case it would be interesting to see the dependence of D meson production in events with different multiplicities.
However a similar problem is observed for B meson production (see (2) and [7,10]) where we deal with a larger scale and so the non-perturbative effects will be heavily suppressed. Therefore, in order to isolate the role of non-perturbative effects, it would be very interesting to see more precise and more differential data, especially for low p t B meson production.
Recall that in spite of the large mass m b , there is the possibility to probe low factorization scales (µ F < m b ) by observing both the B andB mesons and selecting the events with small acoplanarity, see the discussion in Section 5 of [15].