Understanding forward $B$ hadron production

The LHCb collaboration has recently performed a measurement of the production rate of inclusive $B$ hadron production ($pp\to BX$) at both 7 and 13~TeV centre-of-mass (CoM) energies. As part of this measurement, the ratio of these two cross section measurements has been presented differentially in $B$ hadron pseudorapidity within the range of $\eta_B \in [2.0,5.0]$. A large tension ($4\sigma$) is observed for the ratio measurement in the lower pseudorapidity range of $\eta_B \in [2.0,3.0]$, where the data is observed to exceed theoretical predictions, while consistency is found at larger $\eta_B$ values. This behaviour is not expected within perturbative QCD, and can only be achieved by introducing ad-hoc features into the structure of the non-perturbative gluon PDF within the region of $x\in[10^{-3},10^{-4}]$. Specifically, the gluon PDF must grow extremely quickly with decreasing $x$ within this kinematic range, closely followed by a period of decelerated growth. However, such behaviour is highly disfavoured by global fits to proton structure. Further studies of the available LHCb $B$ and $D$ hadron cross section data, available for a range of CoM energies, indicate systematic tension in the (pseudo)rapidity region of $[2.0,2.5]$.


Introduction
The LHCb collaboration has recently presented measurements of inclusive B hadron production in pp collisions at 13 and 7 TeV centre-of-mass (CoM) energies [1], defined through the process pp → BX. The cross section measurements are reported differentially with respect to B hadron pseudorapidity (η B ) within the range η B ∈ [2.0, 5.0], and inclusively with respect to transverse momentum (P B T ). In addition, the ratio of the differential cross section measurements at these two CoM energies has also been presented.
The motivation for considering the ratio of heavy quark cross section measurements is that many sources of (otherwise overwhelming) theoretical and experimental uncertainty, which are highly correlated at different CoM values, partially cancel in the ratio. At the same time, the ratio is still sensitive to the shape of the gluon parton distribution function (PDF) at both small and large values of Björken-x (x) [2,3], since typically different values of x are probed within a fixed kinematic region at different CoM values. Consequently, it is possible to include the heavy quark data at the level of the ratio into a global analyses of proton structure, improving the description of the gluon PDF. This method was recently applied [4] to the double differential D hadron ratio data provided by LHCb [5][6][7].
It therefore comes as a quite a surprise that significant tension is observed for the B hadron ratio data with respect to the corresponding theoretical predictions. In particular, the data is observed to exceed (≈ 4σ) the predictions in the range of η B ∈ [2.0, 3.0], while agreement is found for the more forward region of η B ∈ [3.0, 5.0]. This behaviour is unexpected for the following reasons.
• Firstly, while B and D hadron predictions typically probe different values of x and Q 2 of the input PDFs, there are kinematic regions where the two predictions are highly correlated. No tension is observed for the most precise (13/5 TeV) D hadron ratio measurement in these regions [5].
• Secondly, a striking feature of the B hadron data is that the ratio is observed to decrease with increasing B hadron pseudorapidity, which would indicate the presence of a region of accelerated then decelerated growth of the gluon PDF at values of x ∈ [10 −3 , 10 −4 ] and Q 2 ∼ 50 GeV 2 . This is not a feature of DGLAP evolution, so such a structure would have to be present in the non-perturbative gluon PDF. However, measurements of the heavy quark (charm and beauty) structure functions F qq 2 (x, Q 2 ) at HERA [8] do not find such a feature in this x range, where this sort of effect should be more pronounced since the relevant data is at lower Q 2 values.
The purpose of this work is to perform detailed studies of the available forward B hadron production data to better understand the possible origin of the observed deviation. The remainder of this paper is organised as follows. In § 2, the theoretical set-up for providing B hadron production predictions are discussed, and the kinematics relevant for B hadron production within the LHCb acceptance are studied. In § 3, the available LHCb B hadron cross section data is studied at the level of both the absolute and normalised cross sections at both 7 [9] and 13 TeV [1] CoM energies. After studying the cross section data, the ratio of the 13 and 7 TeV cross section measurements is studied in § 4. In addition to studying the differential ratio as measured by LHCb, a kinematically 'shifted' ratio is introduced which provides direct sensitivity to the growth of the low-x gluon PDF. In § 5, both the theoretical and experimental consistency of the LHCb B hadron ratio data is considered. Firstly, the theoretical consistency of the data is considered by comparing the experimentally extracted values for the growth of the gluon PDF with those obtained with a toy model for PDFs. Secondly, correlations between the predictions for B and D hadron are also considered, and the consistency of the D hadron ratio measurements are also discussed. Finally, some general discussion and conclusions are provided in § 6.
2 Theoretical set-up for forward B hadron production At the LHC, inclusive B hadron production is dominated by the gluon-fusion heavy quark pair production subprocess, and the predictions of the distributions of B hadrons can be obtained by convoluting the partonic cross section for heavy quark pair production with input PDFs and the relevant heavy quark fragmentation functions. The basis for the current state-of-the-art for differential cross section predictions is the next-to-leading order (NLO) partonic cross section [10][11][12][13][14], where predictions can be further improved by matching this massive calculation to a parton shower or a massless calculation. In the following, the theoretical set-up for providing B hadron predictions will be provided. In addition, the partonic kinematics relevant for forward B hadron measurements in the LHCb acceptance are also discussed. While the discussion is focussed towards B hadron production, the predictions for D hadron production proceed in essentially the same way.

