Production of forward heavy-flavour dijets at the LHCb within $k_{T}$-factorization approach

We calculate differential cross sections for $c \bar c$- and $b \bar b$-dijet production in $pp$-scattering at $\sqrt{s} = 13$ TeV in the $k_T$-factorization and hybrid approaches with different unintegrated parton distribution functions (uPDFs). We present distributions in transverse momentum and pseudorapidity of the leading jet, rapidity difference between the jets and the dijet invariant mass. Our results are compared to recent LHCb data on forward production of heavy flavour dijets, measured recently for the first time individually for both, charm and bottom flavours. We found that an agreement between the predictions and the data within the full $k_T$-factorization is strongly related to the modelling of the large-$x$ behaviour of the gluon uPDFs which is usually not well constrained. The problem may be avoided following the hybrid factorization approach. Then a good description of the measured distributions is obtained with the Parton-Branching, the Kimber-Martin-Ryskin, the Kutak-Sapeta and the Jung setA0 CCFM gluon uPDFs. We calculate also differential distributions for the ratio of $c \bar c$ and $b \bar b$ cross sections. In all cases we obtain the ratio close to 1 which is caused by the condition on minimal jet transverse momentum ($p_{T}^{\mathrm{jet}}>20$ GeV) introduced in the experiment, that makes the role of heavy quark mass almost negligible. The LHCb experimental ratio seems a bit larger. We discuss potentially important for the ratio effect of $c$- or $b$-quark gluon radiative corrections related to emission outside of the jet cone. The found effect seems rather small. More refine calculation requires full simulation of $c$- and $b$-jets which goes beyond of the scope of the present paper.


I. INTRODUCTION
The production of Heavy Flavours (HF) is a laboratory for studying effects of perturbative Quantum Chromodynamics (QCD). Very often it is used to "extract" parton distribution functions (PDFs) for gluon in the nucleon. This is usually done within collinearfactorization approximation. There are a few methods for experimental investigation of charm and bottom quark production at hadron colliders. One is the direct procedure based on full reconstruction of all decay products of open D and B mesons. The corresponding hadronic decay products can be used to built invariant mass distributions, permitting direct observation of D or B meson as a peak in the experimental invariant mass spectrum. Different charm mesons D ± , D 0 /D 0 , D ± s were used in this context at the LHC, including ALICE, ATLAS and CMS measurements at mid-rapidities [1-6] and the LHCb measurements at forward directions [7,8]. Similarly B 0 , B ± mesons were studied at ATLAS and CMS experiments [9][10][11].
Studies of open heavy flavour mesons require some additional modelling that takes into account effects related with a transition from quarks to hadrons. The transition is called hadronization or parton fragmentation and can be so far approached only through phenomenological models. In principle, in the case of multi-particle final states the Lund string model [12] and the cluster fragmentation model [13] are often used, however, those methods are originally devoted to mid-rapidities and their application in forward direc- tions is an open question. The hadronization of heavy quarks in non-Monte-Carlo calculations is usually done with the help of fragmentation functions (FF). The latter are similar objects as PDFs and provide the probability for finding a hadron produced from a high energy quark or gluon. Unfortunately, also in this case the FF technique leads to some unambiquities [14,15], especially in forward directions.
The inclusive open heavy flavour meson production at the LHC was studied differentially both in collinear [16][17][18][19][20] and k T -factorization approaches [21][22][23][24][25][26][27] as well as in the dipole model [19,28]. The inclusive distributions describe the experimental ones within theoretical uncertainties (choice of the heavy quark mass, choice of renormalization and factorization scales and choice of parton distributions). Also correlation measurements (meson-antimeson or meson-meson pairs) were performed by the LHCb [29] and CMS [30] collaborations. Those data were nicely explained in the framework of k T -factorization approach [21,25,26,31,32]. Similar studies were also performed with the Next-to-Leading Order (NLO) collinear factorization [33], but there some problems with description of the data were identified that can be solved only by introducing a phenomenological intrinsic transverse momentum or by matching parton-level calculations and parton shower [34]. In any case, theoretical calculations of the open heavy flavour meson production are to some extent biased by not fully understood fragmentation.
Another method is to measure heavy quark/antiquark jets. The reconstruction of jets containing heavy flavour hadrons provides more direct access to the primary heavy flavour parton kinematics than an inclusive measurement of heavy flavour hadrons and allows to study separately production and fragmentation effects. This is nowadays used almost routinely in the LHC era for b-orb-jets. Corresponding measurements have been performed by the ATLAS [35], the CMS [36] and very recently by the ALICE [37] experiments. The latter can reconstruct the b-flavoured jets down to extremely low transverse momenta, i.e. p T ≈ 10 GeV [37]. Recently, the CMS collaboration has measured for the first time jet shapes for b-jets in pp-collisions [38]. Also bb-dijets were measured by the ATLAS [39] and the CMS [40] experiments.
Although b-jet identification algorithms have been deployed for several decades, a task of identification of c-jets is more challenging [41]. Some first trials for inclusive c-jet studies were initiated by the CMS [42], ALICE [43] and LHCb [44] collaborations. Very recently the LHCb collaboration has achieved the identification of charm jets in Run 2 [45]. Last year, the LHCb has reported very unique results of simultaneous measurement of cc and bb-dijets in pp-scattering at 13 TeV [46]. This is the first cc-dijet differential cross-section measurement at the hadron collider.
Here we wish to explore the new LHCb data for cc and bb-dijet production within k Tfactorization approach. This data set is extremely interesting from the point of view of constraining unintegrated gluon densities in a proton. The data points were obtained for 2.2 < η jet1,2 < 4.2, p jet1,2 T > 20 GeV, jet radius R cone = 0.5 and for the difference in the azimuthal angle between the jets ∆ϕ > 1. 5. The specific kinematics set and characteristic mechanisms behind the considered reactions have many benefits in phenomenological treatment. In principle, because of large scales the pQCD methods are fully applicable and almost no heavy quark mass effects are expected. There is no need to use fragmentation functions and the production mechanism is fully dominated by gluon interactions.
Additionally, as will be discussed in the following, the forward direction explored in the LHCb measurement allows for a direct probe of different models of the unitegrated parton distribution functions (uPDFs) for gluon, simultaneously in both small-and large-x regions. All these aspects make the study very interesting from the context of present and future collider experiments on forward heavy flavour production, including those recently being proposed by the Forward Physics Facility (FPF) community [47]. It might also be of a great importance in understanding the dynamics standing behind mechanisms of prompt atmospheric neutrino flux production at present IceCube, Bajkal GVD or other future neutrino observatories. We remind the theoretical formalism for the calculation of the inclusive QQ-pair production in the k T -factorization approach [48], where Q = c, b stands for the charm and beauty quarks, respectively. In this framework the transverse momenta k t 's (virtualities) of both partons entering the hard process are taken into account, both in the matrix elements and in the parton distribution functions. Emission of the off-mass-shell initial state partons is encoded in the transverse-momentum-dependent (unintegrated) parton distribution functions (uPDFs). In the case of heavy-flavoured dijet production, the parton-level cross section is usually calculated via the Leading-Order (LO) g * g * → QQ fusion mechanism of off-shell initial state gluons that is the dominant process at high energies. As will be shown when presenting our numerical results, the q * q * → QQ mechanism remains subleading in the whole kinematical domain considered here. Then the hadron-level differential cross section for the QQ-dijet production, formally at leadingorder, reads: dσ(pp → QQ X) dy 1 dy 2 d 2 p 1,t d 2 p 2,t = d 2 k 1,t π d 2 k 2,t π 1 16π 2 (x 1 x 2 s) 2 |M off−shell g * g * →QQ | 2 (2.1)