General considerations
In the current studies, predictions are provided at NLO accuracy matched to a parton shower (NLO+PS), which is achieved with the POWHEG method [15][16][17] to match the heavy quark pair fixed-order calculation [18] with Pythia8 [19,20]. As a baseline, the default Monash 2013 tune [21] is used throughout. For further details on the various approaches to B (and D) hadron production, the reader is directed to [2,22], where a comparison of predictions obtained at NLO+PS accuracy (including both POWHEG and (a)MC@NLO [23,24] methods) and those obtained with the semi-analytic FONLL approach [25][26][27][28][29] are performed. In addition, information on predictions obtained in the so-called GM-VFNS scheme can be found in [30][31][32][33][34][35][36]. It is worth mentioning that while the calculation of next-to-NLO (NNLO) QCD corrections for massive [37][38][39] (and massless [40]) quark pairs are complete, and results for top quarks distributions have been presented in [41][42][43][44], the application of these results to B (and D) hadron final states is not yet available.
PDFs and α s . For the input PDFs, the n f = 5 fixed flavour number scheme (FFNS) PDF set NNPDF3.0 NLO α s (m Z ) = 0.118 [45] with 1000 replicas is used, and accessed through the LHAPDF6 interface [46]. The internal POWHEG routines are altered to extract α s from the grid provided with the PDFs, as oppose to using the internal α s routines. As discussed in [26], in such a set-up it is necessary to add compensation terms to the evaluation of the differential cross section which account for the mismatch in the running of both α s and PDF evolution with the fixed-order calculation -which is performed in a FFNS with n f = 3(4) for charm (bottom) quark pair production. These compensation terms are implemented in the POWHEG-HVQ library. The benefit of this approach is that the same PDFs are then used for both B and D hadron predictions, and the contributions from the resummed charm quark PDF are included in the B hadron predictions.
Scale variation. The dynamical reference scale (µ 0 ) is set to the transverse mass of the heavy quark in the underlying Born configuration (m T ). Scale variation is then performed by independently varying factorisation and renormalisation scales by a factor of two around the reference scale µ 0 with the constraint 1/2 < µ R /µ F < 2 (a 7-point scale variation).
Input masses. For the input heavy quark pole masses, the following choices for the central value and corresponding uncertainty are made These values are consistent within uncertainties with the recommendations of the HXSWG [47].
Fragmentation. In Pythia8, the heavy quark fragmentation is performed with the Lund-Bowler [48] approach -see for example [49]. The value of the fragmentation fractions, for example f (b → B − ), and distribution of the hadrons depends on the specifics of the particular tune. To investigate the dependence on the tune, distributions are also computed with the 4C tune [50]. In addition, the impact of manually varying the Lund-Bowler B quark fragmentation variable of r b = 0.855 in the default tune is also considered. It should be noted for normalised distributions and cross section ratios, the effects of varying fragmentation settings are negligibly small as compared to scale and PDF uncertainties.
In the most recent LHCb measurement of inclusive B hadron production [1], the measurement is performed for the sum of the (averaged over charge conjugate modes) following exclusive B hadron modes: In the case of Λ 0 b production, a correction factor of δ = 0.25 ± 0.10 was also applied to account for undetected Ω − b and Ξ b baryons. To match this definition, the B hadron final state is also taken as the sum of these four exclusive final states (including a weight of 1.25 for Λ 0 b baryons) and the total sum of these contributions is weighted such that f (b → B) = 1. In essence, Unless distributions are shown for specific B hadrons, this weighted sum is always applied to the B hadron final state.
Total uncertainty. To evaluate the total uncertainty of the 'NLO+PS' predictions (labelled this way in plots), the individual contributions from scale, m b and PDF variations are added in quadrature for both up and down variations as In addition, a more conservative 'Total uncertainty (linear)' will also be occasionally shown. This is computed by adding the scale uncertainty linearly with PDF and m b variations added in quadrature according to

Kinematics
The forward kinematic acceptance of the LHCb detector of η ∈ [2.0, 5.0] provides a unique opportunity to study heavy quark production in a kinematic regime beyond the reach of the central LHC detectors. As the heavy quark pair production process is dominated by gluon-fusion, such studies have the potential to probe the gluon PDF at both extremely small-and large-x values [2-4, 51, 52]. The sensitivity of such measurements is easily understood by considering the PDF sampling of the LO cross section where √ S is the hadronic CoM, and y b,b are the outgoing heavy quark rapidity. For both B and D hadron production, the LHCb detector has the capability to reconstruct hadrons  from p T > 0 (at small-m T ) and a large rapidities (y b ∼ 4.5) which provides sensitivity to low-x. Future measurements of B/D hadron production at large p T and y b also have the potential to probe the large-x gluon PDF [3].
To understand the kinematic region relevant for forward B hadron production, the LO B hadron cross section is shown in Fig. 1, differentially in x 1,2 . In the left plot (7 TeV), the B hadrons are required to be within either the pseudorapidity range η B ∈ [2.0, 2.5] (red) or η B ∈ [4.5, 5.0] (gray). Both of these kinematic regions are accessed in the recent LHCb measurement. As expected from Eq. (2.5), increasing the value of the pseudorapidity requirement simultaneously increases (decreases) the mean value of the x 1 (x 2 ) PDF sampling region. With the requirement of η B ∈ [2.0, 2.5], the mean PDF sampling occurs forx 1 ∼ 1.4·10 −2 andx 2 ∼ 5·10 −4 at a scale of Q 2 ∼ m 2 T ∼ 50 GeV 2 . In the right plot, the B hadron cross section is shown both at 7 and 13 TeV, where the B hadrons are required to be within the range of η B ∈ [2.0, 2.5]. The mean value of the transverse quark mass (m T ) which is probed for these selections is also highlighted. At fixed pseudorapidity, the mean values of the PDF sampling are decreased a factor ofx 13 i ≈ (7/13)x 7 i when increasing √ S from 7→13 TeV. It is worth noting that the region of the gluon PDF which is probed for these kinematic selections is well constrained (to a few %) by HERA DIS data [8,53,54].