The k T -factorization framework
where F g (x 1 , k 2 1,t , µ 2 F ) and F g (x 2 , k 2 2,t , µ 2 F ) are the gluon uPDFs for both colliding hadrons and M off−shell g * g * →QQ is the off-shell matrix element for the hard subprocess. The gluon uPDF depends on gluon longitudinal momentum fraction x, transverse momentum squared i,t + m 2 Q is the quark/antiquark transverse mass. The off-shell matrix elements are known explicitly only in the LO and only for limited types of 2 → 2 QCD processes (see e.g. heavy quark [49], dijet [50], Drell-Yan [51] mechanisms). Recently, higher multiplicity processes (2 → 3 and 2 → 4) with off-shell partons were also calculated at tree level in the context of cc+ jets [22,52,53], γ + c-jet [54], [55], 4-jets [56] and double charm [57] production. Some first steps to calculate NLO corrections in the k T -factorization framework have been made for diphoton production [58,59]. There are ongoing intensive works on construction of the full NLO Monte Carlo generator for off-shell initial state partons that are expected to be finished in near future [60]. Another method for calculation of higher multiplicity final states is to supplement the QCD 2 → 2 processes with parton shower. For the off-shell initial state partons it was done only with the help of full hadron level Monte Carlo generator CASCADE [61]. There, dedicated transverse-momentum dependent initial-state parton showers were introduced using backward evolution, that is not unique and needs to be matched to a given model of uPDFs.
Technically, there is a direct relation between a resummation present in uPDFs in the transverse momentum dependent factorization and a parton shower in the collinear framework. The popular statement is that actually in the k T -factorization approach already at leading-order some part of radiative higher-order corrections can be effectively included via uPDFs, without any additional showering procedure. However, it depends strictly on a theoretical construction of different uPDF models, in which extra emissions of soft and even hard partons can be encoded. In some uPDF models the off-shell gluon can be produced either from a gluon or quark, therefore, in the k T -factorization all channels driven by gg, qq and even by qg initial states are open already at leading-order (in contrast to the collinear factorization). Then, when calculating the heavy flavour production cross section via the g * g * → QQ mechanism one could expect to effectively include contributions related to an additional one or two (or even more) extra partonic emissions which in some sense plays a role of the initial state parton shower.
A. The hybrid model The LHCb measurement of the cc-and bb-dijets is performed within the asymmetric kinematical configuration. Under a general assumption that x 1 ≫ x 2 the cross section for the processes under consideration can be also expressed in the so-called hybrid factorization approach motivated by the works in Refs. [62,63]. In this framework the small-x gluon is taken to be off mass shell, the large-x gluon is treated as collinear and the differ-ential cross section e.g. for pp → QQX via g * g → QQ mechanism reads: where g(x 1 , µ 2 ) is a collinear PDF in one proton and F g * (x 2 , k 2 t , µ 2 ) is the unintegrated gluon distribution in the second one. The dσ g * g→QQ is the hard partonic cross section obtained from a gauge invariant tree-level off-shell amplitude. In the present paper we shall not discuss the validity of the hybrid model on the theoretical level and concentrate only on its phenomenological application in forward production. A derivation of the hybrid factorization from the dilute limit of the Color Glass Condensate approach can be found in Ref. [64].