(Normalised) B hadron cross section data
The purpose of this section is to perform a detailed study of the shapes of forward B hadron data available at both 7 and 13 TeV CoM energies [1,9]. There are two distinct LHCb data sets which will be considered in the following analysis.
• The first corresponds to the cross section measurement performed at both 7 and 13 TeV [1] for B hadrons reconstructed through the semi-leptonic decay modes B → DXµν. The cross section ratio measurements, which will be discussed in the following Section, is performed with this data set. These measurements are presented differentially in B hadron pseudorapidity, and inclusively with respect to transverse momentum. The motivation for considering the semi-leptonic decays is that the relevant branching fractions are well known, which results in a more precise determination of the absolute cross section rate. In contrast, the largest individual source of uncertainty for B hadrons reconstructed through B → J/ψX is associated to the branching fraction uncertainty.
• The second data set corresponds to the 7 TeV measurement of B hadrons reconstructed exactly through the decay B → J/ψX [9]. This measurement is performed for B + , B 0 , B s hadrons (and charge conjugate modes) where all decay products are reconstructed, and both the transverse momentum and rapidity dependence of B mesons are accessed.
Before starting the comparison to data, it is worth mentioning that the experimental precision of these absolute cross section measurements is ≈ 10−20%. In contrast, the NLO accurate predictions for the absolute cross section have large uncertainties of ≈ 50% -for the most part dominated by scale uncertainties. Consequently, a comparison of data to theoretical predictions at the level of the absolute cross section (although still important) is not particularly meaningful, since the overall normalisation of the cross section is uncertain. Instead, as discussed in detail in [2][3][4]52], it is often preferable to consider observables which are less sensitive to these scale uncertainties. The general approach of this Section will be to perform the comparison to data both at the level of the absolute and normalised cross section.

B → Dµν cross section data (7 and 13 TeV)
To begin, the recent forward B hadron cross section measurement [1] is studied, where the B hadrons have been identified through the exclusive semi-leptonic decays B → Dµν. As mentioned in the Introduction, this measurement is performed differentially in η B and inclusively in p B T , and the contributions from the sum of B + , B 0 , B s and Λ b hadrons (averaged over charge conjugate modes) as defined in Eq. (2.2) are included 1 .
The strategy for performing a comparison to this data will be to first normalise the differential cross section data with respect to the integrated fiducial cross section measurement, defined as The fiducial cross section data and corresponding theoretical predictions are summarised in Table 3.1, where consistency (within large theoretical uncertainties) with the predictions is found for both 7 and 13 TeV measurements. No correlation matrix has been provided for this B hadron measurement, and it is therefore assumed that the η B -independent systematic uncertainties (as reported in Table 4 of [1]) are fully correlated between the fiducial and differential data points. For the study of a normalised cross section, it would be beneficial to have access to the experimental bin-by-bin correlations for the cross section measurement.   The motivation for normalising the cross section in this way is that the large scale uncertainties in the absolute cross section are a result of varying the logarithmic scale dependence of the heavy quark transverse mass in the partonic cross section. However, this source of uncertainty primarily affects the overall normalisation of the cross section, and is highly-correlated between the neighbouring (pseudo)rapidity bins of the produced heavy quark. This observable is therefore theoretically more precise, and provides an important test of the shape of available data (rather than being overwhelmed by a normalisation uncertainty).
In Fig. 2 and 3, the distributions for the absolute (left) and normalised (right) LHCb B hadron cross section data is shown for 7 and 13 TeV respectively. For each plot, the predictions and data are also shown normalised to the central theory prediction in the lower panel. The total theoretical uncertainty for the normalised cross section data is below 10% while the absolute cross section uncertainty is close to 50%, demonstrating the above point. This approach also highlights an important feature of the data. For the case of the absolute cross section, the 7 TeV data tend to lie within the (large) theoretical uncertainties while the 13 TeV data tend to lie at the upper end of the theoretical scale uncertainties. At first glance, as the LHCb experiment reported, this may indicate that "The agreement with theoretical expectation is good at 7 TeV, but differs somewhat at 13 TeV". However, as shown by the normalised distributions, this behaviour is not indicated. Actually, perfectly good agreement is found for the shape of the normalised 13 TeV cross section data, while  the shape of the 7 TeV data is not as well described. This statement can be quantified by computing the χ 2 /N dat for the data points with respect to the central theory prediction, an approach which is justified for the normalised distribution as it has small theoretical uncertainties. This comparison gives While the χ 2 norm at 7 TeV is not particularly 'bad', it is substantially worse than that obtained at 13 TeV. The largest deviation is observed in the first bin, where the data is 2.1σ below the central theoretical prediction. It is worth mentioning that such a low value for the χ 2 norm at 13 TeV may indicate that the experimental uncertainties are overestimated. This suggests that it may be important to include bin-by-bin correlations when both normalising the data and computing the χ 2 norm . Finally, it is worth commenting on the behaviour of the absolute cross section at 13 TeV. In this case, it is observed that the absolute cross section tends to be on the upper end (although consistent within uncertainties) of the total NLO uncertainty band, which is dominated by the scale uncertainty. A similar trend has also been observed for D hadron production within the LHCb acceptance at 5, 7 and 13 TeV [5][6][7]. This behaviour is entirely consistent with the observation that the NNLO corrections to the absolute cross section for tt production (which, like cc and bb pair production is also dominated by the gluon-fusion partonic subprocess) at the LHC are large and positive [43].

B → J/ψX cross section data (7 TeV)
In addition to the η B dependent cross section data, a double differential (in p B T and y B ) cross section measurement was also performed at 7 TeV [9], where the B hadrons have been reconstructed through the decay B → J/ψX. It is useful to also consider the consistency of this data, to see if a similar trend is observed for the normalised cross section data. In this case, comparisons are performed for both double and single (p T -integrated) differential cross section data. When considering the rapidity distributions, the following normalisation is applied Like Eq. (3.1), this normalised distribution has the benefit that the uncertainty due to scale variation is highly correlated between numerator and denominator, since both are sensitive to similar values of m T . To construct the experimental distributions, it is assumed that the branching ratio and luminosity uncertainty are fully correlated between bins. With this exception, the experimental uncertainties are added in quadrature as the binby-bin correlations are also not available for this measurement. Both the experimental and theoretical rates for the fiducial cross section σ fid.
y are reported in Table 3.2.
LHCb data [µb] Table 2. Summary of 7 TeV fiducial measurements and predictions for the fiducial B hadron cross section σ fid.
y within the LHCb acceptance. The experimental uncertainties are statistical, systematic (including luminosity) and normalisation due to branching fraction uncertainties. The over normalisation of the theoretical prediction depends on the value of the fragmentation fraction f B used for each B hadron final state.
The comparison to data is shown in Fig. 4, where both the absolute (left) and normalised (right) B hadron rapidity distributions are shown. Again, in the lower panel both the theoretical predictions and data are shown normalised to the central theory prediction. In this case, the shown theoretical prediction corresponds to the B + hadron final state, where a fragmentation fraction of f (b → B u = 0.337) has been applied. Excluding the value of the fragmentation fraction, the individual distributions for B 0 and B s hadrons are extremely similar and are therefore not shown. Finally, in the case of the absolute cross for B s production, the experimental data has been multiplied by a normalisation factor of σ fid.
This normalisation is applied to allow the B s cross section to be compared with the other B hadron final states simultaneously. As demonstrated by this comparison, the individual measurements of all three B hadron are self consistent, and also consistent within uncertainties with the theoretical predictions for both the absolute and normalised cross section. However, there is some tendency for the data to undershoot the predictions in the region y B ∈ [2.0, 2.5].
In addition to the p T integrated distributions, a similar comparison can also be performed for the double differential data. This is done by normalising the data (for each p T bin) to that in the central rapidity bin y B ref. ∈ [3.0, 3.5] [52]. Therefore, for a given y B (i) and p B T (j) bin this observable is defined as   A comparison of selected data and predictions for this observable are shown in Fig. 5, where the double differential B + data (which is most precise) is compared with the corresponding theoretical predictions. The lowest rapidity y B ∈ [2.0, 2.5] region is shown in the left plot, while the neighbouring bin y B ∈ [2.5, 3.0] is shown in the right. In both cases, the predictions and data are normalised to the central value of the data in each bin. In the lower rapidity region of y B ∈ [2.0, 2.5] (and for p B T < 7 GeV), the data tends to systematically lie below the theoretical predictions. Although not shown, this behaviour is also observed for B 0 , and B s hadrons. In contrast, excellent agreement is found (within uncertainties) for the other rapidity bins, shown in the right plot for y B ∈ [2.5, 3.0]. While the tension at low y B ∈ [2.0, 2.5] is rather mild, it is worth mentioning that the experimental uncertainties (which have been constructed) are again likely over estimated since only the luminosity and branching ratio uncertainties are treated as correlated. The agreement with data could be better quantified with the experimental bin-by-bin correlations.
In summary, the 7 TeV B → J/ψX cross section data (both absolute and normalised) are consistent with the theoretical predictions presented differentially in p B T and y B . There is some tendency for the normalised B hadron data (observed for B + , B 0 and B s final states) to undershoot the theoretical predictions in the region y B ∈ [2.0, 2.5] and p B T < 7 GeV. This same trend is observed for the pseudorapidity dependent measurement at 7 TeV (but not at 13 TeV). It will be interesting to see if similar behaviour is observed in a corresponding 13 TeV measurement. In addition, as proposed in [2,3], it would be useful for the ratio of 13 and 7 TeV cross section measurements to be performed (double) differentially in p B T (and y B ).

Ratio of B hadron cross section data
The general motivation for considering a ratio of cross section measurements at different CoM energies is that the theoretical (and many experimental) uncertainties for a specific process are correlated between different CoM energies. Therefore, many sources of uncertainty partially cancel when constructing such a ratio. In some cases, this results in a dramatic reduction in scale uncertainties allowing sensitivity to PDFs, or both experimental and theoretical uncertainties may be reduced to an extent that these measurements can be used for luminosity determination of searches for the effects of physics beyond the Standard Model [56]. As mentioned in the Introduction, this method is particularly useful when considering B (and D) hadron production, as this is a process which is otherwise overwhelmed by large scale uncertainties. At the same time, the rate of the cross section growth with increasing CoM energy provides information on the shape of the gluon PDF at both small-and large-x.
To better understand the behaviour of the B hadron ratio data considered in this Section, it will be useful to introduce the following quantity which effectively describes the logarithmic growth of the gluon PDF with respect to x, and has recently been used to study the asymptotic behaviour of PDFs [57]. This is a useful quantity when considering the ratio of B or D hadron production measurements, since this observable is sensitive to exactly this growth. The computation of α eff. g (x, Q 2 ) for different PDF sets can be performed numerically using the LHAPDF interface, for which the PDF sets are provided as data files on grids in x and Q 2 space. The derivative in Eq. (4.1) can be performed at each x point on the grid by fitting a polynomial to the values of ln xg(x, Q 2 ) obtained for the neighbouring grid points in x. For the results shown in this work, a polynomial of order 3 is fitted to the central x point and the four neighbouring points in either direction. The results of this procedure are shown in Fig. 6, where both the gluon PDF (left) and α eff.
g (x, Q 2 ) (right) are shown for the baseline NNPDF3.0 NLO PDF, as well as the MMHT14 [58] and HERA2.0 [59] gluon PDFs. While not shown here, the effective exponents for the NLO gluon PDF from CJ15 [60], ABM11 [61] and CT14 [62] PDF fits exhibit the same behaviour as those shown. That is, at large-x (x ∼ 0.1) the gluon  Figure 6. Left: the evolved gluon PDF xg at the scale Q 2 = 50 GeV 2 . Right: the effective gluon exponent α eff. g also at the scale Q 2 = 50 GeV 2 . In both cases, the approximate PDF sampling regions of forward B hadron production are highlighted. PDF grows extremely quickly as it is generated by the valence PDF content, while at low-x the logarithmic growth becomes approximately constant. As demonstrated in Fig. 1, both large-and small-x regions are important for describing the forward B hadron ratio data. The remainder of this Section will be dedicated to studying various incarnations of cross section ratios.