The CCFM uPDFs
The CCFM evolution equation for gluon, in the limits of high and low energies (smalland large-x values), is almost equivalent to the BFKL and very similar to the DGLAP evolution, respectively [65]. In order to correctly treat gluon coherence effects it introduces the so-called angular-ordering which is commonly considered as a great advantage of this framework.
In the leading logarthmic approximation, the CCFM equation for unintegrated gluon density F g (x, k 2 t , µ 2 ) can be written as where µ 2 is the evolution (factorization) scale which is further defined by the maximal allowed angle for any gluon emission, The Sudakov and non-Sudakov form factors read: The first term in the CCFM equation is the initial unintegrated gluon density multiplied by the Sudakov form factor. It corresponds to the contribution of non-resolvable branchings between the starting scale µ 2 0 and the current running scale µ 2 . The second term describes the details of the QCD evolution expressed by the convolution of the CCFM gluon splitting function with the gluon density and the Sudakov form factor.
The theta function introduces the angular ordering condition.
The CCFM equation can be solved numerically using the UPDFEVOLV program [66], and the uPDFs for gluon and valence quarks can be obtained for any x, k 2 t and µ 2 values. Since only valence quarks are included in the evolution, the resulting CCFM gluon uPDF can be characterized as nearly 0-flavour scheme density. Therefore, it might be expected to be larger with respect to DGLAP-based unintegrated distributions that represent a variable-flavour-number-scheme (VFNS) framework.
Within the CCFM approach the parton transverse momentum is allowed to be larger than the scale µ 2 . This useful feature enables to effectively take into account higher-order radiative corrections, that correspond to the initial-state real gluon emissions which are resummed into the uPDFs.
The CCFM approach for the gluon uPDF was successfully used in the past to describe B-meson and bb-dijet production data taken by the D0 and CDF at Tevatron [67] as well as by the CMS at the LHC [25].