Fiducial and differential ratio
Before discussing the differential data, it is instructive to first consider the ratio of the fiducial cross section measurements. This observable is defined as where the fiducial cross section σ fid. η B has previously been defined in Eq. (3.1). The experimental measurement and corresponding theoretical prediction are provided below In both cases, the breakdown of the various contributions to the uncertainty are provided. For the theoretical prediction, the scale uncertainties are still the dominant source of uncertainty, since the gluon PDF is predominantly sampled in the region which is constrained to a few % uncertainty. For example, when computing the 7 TeV fiducial cross section at LO, the mean sampling values arex 1 ∼ 3 · 10 −2 andx 2 ∼ 2 · 10 −4 at a scale of m 2 T ∼ 50 GeV 2 . The data is 2.7σ above the central theory prediction, and the predictions and data are consistent within their 2σ CL uncertainties. Although disfavoured by the baseline PDF set (NNPDF3.0), it is in principle still possible to accommodate this behaviour with a more steeply rising gluon PDF at low-x. This can be seen in Fig. 7 where the individual replica predictions obtained with the NNPDF3.0 NLO 1000 PDF replica set are shown.
The handful of outliers which are consistent with the LHCb data have exactly this feature. For example, replica member 200 (which is closest to the data) leads to a prediction of R fid. 13/7 = 2.16. However, it is worth mentioning that in the recent analysis of the forward LHCb D hadron data [4], this same replica member provided an extremely poor description of the normalised D hadron cross section data at 5, 7, and 13 TeV. The cross section ratio measurement by LHCb is also presented differentially with respect to η B , according to The comparison of the theoretical predictions to data for this observable are provided in Fig. 8 (left). In this case, the individual contributions from PDF and and m b uncertainties are also shown, and the more conservative 'linear' combination of uncertainties is provided.
In the lower pseudorapidity region of η B ∈ [2.0, 3.0] the scale uncertainties are dominant, and the PDF uncertainties become more significant at high pseudorapidity as the gluon PDF is probed at smaller values of x (see Fig. 1, right). The behaviour of the theoretical prediction is also easy to understand by examining Fig. 6. For increasing η B values, the ratio becomes more sensitive to the gluon PDF at larger (smaller) x 1 (x 2 ) values. In the low-x region, the logarithmic growth of the gluon PDF is approximately flat which results in an approximately pseudorapidity independent contribution to the ratio. At larger x values, the growth of the gluon PDF accelerates with increasing x, which results in a larger contribution to the ratio with increasing η B . However, this behaviour is clearly not observed in the LHCb data where the ratio is largest in lowest pseudorapidity region of η B ∈ [2.0, 3.0]. In fact, the first data point (η B ∈ [2.0, 2.5]) is 4.3σ above the central theory prediction, and 4.0σ with respect to the conservative upper theoretical uncertainty. The overall agreement with the data is extremely poor and, unlike the fiducial cross section, there are no individual replica PDF members which provide an adequate description of the data. This is shown in Fig. 8 (right), where the χ 2 /N dat of the differential ratio data is computed with respect to each of the  Figure 8. Left: The differential ratio (R 13/7 ) of inclusive B production with respect to pseudorapidity measured within the LHCb acceptance. Right: The χ 2 values obtained for each PDF replica in comparison to the differential ratio data.
1000 replica PDF members. The mean value is χ 2 /N dat = 36/6, and the minimum value (member 200) is χ 2 /N dat = 21/6. As a cross check of the theoretical predictions, it is important to study the perturbative stability of the ratio observable defined in Eq. (4.10). This is shown in Fig. 9, where both data and theoretical predictions are shown normalised to the central theory prediction. In this case, predictions are shown when LO matrix elements (M.E.) are used for the evaluation of the partonic cross section, which is then convoluted with the baseline PDFs (and α s ) evolved at either LO or NLO accuracy -the evolution is performed with the APFEL PDF evolution libraries [63]. This exercise demonstrates that the perturbative corrections, both through the evolution and the partonic cross section, are mild (each below 4%). It is therefore unexpected that NNLO QCD corrections would dramatically alter the theoretical predictions for this ratio observable. Another source of uncertainty not included in the total uncertainty is related to the treatment of the heavy quark fragmentation. The potential impact of this uncertainty has been assessed by varying the Lund-Bowler bquark fragmentation variable r b within the range of r b ∈ [0.67, 1.00], and by additionally showering events with the non-default Pythia8 Tune 4C. Further to this, the POWHEG events have also been showered with the Herwig7.0 PS [64,65] 2 . In all cases, the resultant ratio predictions differ from the central prediction by less than 3% within the region of η B ∈ [2.0, 5.0], which justifies not including this contribution in the total uncertainty.
To understand the origin of the tension observed in data, it will be useful to define kinematically shifted ratio observables which will be considered in the remainder of this Section.

Rapidity shifted differential ratio
As discussed in the previous Subsection, the behaviour of the differential ratio defined in Eq. (4.10) depends on both the behaviour of the gluon PDF at small-and large-x values.  This is because increasing √ S results in a shift of the PDF sampling in both x regions according tox 13 1,2 ∼ (7/13)x 7 1,2 . It is however possible to construct a ratio where either the small or large-x PDF regions are aligned. This can be achieved by introducing a rapidity shift between the kinematic region for which the numerator and denominator of the ratio are evaluated at. In heavy quark pair production, the PDF sampling depends on the outgoing rapidities of both heavy quarks -see Eq. (2.5). However, for a given value of the b quark rapidity y b , theb quark rapidity is symmetrically 3 distributed around y b such that on average yb = y b . Therefore, an alignment of the mean x sampling regions can be achieved by introducing the shift ∆y = ln 13 TeV 7 TeV = 0.62 . With this shift, one can specifically align (separate) thex 1 (x 2 ) sampling regions by introducing the observable where the rapidity shift is introduced in the numerator through dy B = dy B + ∆y. An example of this alignment is shown in Fig. 10 (left), where the LO B hadron cross section at 7 TeV is shown as a function of x 1,2 , integrated within the region of y B ∈ [2.0, 2.5]. The same cross section is shown at 13 TeV with the shifted integration region of y B ∈ [2.0 + ∆y, 2.5 + ∆y], demonstrating the alignment of the large-x regions. The benefits of introducing this shifted ratio are that the dependence on large-x region is eliminated in favour of sensitivity to the low-x region, since the low-x sampling regions are separated by a factor ofx 13 2 ≈ (7/13) 2x7 2 . At the same time, the theoretical uncertainties due to scale and m b variation are also reduced, since very similar values of m T are probed when evaluating the partonic cross section. This can be seen by examining the NLO+PS accurate predictions for the observable R 13/7 , which are provided in Fig. 10 (right). Although no data is currently available for this ratio (which requires the shifted kinematics), future analyses of the LHCb data would have access to this observable in the region of y B ∈ [2.0, 4.0]. Such a measurement would be very useful for understanding the tension observed in the η B dependent measurement.
An important feature of the shifted ratio observable is that the partonic kinematics which enter the evaluation of the partonic cross section become highly aligned. Consequently, the kinematic dependence of the ratio on partonic cross section is extremely mild, and this observable is essentially only sensitive to the growth of the low-x gluon PDF. This can be demonstrated by studying the variable  Fig. 10 (left). Explicitly,  approximation demonstrate that the shifted ratio is indeed directly sensitive to the growth of the low-x gluon PDF. Beyond LO, the dependence on the choice of the unphysical scale which enters the evaluation of the PDFs is evidently reduced by the mass factorisation terms present in the partonic cross section.
Of course, a comparison to data for the observable R 13/7 should be performed with respect to the most theoretically precise predictions (currently NLO), where the dependence on the choice of the unphysical scale is minimal. However, the reason for introducing the approximate relation R appr. 13/7 is to first demonstrate that R 13/7 is indeed directly sensitive to the growth of the low-x gluon PDF, but to also allow a qualitative study of the behaviour of the R 13/7 observable in terms of the quantity α eff. g (Q 2 , x). The point is that R appr.
13/7 (y B ) measures the growth of the gluon PDF across a given x range of x ∈ [x 7 TeV 2 (y B ),x 13 TeV 2 (y B )]. This is akin to the quantity α eff. g (Q 2 , x) which was introduced in Eq. (4.1) to study the logarithmic growth of the gluon PDF with respect to x, meaning that it is possible to approximately extract α eff. g from a differential measurement of R 13/7 (y B where the x dependence of α eff. g can be reconstructed in a similar fashion to what was done for the R appr. in Fig. 11 (right), and are compared to the quantity α eff.
g (x, Q 2 ) obtained directly from the PDFs and computed for the scale choices Q 2 = 12.5, 50.0, 200.0 GeV 2 . This method clearly allows the qualitative behaviour of α eff. g to be extracted from a measurement of R 13/7 (y B ).

Pseudorapidity shifted differential ratio
The LHCb measurement [1] is performed differentially in pseudorapidity bins of width ∆η B = 0.5, and it is therefore not possible to perform the alignment of the x regions as discussed above. However, a partial alignment can be performed by constructing a pseudorapidity shifted ratio according to (4.10) The success of this partial alignment is shown in Fig. 12, where the LO B hadron cross section is shown with respect to x 1,2 , integrated within the region of η B ∈ [2.0, 2.5]. The same cross section is shown at 13 TeV with the shifted integration region of η B ∈ [2.0 + ∆η B , 2.5+∆η B ]. The mis-match in the x 1 PDF sampling region is approximatelyx 13 1 (η B ) = 0.9x 7 1 (η B ), and results in the shifted pseudorapidity ratio having minor dependence on the behaviour of the large-x gluon. In the region of η B ∈ [2.0, 4.0], this mis-match is estimated to account for a flat correction factor to the ratio of 1.05. This 'correction factor' is obtained at LO by computing the valuesx 13 1 (η B ) andx 7 1 (η B ) for each pseudorapidity bin, and by then evaluating xg x 13 1 (η B ) /xg x 7 1 (η B ) . With this exception, the behaviour of this ratio (like the rapidity shifted ratio) is driven by the growth logarithmic growth of the gluon at low-x which is approximately flat below x ∼ 10 −3 -see Fig 6 (right).
A comparison of the LHCb data and the corresponding predictions of the pseudorapidity shifted ratio R 13/7 are shown in Fig. 13 (left). To obtain the experimental uncertainties, it is assumed that the same strength of correlation between η-independent systematics quoted for the 'non-shifted' ratio in Table 4 of [1] also applies to the shifted ratio. In this case, the data is again observed to exceed the theoretical predictions in the low pseudo-rapidity region of η B ∈ [2.0, 3.0], and there is a clear trend for the ratio to decrease with increasing η B .  To understand this behaviour in terms of the low-x gluon PDF, one can again consider the approximate relation between the shifted ratio and α eff. g introduced in Eq. (4.9). In Fig. 13 (right), the LHCb data has been extracted in a similar fashion to what was done for the rapidity shifted ratio in Fig. 11. In this case, the x dependence which enters the denominator of Eq. (4.9) throughx 13 TeV 2 (η B ) andx 7 TeV 2 (η B ) is extracted numerically in each pseudorapidity bin at LO. In addition, the flat 'correction factor' of 1.05 which accounts for the slight mis-alignment of the large-x region is also applied. For reference, this method is also applied to the NLO prediction in exactly the same way. The LHCb data clearly prefers large negative values of α eff. g around x ∼ 3 · 10 −4 , corresponding to an extremely fast growing gluon PDF, followed by a fast deceleration in the growth of the gluon PDF at lower-x values. The experimental and theoretical consistency of this behaviour will be discussed in the following Section.
Before continuing, it is important to emphasise that the extraction of α eff. g in this way is an approximation based on LO kinematics of the heavy quark production process. Nevertheless, this approach is still extremely useful for studying the qualitative features of the data. In this case, demonstrating that a significant change in the behaviour of the low-x gluon PDF is necessary to accommodate the data.

On the (in)consistency of the B hadron data
In the previous Section it was argued that, due to the kinematic alignment of the largex regions present for the shifted ratio observable, the large deviation observed in data necessarily points to a significant modification of the behaviour of the low-x gluon PDF. The purpose of this Section is to discuss both the theoretical and experimental consistency of this behaviour.

Theoretical consistency
From inspection of Fig. 13 (right), the LHCb data clearly prefers large negative values of α eff. g in the region of x ∈ [10 −4 , 10 −3 ] and Q 2 ∼ 50 GeV 2 which are inconsistent with the those obtained from global PDF fits. This is a region where the shape of the gluon PDF is governed by a combination of both perturbative effects (through DGLAP evolution) and non-perturbative effects (through the input PDFs at Q 0 ∼ 1 GeV). To investigate the origin of the tension between the values of α eff. g extracted from the LHCb data and those obtained from global analyses of proton structure (see Fig. 6), it is useful to introduce a toy model PDF set. Such an exercise is useful for understanding the perturbative behaviour of α eff.
g (x, Q 2 ) based on a simple model for the structure of the non-perturbative inputs PDFs.
To do so, such a model is introduced into APFEL [63] at the scale Q 0 = 1 GeV, based upon the following parametrisation of the input PDFs where q(x) = u(x) + d(x). The form of the input PDFs is motivated by non-perturbative QCD considerations [69,70], where Regge theory predicts the low-x behaviour x α , and Brodsky-Farrar quark counting rules predict the large-x behaviour (1−x) β . This functional form, now superseded by much more flexible parameterisations, has long been the starting point of PDF parameterisations. In the current model there are a total of 9 free parameters, two of which are fixed by the sum rules 3) The first sum rule is used to fix the normalisation of the valence content (N q ), and the second the normalisation of the gluon PDF (N g ). The exponents of the valence quark content are fixed to α q = 0.5 and β q = 3.0, and it is found that altering these values has little impact on the qualitative behaviour of the low-x gluon PDF. As a benchmark, it is assumed the sea quark and gluon distributions have identical shapes governed by α g = α S = −0.2 and β g = β S = 5.0, with the normalisation N S = 3/4N g . Several variations of the benchmark model are then considered by enforcing a vanishing sea or gluon content at the starting scale Q 0 . Practically, these scenarios are achieved by setting either N S = 0 or N g = 0, and correspond to generating the sea content or gluon PDF only perturbatively. In addition, variations of the component α g are also considered, which modify both the shape of the gluon PDF at low-x and also the normalisation (through Table 3. Summary of the various input parameters for the considered toy model PDF sets. the momentum sum rule). The default choices for these parameters and the considered variations are provided in Table 3.
The perturbative behaviour of α eff. g (x, Q 2 ) is the examined in this model by evolving the PDFs at NLO QCD accuracy using the APFEL evolution routines. In a similar fashion to how α eff.
g (x, Q 2 ) was computed for the LHAPDF grid files, the values of xg(x, Q 2 ) obtained for the toy models are tabulated on a grid in x space and the derivate defined in Eq. (4.1) is performed numerically. The results of this study are shown in Fig. 14 at the scale Q 2 = 50 GeV 2 , and are compared to the extracted values from the LHCb data. In addition to the predictions from the toy model variations, an analytic prediction based on the double asymptotic scaling of PDFs (DAS) [71][72][73] is also shown. The logarithmic growth of the gluon PDF in this case is provided analytically, based on an approximation valid in the double limit of large-Q 2 and low-x, according to The same general features are found for α eff. g (x, Q 2 ) in all cases. Firstly, the gluon PDF grows extremely quickly in the region x ∼ 0.1 as it is seeded by is valence-like PDF content at large-x. This growth then decelerates with decreasing x and eventually tends to a constant value, which depends on the choice input parameter α g . The general behaviour of α eff.
g (x, Q 2 ) for the toy model, based on the simplified parameterisation in Eq. 5.1, is therefore governed by a combination of the presence of valence-like PDFs at large-x and DGLAP evolution effects. These are the same features which are also observed in global PDF fits -see Fig. 6.
Based on these studies, it would seem the only way to accommodate the values of α eff.
g (x, Q 2 ) preferred by the LHCb data, is to introduce extremely ad-hoc behaviour in the non-perturbative gluon PDF in the range of x ∈ [10 −4 , 10 −3 ]. The reason is that the behaviour of α eff.
g (x, Q 2 ) at large-x is a general consequence of the valence-like content within the proton, which is well established. Therefore, to reach values of α eff.
g (x ∼ 5 · 10 −4 , Q 2 = 50 GeV 2 ) ≈ −0.7 as indicated by the LHCb data, it is necessary to introduce a region of accelerated growth in the non-perturbative gluon PDF around x ∼ 10 −3 . This period of accelerated growth must then be closely followed by a period of decelerated growth to accommodate the values of the ratio obtained at larger η B values. Introducing such a feature into the definition of g(x ∼ 10 −3 , Q 0 ) is in principle possible, since we should really be agnostic about the shape of non-perturbative object. However, the cost of doing so would be to drastically change the predictions of many collider observables. As an example, consider the prediction of the inclusive charm and bottom Structure Functions F qq 2 (x, Q 2 ), q = c, b, which are an ingredient of the cross section prediction for charm and bottom production in DIS. The LO prediction for this quantity is directly proportional to the gluon PDF, and is obtained through the convolution of the gluon PDF with heavy quark coefficient function. Measurements of both charm and bottom quark structure functions have been performed in the range of x ∼ [10 −4 , 10 −3 ] at Q 2 values of 6.5, 12.0, 25.0 GeV 2 [8].
No evidence for a steeply rising non-perturbative gluon PDF, which would result a sharp rise of both F cc 2 (x, Q 2 ) and F bb 2 (x, Q 2 ), is observed.