The PB uPDFs
The Parton Branching (PB) method, introduced in Refs. [68,69], provides an iterative solution for the evolution of both collinear and transverse momentum dependent parton distributions. Within this novel method the splitting kinematics at each branching vertex stays under full control during the QCD evolution. Here, soft-gluon emission in the region z → 1 and transverse momentum recoils in the parton branchings along the QCD cascade are taken into account simultaneously. Therefore, the PB approach allows for a natural determination of the uPDFs, as the transverse momentum at every branching vertex is known. It agrees with the usual methods to solve the DGLAP equations, but provides in addition a possibility to apply angular ordering instead of the standard ordering in virtuality.
Within the PB method, a soft-gluon resolution scale parameter z M is introduced into the QCD evolution equations that distinguish between non-resolvable and resolvable emissions. These two types of emissions are further treated with the help of the Sudakov form factors and with the help of resolvable splitting probabilities P (R) ba (α s , z), respectively. Here a, b are flavor indices, α s is the strong coupling at a scale being a function of µ ′2 , z is the longitudinal momentum splitting variable, and z M < 1 is the soft-gluon resolution parameter.
Then, by connecting the evolution variable µ in the splitting process b → ac with the angle Θ of the momentum of particle c with respect to the beam direction, the known angular ordering relation µ = |q t,c |/(1 − z) is obtained, that ensures quantum coherence of softly radiated partons.
The PB evolution equations with angular ordering condition for unintegrated parton densities F a (x, k t , µ 2 ) are given by [68] F a (x, k t , Here, the starting disitribution for the uPDF evolution is taken in the factorized form as a product of collinear PDF fitted to the precise DIS data and an intrinsic transverse momentum distribution in a simple gaussian form. Unlike the CCFM parton distributions, the PB densities have the strong normalization property: and therefore well reproduce modern collinear PDFs after integrating out their k tdependence. The PB uPDFs, including gluon and quarks, can be calculated by an iterative Monte-Carlo method and are characterized by a steep drop of the parton densities at k 2 t > µ 2 , again in contrast to the CCFM unintegrated distributions 1 . There are two basic sets of the parton-branching uPDFs: PB-NLO-2018-set1 and PB-NLO-2018-set2, that correspond to different choice of the parameters of the initial distributions [69]. Both of them, including uncertainties are available in TMDLIB [71] and both of them were used previously with success to describe the LHCb data on forward production of open charm meson [22,23].
According to this approach the unintegrated gluon distribution is given by the following formula This formula makes sense for k t > µ 0 , where µ 0 ∼ 1 GeV is the minimum scale for which DGLAP evolution of the conventional collinear gluon PDF, g(x, µ 2 ), is valid. A similar expression can be written also for quarks.
The virtual (loop) contributions may be resummed to all orders by the Sudakov form factor 1 Very recently, only a first attempt to incorporate CCFM effects into the PB method has been done [70].
which gives the probability of evolving from a scale k t to a scale µ without parton emission. The exponent of the gluon Sudakov form factor can be simplified using the following identity: P qg (1 − z) = P qg (z). Then the gluon Sudakov form factor reads , (2.12) where n F is the quark-antiquark active number of flavours into which the gluon may split. Due to the presence of the Sudakov form factor in the KMR prescription only last emission generates transverse momentum of the gluons initiating hard scattering.
In the above equation the variable ∆ introduces a restriction of the phase space for gluon emission and is crucial for the final shape and characteristics of the unintegrated density. In Ref. [72] the cutoff ∆ was set in accordance with the strong ordering (SO) in transverse momenta of the real parton emission in the DGLAP evolution, This corresponds to the orginal KMR prescription where one always has k 2 t < µ 2 F restriction and the Sudakov form-factor always satisfies the T g (k 2 t , µ 2 ) < 1 condition. The prescription for the cutoff ∆ was further modified in Ref. [73,74] to account for the angular ordering (AO) in parton emissions in the spirit of the CCFM evolution, (2.14) This modification leads to a bigger upper limit for k t than in the DGLAP scheme and opens the k 2 t > µ 2 F region. In this extra kinematical regime one gets T g (k 2 t , µ 2 ) > 1, which contradicts its interpretation as a probability of no real emission. Thus, the Sudakov form factor is usually set to be equal to one in that domain. For transparency, here the modified KMR model will be referred to as the Martin-Ryskin-Watt (MRW) model [74].
Different definitions of the ordering cut-off lead to significant differences between the two models. In the KMR model the k 2 t > µ 2 F region is forbidden while in the MRW case the k 2 t > µ 2 F contributions are directly allowed (see e.g. a detailed discussion in Ref. [79]). In the MRW model both in quark and gluon densities large k t -tails appear, in contrast to the KMR case. In the numerical calculations below we used the MMHT2014 [80] collinear PDFs to calculate both, the KMR and the MRW unintegrated densities.

The Kutak-Sapeta uPDFs
An alternative approach to those presented above has been applied in the Kutak-Sapeta (KS) model [63,81] for the gluon uPDF. There the unintegrated gluon density is obtained from the unified framework of the Balitsky-Kovchegov (BK) [82][83][84] and DGLAP evolution equations and then fitted to combined HERA data. This framework is a continuation of the model presented some time ago in Ref. [85] for unified BFKL-DGLAP evolution and it extends the ideas behind by taking into account the nonlinear (or gluon saturation) effects in the QCD evolution. In order to account for some effects related to the saturation of gluons, the unified BFKL-DGLAP evolution equation is supplemented with the nonlinear term in the BK form. Thus, the authors obtained the so-called modified BK equation in which the unified linear part deals with partial resummation of the NLLA corrections and the nonlinear term is taken in the basic LLA approach.
According to this approach the improved nonlinear equation for the unintegrated gluon density, written in momentum space, reads as follows: where ρ is the radius of the hadronic target, and F (0) g (x, k 2 ) is the starting distribution. The linear part of the equation is given by the BFKL kernel while the nonlinear part is proportional to the triple pomeron vertex which allows for the recombination of gluons.
The basic model of the KS gluon uPDF [63] has been further extended by introducing a factorization scale dependence [81] of the originally scale-independent density F g (x, k 2 ).
The scale dependence is obtained as follows: and the Sudakov form factor assumes the form: where ∆ = µ µ+k and P a ′ a is a splitting function with subscripts a ′ a specifying the type of transition.
In the numerical calculations below we applied both the KS-linear and KS-nonlinear gluon uPDF sets as implemented in the TMDlib. Phenomenological applications of these densities are limited to the x < 10 −2 kinematic domain and therefore they can be applied only when considering small-x effects. The KS gluon uPDFs are frequently used within the hybrid approach and are found very useful especially in phenomenological studies of forward particle production that is taking place in highly asymmetric kinematical configurations. Recently, the densities have been examined in the context of forward jet, forward-forward dijet as well as central-forward dijet production at the LHC (see e.g. Refs. [86][87][88][89]).

III. NUMERICAL RESULTS
In this section we present a variety of numerical results of the theoretical models described above for inclusive production of the cc-and bb-dijets at the LHC. The theoretical predictions are confronted with corresponding experimental data from the LHCb experiment [46] collected recently at √ s = 13 TeV. In the first step we show results obtained within the exact k T -factorization framework and then we carefully discuss the kinematics behind the processes under consideration. Next, we focus on the predictions based on the hybrid model showing explicitly the role of the large-x behaviour of gluon densities for a satisfactory description of the LHCb data and understanding the QCD dynamics behind. In the end we discuss the scale uncertainty of our predictions and R = cc bb dijet cross section ratio as a function of a few kinematical variables. For completeness of our study we also present a direct comparison to the corresponding results of calculations based on the collinear framework.

A. Predictions of the standard k T -factorization framework
As a default set in the numerical calculations below we take the renormaliza- it n (averaged transverse mass of the given final state) and the charm and bottom quark masses m c = 1.5 GeV and m b = 4.75, respectively. The strong-coupling constant α s (µ 2 R ) at NLO is taken from the MMHT2014 PDF routines.
For the CCFM uPDFs we always set the factorization scale to µ 2 = M 2 QQ + p QQ T , where M QQ is the invariant mass of the QQ system (or energy of the scattering subprocess) and p QQ T is the transverse momentum of QQ-pair (or the incoming off-shell gluon pair). This has to be applied as a consequence of the CCFM evolution algorithm.
Here and in the following only the dominant at high energies pQCD gluon-gluon fusion g * g * → QQ mechanism is taken into account. We have checked numerically that the annihilation q * q * → QQ mechanism is negligible here and can be safely neglected.

Charm dijets
We start with presentation of numerical results for forward production of the cc-dijets and F 2 data, however, even for the latter case, here the discrepancy between the predictions and the LHCb measurement is huge. Surprisingly, the Jung setA0 gluon density that a bit older set of the CCFM gluon uPDF leads to a much better description of the experimental data and seems to only slightly overestimates the measured distributions.   In Fig. 4 we show similar k T -factorization results but here we use the KMR (dashed histograms) and the MRW (solid histograms) gluon uPDFs. Both models seem to reasonably well describe the LHCb data. The histograms that represent the KMR predictions show a slightly better agreement with the data. It is visible especially for the leading jet η and ∆y * distributions which are displayed in linear scale. However, the effect is not significant for the overall picture and does not strongly disfavour the MRW gluon density here.
Next, we also examined the Parton-Branching gluon uPDFs against the LHCb data. In   gluon densities. In opposite to the CCFM and the KMR/MRW cases, here a slight tendency to underestimate the data points appears. However, both of the gluon densities provide a very good description of the measured distributions and both lead to a very similar results. A small missing strength is observed only at smaller values of the leading jet η and some differences between the two predictions are found only at large leading jet p T 's and large dijet invariant masses M cc-dijet . The PB-NLO-set2 density seems to reproduce slopes of the distributions a bit better than the PB-NLO-set1 one. The two models of the gluon uPDF differ from each other in a value of the starting evolution scale as well as in a set of the argument of α S . In the set1 case, the integrated parton density and the initial parameters are the same as the ones obtained by HERAPDF2.0 PDF. In the set2 case, the parameters are slightly modified and thus the integrated parton distributions do not fully reproduce the collinear ones. However, in both cases a reasonably good fit to HERA data was obtained.  histograms for charm and bottom dijets is extremely small. We will discuss this issue in more detail in the following, when presenting results for the charm to bottom dijet ratio Taking into account all the statements above, our basic findings about the application of different uPDF models in description of the LHCb cc-dijet data applies also here. The main conclusion from a comparison of our predictions now with the bb-dijet data is that the CCFM and the MRW/KMR gluon densities lead to a larger discrepancy between predictions and the data, overshooting the bb-dijet data points even more stronger than in the cc-dijet case. On the other hand the Parton-Branching densities describe the bottom data better than the charm one, where we have found a small missing strength, but here this is not the case.

B. Kinematics of the process probed by the LHCb experiment
In the previous section we have shown that within the k T -factorization approach some of popular models of the gluon uPDF have difficulties with a reasonable description of the LHCb data on forward production of heavy flavoured dijets. In order to identify the possible reasons of the failure we wish to carefully illustrate the kinematics behind the considered process as probed by the LHCb experiment.
In Fig. 9 we present the double differential cross sections for forward production of cc-dijets in pp-scattering at √ s = 13 TeV probed in the LHCb experiment as a function of the longitudinal momentum fractions log 10 (x 1 ) × log 10 (x 2 ) (left panels) as well as the initial gluon transverse momenta k t1 × k t2 (right panels   the maximum of the cross section around 10 −1 and x 2 is much smaller with the maximum in the 10 −4 < x 2 < 10 −3 range. At the same time, the large-x gluon has rather small transverse momenta k t1 5 GeV, while the small-x gluon takes much larger values of the k t2 with the maximal range depending on the model of the gluon uPDF. The largest tail in k t2 is observed for the case of the CCFM density, while the Parton-Branching model leads to the smallest one. The correlation between x and k t is plotted explicitly in Fig. 10. It is also interesting to see how the initial transverse momenta of incoming gluons k t contribute to the leading jet p T . The mutual relation is presented in Fig. 11. The leading jet p T distribution in the range of 20 < p T < 50 GeV is driven by both gluons incident transverse momenta. Above 50 GeV mostly the large-x gluon contributes to the spectrum.  Fig. 9 but here as a function of the initial gluon transverse momenta and the longitudinal momentum fractionsk t1 × log 10 (x 1 ) (left panels) as well as k t2 × log 10 (x 2 ) (right panels). the large-x behaviour is challenging and might not be under full theoretical control. At large-x one might expect that a given uPDF model should closely reproduce well known collinear PDFs when integrating out the k t -dependence. Here we wish to take a closer look on this issue.
In Fig. 12 we present the gluon uPDFs integrated over the transverse momentum k t as a function of the longitudinal momentum fraction x at a given scale µ = 25 GeV, rel-  Fig. 9 but here as a function of the leading c-jet transverse momentum and the initial gluon transverse momentap T × k t1 (left panels) as well as p T × k t2 (right panels). evant for the process under consideration. The left and right panels correspond to the Parton-Branching and the CCFM gluon uPDFs, respectively. As a reference the collinear MMHT2014nlo68cl gluon PDF is shown. We found that the Parton-Branching uPDF reproduces very strictly the collinear PDF, while for the CCFM uPDFs this is not the case.
In Fig. 13 we show in addition the gluon uPDFs integrated over the transverse mo- gluon uPDFs, respectively. The top panels are obtained for x = 10 −1 and the bottom ones for x = 10 −4 . As a reference the collinear MMHT2014nlo68cl gluon PDF is also shown.
Here we see that for large-x gluon the relation between the integrated uPDFs and the collinear PDF does not change with the scale µ. On the other hand, the integrated gluon uPDFs grow faster with the scale than their collinear counterparts. This effect is stronger especially in the case of the JH2013 CCFM uPDF sets.
Having shown somewhat technical problems within full k T -factorization approach in the following we shall consider also the hybrid factorization approach.

D. Predictions of the hybrid factorization framework
All the calculations below preformed within the hybrid model are done using the collinear MMHT2014nlo68cl gluon PDF on large-x side. We have checked numerically that at large-x there is no significant difference between different sets of the gluon PDF at NLO, given in the literature by different PDF groups. gluon PDF is shown in addition.

Charm dijets
The results of the hybrid model will be shown in an analogous manner as for the full k T -factorization. We start with forward production of the cc-dijets in pp-scattering at √ s = 13 TeV. In Fig. 14 we show the corresponding differential cross sections as a func-tion of the leading jet η (top left panels), the rapidity difference ∆y * (top right panels), the leading jet p T (bottom left panels) and the dijet invariant mass M cc-dijet (bottom right panels). The theoretical histograms are obtained for the interaction of off-shell small- x and on-shell large-x gluon distributions using three different sets of the CCFM uPDFs: JH2013set1 (dashed), JH2013set2 (solid), and Jung-setA0 (dotted). We observe that within the hybrid approach the discrepancies between predictions and the experimental data decrease with respect to the case of the full k T -factorization calculations. Within the hybrid factorization and the Jung setA0 gluon uPDF we get an excellent description of the experimental distributions. Also here the JH2013 sets of the CCFM gluon uPDFs overestimate the data points by less than a factor of 2 only, which is also a significant improvement.
A similar effect is observed also in the case of the KMR/MRW uPDF. As can be seen in Fig. 15, using the hybrid approach leads to better agreement with the data and removes a slight tendency to overshoot the data, reported previously in the full k T -factorization case.
On the other hand, the hybrid model results obtained with the Parton-Branching uPDFs show an opposite effect (see Fig. 16). Here we observe a small enhancement of the cross section but it works in a good direction and improves the agreement with the experimental data, resulting in an excellent description of all four measured differential distributions.

Bottom dijets
Now let us go to bottom dijet production within the hybrid approach. Here we repeat the analysis from the previous subsection but for the bb-dijet LHCb data (see Figs. 18,19,20,and 21). As was already mentioned our predictions are almost independent of the heavy quark mass so the theoretical distributions are almost identical to those from the previous subsection, however, the experimental data points are changed. The Jung setA0, the KMR, and the KS-linear gluon uPDFs provide a very good description of the charm dijet data within the hybrid model, here showing a visible tendency to overshoot the bottom dijet data. Again the PB-NLO-set1 uPDF leads to the best agreement with the data, however, here some problem appears for the very first bin in the leading jet transverse momentum where the theoretical cross section slightly overestimates the experimental point. The discrepancy affects also the leading jet η and the rapidity difference To summarize this section, we conclude that the best agreement between the LHCb heavy flavoured dijet data and the predictions of the hybrid model can be obtained for the PB-NLO-set1 uPDF. However, the results of the model show a small tendency to underestimate and overestimate the experimental distributions for the cc-and the bbdijets, respectively. It may suggest, that within our model the production cross section ratio R = cc bb is not well reproduced.

E. The charm/bottom dijet ratio
The last conclusion above can be directly examined beacuse the LHCb collaboration has also presented the corresponding distributions of the ratio R = cc bb as a function of the leading jet η, the rapidity difference ∆y * , the leading jet p T and the dijet invariant mass M QQ-dijet . A comparison of our predictions to the experimental distributions is shown in Fig. 22. The theoretical ratio is identical for all models of the gluon uPDF used here. The predicted ratios are slightly above 1 and are almost independent of the considered kinematic variables. The experimental results correspond to a rather larger values. In general, we obtain a theory/data agreement, however, we reached only the lower experimental limits defined by large total experimental uncertainties. We will continue the discussion of the charm to bottom ratio in the following. Here the central set of the scales corresponds to the µ 2 = m 2 t case. As we can see the sensitivity is sizeable but not huge and seems to be decreasing with the leading jet transverse momentum.
In Figs. 25 and 26 we present in addition how our results depend on the definition of the central set for the scales. The sensitivity has a similar size as above and one has to keep in mind that different choices of the central sets might slightly modify the overall picture. For example, using the µ 2 = M 2 QQ we get slightly smaller cross sections which means that within this choice we get better agreement with the data for bottom than for charm production -in opposite to conclusions from our default scale set.
What we find interesting here, is that playing with the definition of the central scale does not improve our description of the charm to bottom production cross section ratios (see Fig. 27).

G. A comparison with the results of the collinear approach
Finally, we wish to compare our predictions based on the k T -factorization framework with those obtained according to the collinear approximation. In Figs. 28, 29, and 30 we plot our hybrid model results for the PB-NLO-set1 uPDF (solid bands and histograms) and results of two different collinear approach calculations: MADGRAPH5 aMC@NLO [90] (dashed bands) and PYTHIA8 [91] (dotted histograms), both taken from Ref. [46].
The two models correspond to NLO and LO matrix element calculations, respectively, supplemented further with dedicated parton-showers.
The main conclusion is that all three models presented here lead to a rather consistent

IV. A COMMENT ON FINITE JET SIZE EFFECTS
So far we have calculated distributions at the parton level. The agreement with experimental data is, however, quite good, both for cc and bb dijets. In principle, one could worry about jet size effects. The LHCb collaboration chosen the jet cone size R cone = 0.5 in their analysis [46]. Could this affect our partonic results?
In general, a thorough answer to this question requires modelling of c and b jets which is rather complicated and goes beyond the scope of the present paper. Instead of following this path, we will try to estimate a possible effect as was done for one jet case in [92].
It is energy of the jet which decides about the shape of the jet [92]. The larger energy of the jet the smaller finite cone size effects. In Fig. 31 we show distribution in jet energies for the LHCb kinematics. The typical jet energies are relatively large E 1 , E 2 > 200 GeV. In the approach from [92] the basic ingredient is dE g /dxdk 2 where E g is energy of the emitted gluon, x is its longitudinal momentum fraction and k 2 is transverse momentum squared of the gluon with respect to the heavy quark (c or b in our case). This quantity depends on the mass of the emitting quark. In principle, one can expect even a spectacular effect known under the name dead cone [93]. The dead cone effect was observed experimentally recently for the first time by the ALICE collaboration [94]. Energy loss due to emission outside of the jet cone can be then written as where θ open (R cone ) is jet openning angle corresponding to a given jet radius R cone . The distribution dE g dxdk 2 is obtained [92] from a generalization of the Gunion-Bertsch formula [95] for the gluon number distribution dn g dxdk 2 .
At high energies, as in our case, such effects are completely unimportant as dead cone angle is of the order of a small fraction of one degree. In our case, it is the energy (transverse momentum) which escapes from experimentally defined jet cones. We have estimated that in our case the relative energy loss ∆E/E, including range of jet energies and R cone = 0.5, is of the order of a few percent. It is reasonable to expect the same is true for ∆p T /p T . The main effect is due to the fact that the cuts are imposed on the measured transverse momenta, not momenta of heavy quarks/antiquarks which are of course bigger than the measured ones. In our case we have to apply this procedure to both measured jets. This leads to a damping of the cross section by approximately a few percent. We have checked that the damping is practically the same for c/c and b/b quarks/antiquarks. In our energy range the quark mass effect is negligible. We conclude that jet cone size effect is most probably not responsible for the charm-to-bottom ratios discussed in our paper. We think that the effect of the shape of the heavy quark jets requires further studies, also in the present context. In Fig. 32 we show results for different (independent off transverse momentum and pseudorapidity) values of the leakage of the energy outside of the jet cone radius R cone = 0.5. Here we have assumed in addition ∆p T p T ≈ ∆E E for each of the two jets. Although the assumed leakage (1%, 3%, 6%) is rather small the resulting effect on the cross section normalization is rather sizeable.

V. CONCLUSIONS
In the present paper we have studied the cc-and bb-dijet production in pp-scattering for √ s = 13 TeV and LHCb acceptance limitations on jet pseudorapidities 2.2 < η jet1,2 < 4.2, transverse momenta p jet1,2 T > 20 GeV, jet cone R cone = 0.5 and on the azimuthal angle between the jets ∆ϕ > 1.5. The differential cross sections were calculated within k Tfactorization and hybrid approaches. We have used different unintegrated parton distribution functions for gluon, including the PB, the KMR/MRW and various CCFM models. We have calculated distributions in leading-jet transverse momentum p T , leadingjet pseudorapidity η, rapidity difference bewteen jets ∆y * and in dijet invariant mass M QQ-dijet . We found that an agreement between the predictions and the data within the full k T -factorization is strongly related to the modelling of the large-x behaviour of the gluon uPDFs which is usually not well constrained. For this case only the PB-NLO-set1 gluon uPDF leads to a satisfactory description of the LHCb data. Rest of the models seem to visibly overestimate the experimental data points. This was understood as due to incorrect effective integrated gluon distribution obtained from the unintegrated one for large-x.
The problem may be avoided by following the hybrid factorization. Then a good description of the measured distributions is obtained with the PB-NLO-set1, the KMR-MMHT2014nnlo, the Kutak-Sapeta and the Jung setA0 CCFM gluon uPDFs for both cc and bb production taking into account experiemental uncertainties. Only the most recent JH-2013-set1 and JH-2013-set2 CCFM gluon uPDFs seems to overestimate the data sets even within the hybrid approach. In any case we observe a small tendency of our predictions to slightly overestimate the data for bb-dijets while the corresponding data for cc-dijets can be perfectly described within the same gluon uPDFs. In general, taking into account theoretical uncertainties our hybrid model predictions are consistent with collinear results of the PYTHIA8 and the aMC@NLO frameworks.
We have presented also the ratios of cc to bb cross sections. We get cross sections ratio very close to 1. This could be understood by the fact that for large transverse momenta the effect of quark mass is rather small. Slightly larger, but with large systematic error bars, effect was obtained in the experimental extraction. It is difficult to understand in the moment the possible disagreement. In principle, it can be due to inappropriate extraction of c/c jets (e.g. by admixture of light quarks and/or b/b jets). We have also discussed the effect of damping of the cross section within a simple model due to finite jet radius.
It may be expected that statistically a part of the parton (gluon) energy escapes outside of the jet cones. Within this, a bit naive approach, one may expect sizeable corrections.
A much better approach would be to take into account internal structure of c-and b-jets due to parton shower. This definitely goes beyond the scope of the present paper.
In our k T -factorization approach the production of c andc as well as b andb is identical (symmetric). Asymmetry effects ( dσ(c) dξ = dσ(c) dξ and dσ(b) dξ = dσ(b) dξ , where ξ represents schematically y and p T ) were discussed in Ref. [96]. The asymmetry effects found there are rather small, of the order of 1 %. Also the contribution of electroweak processes, as discussed in Ref. [96], is rather small. So we think these interesting effects are not crucial for our studies.