Experimental consistency with D hadron data
An another important consistency check can be performed by drawing comparison to the available forward D hadron data. The motivation for performing this check is that the theoretical framework for providing B and D hadron predictions is equivalent. In addition, there is LHCb data for D hadron production in a kinematic regime which is highly correlated with that of B hadron production 5 . Therefore, the consistency between D hadron predictions and data provides an important cross check of the B hadron results. Measurements of forward D hadron production have been presented at 5, 7, and 13 TeV [5][6][7], and as part of these measurements the ratio of double differential D hadron production at 13 TeV with respect to 5 and 7 TeV has been presented. These ratio measurements are available within the kinematic range of y D ∈ [2.0, 4.5] and p D T ∈ [0, 8] GeV. In the following, comparisons of the D hadron data are performed at the level of the double differential ratio according to As the tension in the B hadron data is observed in the lower pseudorapidity bins, the focus of the D hadron studies is also towards the low rapidity region of y D ∈ [2.0, 2.5]. In addition, particular attention should be paid to the region of p D T ∈ [6,8] GeV, where similar values of m T are probed with respect to the B hadron predictions. The first comparison is provided in the upper panel of Fig. 15 where the data and theoretical predictions for R D 13/5 are provided, normalised to the central value of the data. In the lower panel, the correlation of the D hadron ratio predictions with those of the B hadron ratio predictions (for specific choices of pseudorapidity bins) are shown. As expected, the correlation between these predictions is strongest in the high p T range, amounting to 0.4 and 0.7 for the B hadron ratio in the pseudorapidity bins of η B ∈ [2.0, 2.5] and η B ∈ [2.5, 3.0] respectively. As shown in Fig. 8 (left), these are the two pseudorapidity bins for which the tension in data is observed, however the predictions of R D 13/5 within this region are entirely consistent with the data. The same comparison is also performed for the experimentally less precise R D

13/7
data, and is shown in Fig. 16. In this case, the measured ratio systematically exceeds the theoretical predictions. This feature is exactly the same as that observed in the ratio of B hadron cross section measurements at 13 and 7 TeV, although the deviation is less significant in this case.
This situation is quite perplexing, as no deviation is found in the ratio of the 13 and 5 TeV cross section data, as shown in Fig. 15, which is expected be more sensitive to changes in the shape of the gluon PDF both at small-and large-x values. The fact that no deviation  Figure 16. Same as Fig. 15 (upper), now for the ratio of D hadron cross section data at 13 and 7 TeV.
is observed in this case, suggests that the D hadron data is not self consistent. Another way of viewing the tension in D hadron data is to construct the ratio R D 7/5 from the available cross section data. This is done by adding the experimental uncertainties in quadrature (a direct measurement of this ratio was not presented), and the results of this combination are shown in Fig. 17 for the rapidity region of y D ∈ [2.0, 2.5]. The experimental results for the ratio are generally below 1.0, which indicates that the differential cross section decreases with increasing CoM energy. These results are not in line with expectations based on perturbative QCD, where the evolved gluon PDF is expected to grow with decreasing x.

Discussion and conclusions
The large discrepancy observed in the ratio of forward B hadron production at 13 and 7 TeV has motivated a detailed of study of the available LHCb data, both through normalised cross section observables and various cross section ratios. It has been argued that, due to alignment of the PDF sampling regions, the tension present in the shifted pseudorapidity observable shown in Fig. 13 (left) must be attributed to solely to the behaviour of the low-x gluon PDF. In fact, it is possible to approximately relate the deviation observed in data to the logarithmic growth of the gluon PDF at low-x, described by the quantity α eff. g as defined in Eq. (4.1). After constructing this shifted ratio with the LHCb data, it is shown that the extracted values of α eff.
g are not consistent with the expectations from global PDF fits -see for example Fig. 13 (right). The reason for this tension is that the LHCb data indicates the presence of a region of accelerated growth of the gluon PDF, closely followed by a period of deceleration, within the kinematic range of x ∈ [10 −3 , 10 −4 ]. The only way to theoretically accommodate such behaviour is to introduce this structure into the nonperturbative gluon PDF, since this sort of feature is not generated by DGLAP evolution. However, introducing this behaviour would also lead to an extremely fast growth of the heavy quark structure functions F QQ 2 (x, Q 2 ) within this x-range, which is ruled out by measurements at HERA [8].
Studies of the normalised cross section data at 7 and 13 TeV, as shown in Fig. 2 and 3 respectively, do not conclusively indicate a problem with a particular data set. The overall agreement with the data sets is reasonable, as quantified in Eq. (3.2), while it is noted that there is local tension in the lowest pseudorapidity bin in the 7 TeV measurement of 2.1σ. No such tension is observed for the 13 TeV measurement. A further consistency check of the LHCb B hadron data is performed by studying the ratios of forward D hadron production available at 5, 7 and 13 TeV. This study was focussed on the rapidity region of y D ∈ [2.0, 2.5], corresponding to the region where tension is observed for the B hadron data. The measurement of of the 13/5 TeV D hadron cross section ratio (which is experimentally most precise), is found to be fully consistent with the theoretical predictions. Of the available D hadron ratio data, this observable is also the most correlated with the B hadron ratio, and is expected to be most sensitive to the shape of the low-x gluon PDF. However, a comparison to both D hadron ratios of 13/7 TeV and 7/5 TeV indicate a similar tension to what is observed for B hadron ratio. That is, this ratio exceeds the theoretical expectations in the lower (pseudo)rapidity region of [2.0, 2.5]. In particular, the reported cross section for the 5 TeV cross section measurement in this rapidity region is larger than that measured at 7 TeV.
To summarise, systematic tension is observed in the LHCb cross section measurements of B and D hadrons in the (pseudo)rapidity region of [2.0, 2.5]. Based on consistency checks of the data (through ratios and normalised distributions), this appears to be caused by a (pseudo)rapidity dependent efficiency correction which affects either 7 or both 5 and 13 TeV cross section measurements. If indeed this is the case, then the results analyses quantifying the impact of the LHCb B/D hadron data on proton structure [2,4,52] may also be affected. It is worth pointing out that the PDF constraints from this data are strongest in the large rapidity region, which seems to be a region which is least affected. Therefore, it is not likely that the results of these analyses would qualitatively change.
An extraction of the low-x gluon PDF obtained from analyses of forward B and D hadron requires reliable data. Given that a detailed understanding of both the magnitude and shape of the gluon PDF below x ∼ 10 −5 has important consequences for a range of physics processes such as LHC (and future collider) phenomenology [21,74,75], the predictions of atmospheric charm production [76,77], and the Ultra High Energy neutrinonucleon cross section [4,78]. It is therefore vital that LHCb re-investigate the measurements of forward B and D hadron production.