Consistent treatment of charm production in higher-orders at tree-level within $k_T$-factorization approach

We discuss production of $c \bar c$-pairs within $k_T$-factorization approach (off-shell initial partons) with unintegrated parton distribution functions (uPDFs). We present a consistent prescription which merges the standard leading-order (LO) $k_T$-factorization calculations for this process with tree-level next-to-leading order (NLO) and next-to-next-to-leading order (NNLO) matrix elements. For the first time we include in this framework 2 $\to$ 3 and 2 $\to$ 4 processes with extra partonic emissions for single particle distributions as well as for correlation observables. The use of the KMR uPDF leads to a good description of the existing charm ($D$-meson) data already at the leading-order. On the other hand, a new Parton-Branching (PB) uPDF strongly underestimates the same experimental data. A direct inclusion of the higher-orders at tree-level leads to an overestimation of the data, especially for the KMR uPDF. This suggests a significant double-counting. We propose a simple method how to avoid the double-counting. Our procedure leads to a much better description of the experimental data when including the higher-order contributions. Then with the KMR uPDF we get similar results (both for single particle and correlation observables) as for the standard calculations of the 2 $\to$ 2 processes. For the PB uPDF inclusion of the higher-orders considerably improves description of the experimental data. We conclude that the LO calculation with the KMR uPDF effectively includes the higher-orders which is not the case for the PB uPDF.


I. INTRODUCTION
The production of heavy flavours is known to be a good example of perturbative QCD calculations -the quark mass sets already a sizeable scale. Charm quark is the lightest heavy quark where one believes in the pQCD treatment. Leading-order collinear approach gives much too small cross section, compared to the experimental cross section for charmed meson production. Clearly higher-order corrections are needed. Next-toleading order approach was developed for inclusive variables only (single charm distribution). In general, one is interested not only in single charm distributions but also in correlation observables. Some studies of correlation observables were done e.g. in Refs. [1,2] within k T -factorization approach.
In the present paper we will discuss both single particle distributions (distributions in rapidity or charm transverse momentum) and correlation observables (distribution in azimuthal angle between cc, invariant mass distribution, transverse momentum of the cc pair, difference in rapidity between c andc).
We propose and discuss a consistent prescription which merges the standard leadingorder (LO) k T -factorization calculations for this process with tree-level next-to-leading order (NLO) and next-to-next-to-leading order (NNLO) matrix elements. The applied procedure was originally proposed in the context of BB pair production within the Parton-Reggeization-Approach (PRA) in Ref. [3]. There, the LO calculations for off-shell initial state partons were supplemented by the NLO corrections from the emission of one additional hard gluon. The consistent merging procedure constructed there can be implemented for the calculation of charm production within the k T -factorization approach. In this paper, we basically follow the ideas presented in Ref. [3] but we extend those studies to the case of NNLO corrections from the emission of two additional hard gluons.
Our new scheme for the calculations provides a possibility to study the charm cross section differentially beyond the NLO collinear approaches, FONLL [4] and GM-VFNS [5], commonly used as the state of the art in this context. We expect the 2 → 4 contributions, that are missing there, to be of the special importance for the large transverse momenta of charm particles studied at the LHC. We recall the standard theoretical formalism for the calculations of the cc-pair production in the framework of the k T -factorization [6]. This approach is commonly known to be very efficient not only for inclusive particle distributions but also for studies of kinematical correlations. According to this approach, the transverse momenta k t 's (virtualities) of both partons entering the hard process are taken into account. In the case of charm (or in general heavy) flavour production the parton-level cross section is usually calculated via the 2 → 2 leading-order g * g * → cc fusion mechanism of off-shell initial state gluons that is dominant processs at high energies. Emission of the initial state partons is encoded in the so-called unintegrated parton distribution functions (uPDFs). Then the hadron-level differential cross section at the tree-level for the cc-pair production reads: dσ(pp → cc 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 * →cc | 2 (2.1) 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 * →cc is the off-shell matrix element for the hard subprocess. The gluon uPDF depends on gluon longitudinal momentum fraction x, transverse momentum squared k 2 t of the gluons entering the hard process, and in general also on a (factorization) scale of the hard process µ 2 F . The extra integration is over transverse momenta of the initial partons. Here, one keeps exact kinematics from the very beginning and additional hard dynamics coming from transverse momenta of incident partons. Explicit treatment of the transverse momenta makes the approach very efficient in studies of correlation observables. The two-dimensional Dirac delta function assures momentum conservation. The unintegrated (transverse momentum dependent) gluon distributions must be evaluated at: where m i,t = p 2 i,t + m 2 c is the quark/antiquark transverse mass. In the case of charm quark production at the LHC energies, especially in the forward rapidity region, one tests very small gluon longitudinal momentum fractions x < 10 −5 [1].
The off-shell matrix elements are known explicitly only in the LO and only for limited types of QCD 2 → 2 processes (see e.g. heavy quarks [7], dijet [8], Drell-Yan [9]). The calculation of higher-order corrections in the k T -factorization is much more complicated than in the collinear approximation approach. Some first steps to calculate NLO corrections in the k T -factorization framework have been tried only very recently for diphoton production [10,11]. 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 [12]. 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 CAS-CADE [13]. However, this method can be consistently used only with dedicated models of uPDFs.
On the other hand, the popular statement is that actually in the k T -factorization approach with 2 → 2 tree-level off-shell matrix elements some part of real higher-order corrections can be effectively included. This is due to possible extra emissions of soft and even hard partons encoded in the uPDFs. In this sense, when calculating the charm production cross section via the g * g * → cc mechanism one could expect to include some contributions related to an additional one or even two extra partonic emissions, effectively taking into account, as an example, the gcc and ggcc final states. However, the presence and a size of the extra emissions strongly depends on the internal construction of the unintegrated parton distributions. The extra emissions (from the uPDFs) are expected to be very important, especially for studies of kinematical correlations. The correlation observables, such as the azimuthal angle difference of c andc or transverse momentum of the produced system are very useful for testing transverse momenta of initial partons and may be helpful in limiting uPDFs uncertainties and in understanding their evolution.
Some time ago we showed that in the case of charm production at the LHC, within the above formalism only the Kimber-Martin-Ryskin (KMR) uPDF leads to a reasonable decription of the experimental data for D-meson and DD-pair production [1]. This was further confirmed by other authors [2]. As also discussed in Ref. [14] the k T -factorization approach at leading-order with the KMR uPDF leads to results well consistent with collinear NLO approach. The KMR uPDF is known to allow by its construction for a large contribution from the k T > µ F kinematical regime. Effectively, this extra emission of hard partons (gluons) from the uPDF corresponds to higher-order contributions. As reported in Ref.
[1] the rest of the commonly-used models of the uPDFs from the literature is rather missing those contributions and significantly underestimates the experimental data on charm production at the LHC.
In the numerical calculation below, based on the standard k T -factorization framework, we apply the KMR uPDF in its original form. The KMR distributions are calculated from the MMHT2014 [15] and CT14 [16] gluon PDFs. As a default set of the calculations we use the renormalization and factorization scales iT n and charm quark mass m c = 1.5 GeV. The uncertainties related to the choice of the collinear PDF and of the renormalization/factorization scales will be discussed when presenting numerical results.
The transition of charm quarks to open charm mesons is done in the framework of the independent parton fragmentation picture (see e.g. Ref. [17]). Here we follow the stan- where p t,c = p t,D z and z is the fraction of longitudinal momentum of charm quark c carried by a meson D. In the numerical calculations we take the Peterson fragmentation function [18], often used in the context of hadronization of heavy flavours. Then, the hadronic cross section is normalized by the relevant charm fragmentation fractions for a given type of D meson [19]. This is a direct analogy to the issue of merging hard emissions from higher-order matrix elements with soft emissions from the parton showers [20]. Due to the lack of the full NLO and/or NNLO framework of the k T -factorization, within the present methods this can be done only at tree-level. In the proposed scheme, we include and sum up the dominant 2 → 2 , 2 → 3 and even 2 → 4 contributions to heavy quark-antiquark pair production under a special conditions introduced to avoid a possible double-counting. In this model, the higher-orders with hard extra emissions come from the higher-order matrix elements, while only the softer extra emissions are included via the uPDF. Therefore, within this method for studies of heavy flavour production one could apply different models of uPDFs that do not include in their evolution a sufficiently hard extra emissions.
In this model we calculate the g * g * → cc, g * g * → gcc and g * g * → ggcc mechanisms for off-shell initial state partons. We have check numerically, that for the LHC energy the channels driven by gluon-gluon fusion are the dominant ones for each of the considered reactions. For the present studies (high-energy collisions) the contributions from the quark-induced processes can be safely neglected. The numerical calculations for the considered higher-order contributions (see Fig. 2) are performed in the framework of the k T -factorization approach within the methods adopted in the KaTie Monte Carlo generator [21]. The off-shell matrix elements for higher final state parton multiplicities, at the tree-level can be calculated analytically applying well defined Feynman rules [22] or recursive methods, like generalised BCFW recursion [23], or numerically with the help of methods of numerical BCFW recursion [24]. The latter method was already applied for 2 → 3 production mechanisms in the case of cc + jet [14] and even for 2 → 4 processes in the case of cccc [25], four-jet [26] and cc + 2jets [27] final states within the KaTie enviroment.
In general, the cross secton for pp → g(g)cc X reaction in the k T -factorization approach can be written as Then, the elementary cross section from the above can be written somewhat formally as: with n = 3 and n = 4 for g * g * → gcc and g * g * → ggcc, respectively, where E l and p l are energies and momenta of final state gluon(s) and charm quarks. Above only dependence of the matrix element on four-vectors of incident partons k 1 and k 2 is made explicit. In general, all four-momenta associated with partonic legs enter. Also in this case, the matrix element takes into account that both gluons entering the hard process are off-shell with virtualities k 2 1 = −k 2 1t and k 2 2 = −k 2 2t . Here, as a default choice we set the renormalization/factorization scale to be equal to the averaged sum of the transverse mass squared of the final state particles µ 2 = ∑ n i=1 m 2 iT n , where m iT = m 2 i + p 2 iT with n = 3 and n = 4 for 2 → 3 and 2 → 4 cases, respectively.
Calculating the minijets at tree-level requires some technical methods for regularization of the cross section. For this purpose, we follow the method adopted e.g. in PYTHIA [28] for the calculations of the 2 → 2 pQCD processes with light quarks and gluons in the final states. This procedue was also applied recently in the context of D-meson production via unfavoured fragmentation [29] or J/ψ-meson production in the color-evaporation model [30]. Here, we introduce a special suppression factor F sup (p T ) = p 4 T /(p 2 T0 + p 2 T ) 2 for the g * g * → gcc and g * g * → ggcc cross sections with p T being the outgoing minijet transverse momentum and with p T0 being a free parameter. This parameter could, in principle, be fitted to total charm cross section measured experimentally or calculated in the NLO/NNLO collinear calculations. Usually, the values p T0 = 1 − 3 GeV are taken in phenomenological applications. As a default set in our calculations here we use p T0 = 1 GeV. Within the referred range, we expect only a small sensitivity of the calculated charm quark distributions on the value of this parameter. The uncertainties could be visible only at very small transverse momenta of charm quark, i.e. p c T < 3 − 4 GeV. This is the region where still the leading-order mechanism should dominate. At larger charm quark transverse momenta, where higher-order contributions should play the most important role, our calculations are not sensitive to the choice of the regularization parameter. The calculated transverse momentum distributions of D-meson can be treated as the approximate attempt to study the differential charm cross section beyond the NLO.
Within the proposed scheme we sum together the three contributions g * g * → cc, gcc, and ggcc. It is known, when mixing different final state multiplicities that a problem of possible double counting appears. The double-counting effects shall also appear in the case under consideration. Their consistent treatment is not an easy task since there is a lack of well established theoretical techniques. However, very recently the consistent prescription for merging leading-and next-to-leading-order calculations for off-shell initial state partons were established [3]. Here, for the first time a similar procedure is adopted also for the case of NNLO corrections. According to this approach, we introduce for each of the three reactions, a set of the double-counting-exclusion (DCE) cuts. For the three reactions, transverse momenta of (mini)jets from the uPDF are constrained to be subleading. A similar constrain was also used in the PRA studies of dijet azimuthal decorrelations [8]. Therefore, the proposed DCE conditions set the following extra limitations on transverse momenta of incident partons: 1. k T < µ F for g * g * → cc, where µ F is the factorization scale, 2. k T < p T of the minijet for g * g * → gcc, 3. k T < p min T of the two minijets for g * g * → ggcc.
As was shown in Ref. [3] the above kinematic cuts shall provide a clear separation of events that correspond to the 2 → 2, 2 → 3 and 2 → 4 reactions. The first condition excludes the possible extra hard emissions from the uPDF in the 2 → 2 case, that are not under full theoretical/kinematical control. It reduces this contribution rather to the leading-order collinear calculations with c andc being the leading (mini)jets. The second condition assures that the hardest (mini)jet in the 2 → 3 event comes always from the hard matrix element. It removes the contributions that correspond rather to the mechanisms explicitly present in the 2 → 4 calculations. Similarly, the third condition assures that the two hardest (mini)jets in the 2 → 4 event also do not originate from the uPDF.
Including one or two hardest minijets from the uPDF in the 2 → 4 case in association with soft minijets from the matrix element may also lead to contributions already present in the case of the 2 → 3 processes. The additional hard emissions associated with charm determine the kinematics of the c andc and their correlations. In this context, having them at the level of hard matrix elements seems to be more accurate. Within the presented framework the leading-and higher-order contributions can be consistently taken into account together without additional double-counting subtractions to describe e.g. correlation observables.

Kimber-Martin-Ryskin uPDF
As a default, in the calculations below we use the leading-order Kimber-Martin-Ryskin (KMR) approach [32,33] with the angular (or rapidity) ordering constraints imposed, which comes from inclusion of coherence effects in gluon emission. 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.
The virtual (loop) contributions may be resummed to all orders by the Sudakov form 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 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 incoming gluons. Here, the variable ∆ = k t /(k t + µ) introduces a restriction of the phase space for gluon emission due to the angular-ordering condition. This constraints translate into the permission for hard emissions from the uPDF, that correspond to the k t > µ kinematical regime.
Taking all together, the precise expression for the unintegrated gluon distribution This prescription was found to be consistent with the Multi-Regge-Kinematics limit of the QCD amplitudes [3].
In numerical calculations below we apply different sets of the KMR gluon uPDF, obtained from different collinear PDFs, including LO, NLO, and even NNLO fits. As was discussed in Ref. [34], the LO KMR model together with NLO PDFs leads to gluon distributions compatible with their counterparts calculated within full NLO KMR approach.
Thus, in phenomenological studies one can safely neglect effects related to the higherorder perturbative-splitting functions and concentrate only on the collinear PDF input.

Parton-Branching uPDF
The Parton Branching (PB) method, introduced in Refs. [31,35], 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 [31] (2.10) These equations can be solved by an iterative Monte Carlo method. In this method every resolvable branching is reconstructed explicitly and the full kinematics at each branching is taken into account. 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.
There are two sets available 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 [35]. Both of them are based on the HERAPDF2.0 collinear parton densities at NLO. In the numerical calculations below we use the PB-NLO-2018-set1 uPDF.

Comparison of the uPDF distributions
Let us present now a numerical comparison between the two uPDF models used in the calculations below. Here, we compare the KMR-MMHT2014lo and PB-NLO-set1 gluon unintegrated distributions. In Fig. 3 we show transverse momentum dependence of the uPDFs for two different values of longitudinal momentum fractions: x = 0.01 (upper panels) and x = 0.0001 (lower panels) as well as for three different values of the factorization scale: µ = 3, 10, 100 GeV (left, middle and right panels). We observe significant differences between the two models of uPDF at both, very small and very large transverse momenta of gluons. The differences for k t 1 GeV come from miscellaneous treatment of the non-perturbative regime in the considered uPDFs, which in both cases is rather uncertain. However, the main visible difference appears in the region of large transverse momenta. The KMR uPDF (solid lines) has long tails which is a consequence of k t > µ contributions allowed for the gluon emissions. In contrast, the PB-NLO-set1 uPDF (dashed lines) is strongly suppressed in this kinematical regime. As will be discussed in the following, this behaviour of the two uPDFs has a crucial meaning for valueable predictions of charm hadroproduction at the LHC.

A. The standard k T -factorization calculations of the D-meson cross section including only
g * g * → cc mechanism with the KMR uPDFs We start presentation of our numerical results with the inclusive distributions of Dmeson and with correlation observables for DD meson-antimeson pair production. We compare our theoretical predictions with the LHCb open charm data from pp-scattering [37,38]. Here, we follow the standard k T -factorization approach and calculate the cross section for cc-pair production by taking into account the g * g * → cc mechanism. We use the KMR gluon uPDF in the original form, that allows for extra hidden hard emissions at the very last step of its evolution, i.e. including contributions from the k T > µ F kinematical regime. In this way, a part of real higher-order corrections is effectively taken into account in the calculations. This was already discussed in the case of charm production, e.g. in Ref.
[1]. Here, we wish to extend those studies by discussing some important details of the calculation, relevant to estimate overall uncertainties of the model. TeV together with the LHCb data [37]. Here, we compare the results for different collinear PDFs used for calculating the KMR uPDFs.
In Fig. 4 we show the transverse momentum distributions for different rapidity bins (left panel) and the rapidity distribution (right panel) of the charged D-meson measured by the LHCb experiment in the kinematical range: 0 < p T < 8 GeV and 2 < y < 4.5.
We compare results obtained with two different sets of the KMR gluon uPDF: KMR-MMHT2014lo (solid) and KMR-CT14lo (dotted). Both of them lead to a reasonable description of the LHCb data for larger transverse momenta, however, we observe some vis- Here, instead of varying the default set to produce the scale uncertainty band, we also consider two different sets: iT n . A comparison of corresponding results for different scales is shown in Fig. 6. These three sets of the scales lead to visible differences only for very small meson p T 's. The overall uncertainty related to the scales is of the same order as those discussed above.    and rapidity distance (bottom-right) distributions for charged DD meson-antimeson pair production at √ s = 7 TeV together with the LHCb data [38]. Here, we compare results for different collinear gluon PDFs used in calculating the KMR uPDF. panels) and rapidity difference ∆Y = |y D − yD| (bottom right panels). Again, we show uncertainties due to the choice of collinear gluon PDFs (Fig. 7) and due to the choice of scales (Fig. 8). Here the LHCb correlation data [38] are not absolutely normalized and Summarizing this subsection, we conclude that within the typical pQCD uncertainties we are able to get a satisfactory description of the LHCb charm data. The statement is valid for both, the absolutely normalized inclusive D-meson distributions as well as for the shapes of the DD correlation observables. The framework of the k T -factorization together with the KMR gluon uPDF allows to describe the LHCb charm data already within the leading-order g * g * → cc mechanism. This is completely opposite to the calculations within the collinear-approximation. There, only the NLO framework is able to obtain the same level of quality of the description of the LHC heavy flavour data [4,5]. This clearly shows that within the k T -factorization approach we effectively include higher-order contributions. However, the fact and the size of the effective resummation strictly depends on the construction of the used uPDF. The KMR model is unique and very useful in this context. It allows even for two extra emissions of hard partons from the uPDFs, that cor-respond to the gg → gcc and the gg → ggcc contributions. As we have shown for charm production at the LHC this model seems to work very well, however, the overall picture is more complicated.
The role of the extra emissions from the KMR uPDF can dramatically change when going to the lower energies. The emissions are not under full kinematical control and in the case of some processes, e.g. for dijets production, they may even lead to a problematic double counting. Moreover, other models of the uPDFs from the literature do not contain such large contributions from the k T > µ F regime and are useless at leading-order calculations for processes where higher-orders are of the special importance. Due to the lack of the full NLO/NNLO formalism with off-shell initial state partons we propose a simplified scheme for the calculation of heavy flavour cross section in the k T -factorization with higher-order mechanisms taken into account at the tree level.

The Kimber-Martin-Ryskin uPDF with limited hard emissions
The idea of the proposed scheme is to exclude the extra hard emissions from the uPDF and include the higher-order contributions g * g * → gcc and g * g * → ggcc explicitly at the level of hard matrix elements.
First of all we wish to show the importance of the k T > µ F contributions in the case of the KMR uPDF for the leading-order g * g * → cc mechanism. In Fig. 9 we present the c-quark transverse momentum (left) and ϕ cc azimuthal angle (right) distributions for ccpair production at √ s = 7 TeV. The solid histograms correspond to the standard KMR calculations with the k T > µ F limitation included and the dashed histograms are for the calculations with excluded contributions from the k T > µ F region. We observe a significant differences between the both results. The k T > µ F contribution is very important for the whole considered distribution of the transverse momenta and concentrated especially at small azimuthal angles. The k T < µ F limitation of the KMR uPDF that allows only for soft extra emissions and as a consequence significantly reduces the predicted cross section. and g * g * → ggcc (right) mechanisms.
In Fig. 10 we present a similar analysis as the above one for the higher-order components. Here we consider the role of the k T > µ F contribution in the KMR uPDF both for the g * g * → gcc (left panel) and g * g * → ggcc (right panel) mechanism. Also in the case of the higher-order processes the k T > µ F kinematical region in the KMR uPDF significantly contributes to the charm quark production cross sections in the whole range of considered p T 's. As we already argued, in the case of higher-order processes the k T < µ F limitation is not enough to fully avoid double-counting effects when summing up the leading-and higher-order contributions. Therefore, we also plot here the contributions that correspond to the case of the proposed double-counting-exclusion (DCE) cuts (see the dotted histograms). The effects related to these cuts are very important for both, the 2 → 3 and 2 → 4 processes. In the latter case the DCE cut significantly reduces the basic cross section by about one order of magnitude. results for summed contributions from the g * g * → cc, g * g * → gcc and g * g * → ggcc mechanisms with extra conditions. Now we wish to present the results of the proposed new scheme that includes higherorder corrections explicitly in comparison to the standard (leading-order) KMR calculations. In Fig. 11 we show the standard 2 → 2 KMR calculations with the k T > µ F included (dashed histograms) and the results obtained within the proposed scheme for 2 → 2 + 3 + 4 calculations. For the latter case, here we show both results, with only the k T < µ F limitations (dotted histograms) and with the DCE cuts (solid histograms). We clearly see that the DCE cuts are necessary for 2 → 2 + 3 + 4 calculations to reproduce the successful standard 2 → 2 KMR calculations. The calculations with the k T < µ F limitations would also lead to a significant overestimation of the LHC charm data. The calculations within the 2 → 2 + 3 + 4 scheme with the DCE cuts almost coincides with the standard 2 → 2 calculations in the broad range of the considered transverse momenta of charm quark. Some discrepancy appears only at small p T 's. The reason could be a different collinear PDFs used in both calculations. It is not clear for the 2 → 2 case whether the LO, NLO or NNLO PDFs should be used, so there we keep the LO PDF as a default while in the case of the 2 → 2 + 3 + 4 we assume that the NNLO PDFs are the most appropriate. Our new scheme also leads to a very similar azimuthal angle distribution as in the standard 2 → 2 calculations.

The Parton-Branching uPDF
In this subsection we basically repeat the above studies for the KMR uPDF but here we apply the Parton-Branching uPDFs [31,35]. As we observe from Figs. 12 and 13 in the case of the PB uPDFs the k T > µ F contributions are very small for the 2 → 2 mechanism and almost negligible for the 2 → 3 and 2 → 4 higher-orders (see almost coinciding solid and dashed histograms in Fig. 13). Therefore, it is imposible to decribe the LHCb open charm data within the PB uPDFs when considering only the g * g * → cc mechanism. Here, the effects of the DCE cuts are still sizeable, however, much smaller than in the case of the KMR uPDF. In Fig. 14, for a more general comparison, we show the 2 → 2 + 3 + 4 results with the DCE cuts for the two uPDFs together on the same plots. In the case of the 2 → 2 mechanism some differences bewteen the two results are obtained. For the 2 → 3 and 2 → 4 mechanisms, when the DCE cuts are imposed the results almost coincide.
In Fig. 15 we again show the standard 2 → 2 KMR calculations with the k T > µ F included (dashed histograms) and the results obtained within the proposed scheme for the 2 → 2 + 3 + 4 calculations but with the KMR (dashed histograms) and the PB uPDFs (solid histograms) in addition. One can observe that the 2 → 2 + 3 + 4 calculations for the two different uPDFs lead to very similar results. The proposed procedure is the only scheme which allows for a reasonable prediction for heavy flavour production within the PB uPDFs. Further improvement can be done only by the full NLO/NNLO k T -factorization calculations. with the extra double-counting-exclusion cuts for the g * g * → gcc (left) and g * g * → ggcc (right) mechanisms.

Numerical representation of the double-counting exclusion cuts
Here, we wish to illustrate whether the proposed double-counting exclusion cuts can really result in separation of the leading-and higher-order contributions. In Fig. 16 we plot the two-dimensional distributions as a function of the leading gluon jet p T and averaged transverse momentum of the charm quark and antiquark in the 2 → 2 (left panels), 2 → 3 (middle panels) and 2 → 4 (right panels) events that could be helpful in schematic illustration of the complementarity of phase spaces for the leading-and higher-order contributions. The top panels correspond to the direct calculations without the DCE cuts while the bottom panels correspond to the calculations with the DCE cuts included. As we observe, the DCE cuts remove from the 2 → 2 calculations the contributions of the 2 → 3 and 2 → 4 type, as well as suppress the 2 → 3 and 2 → 4 components in the region populated by the 2 → 2 mechanism. Although, the separation is not sharp which may be related to the chosen scales, the main tendency of the applied procedure is clear and seems to support the applied procedure. Similar conclusions were drawn in Ref. [2] in the case of bb-pair production. results for summed contributions from g * g * → cc, g * g * → gcc and g * g * → ggcc mechanisms with the extra conditions.

Double-counting and two-dimensional distributions
We wish to discuss also how our prescription devoted to avoid double counting works for example for a two-dimensional distribution in (M cc , φ cc ). In Fig. 17  What about higher orders? In Fig. 18 we show similar distributions for 2 → 3 (upper panels) and 2 → 4 (bottom panels). We start with the PB-NLO-set1 distributions (left and middle panels). We see that the removed (for the KMR) regions reappear in the higherorder corrections. The left panels are results of direct calculation, whereas the middle panels include the DCE cuts. The direct calculations (left panels) lead to a significant contributions for back-to-back configurations already included in the 2 → 2 processes.
We observe that our DCE cuts allow to avoid double counting. In the right panels we show for comparison results with the DCE cuts but for the KMR uPDF.
We see that our DCE cut fulfills the necessary requirements supporting its practical In the latter case a kind of a separation of the leading-and higher-order contributions is obtained.

Comparison to the calculations based on the collinear approach
Here, we wish to demonstrate how the 2 → 2 + 3 calculations at the tree-level with the DCE cuts differs from the full NLO approach. This can be done only in the collinear approximation because of the lack of the NLO k T -factorization framework. In the left panel of Fig. 19 we compare results for the transverse momentum distribution of charm quark obtained in the full NLO framework (solid line) and calculated by summing the 2 → 2 and 2 → 3 contributions at the tree-level (solid histogram). The two approaches almost coincide in the broad range of considered p T 's. Significant differences appear only at very small transverse momenta where the effects related to the loop-corrections are expected to be of special importance. Here, we plot in addition the 2 → 2 + 3 + 4 contribution (dotted histogram). We observe a huge contribution of the NNLO-type to the transverse momentum distribution of charm quarks. It is not taken into account in the state-of-art calculations of the FONLL and GM-VFNS frameworks. Here, within our more pragmatic model we only wish to pay attention to the importance of the NNLO corrections to differential distributions of heavy quarks. Definite conclusions about their size are strongly limited since the NNLO collinear framework is not available for differential distributions of charm and bottom quarks.
For completeness, in the right panel of Fig. 19 we compare the 2 → 2 + 3 + 4 contributions calculated in both, the collinear and k T -factorization tree-level approach with the DCE cuts. In both cases we get the description of the experimental data of the same quality as in the case of the standard 2 → 2 k T -factorization calculations with the KMR uPDF.  [37]. Here, we show the PB-NLO-set1 and the KMR-MMHT2014nnlo results for summed contributions of g * g * → cc, g * g * → gcc and g * g * → ggcc mechanisms with the extra conditions.  [38]. Here, we show results for the PB-NLO-set1 and KMR-MMHT2014nnlo results for summed contributions from g * g * → cc, g * g * → gcc and g * g * → ggcc mechanisms with the extra conditions.

IV. CONCLUSIONS
In the present paper we have considered charm production at the LHC within the k T -factorization approach beyond the standard leading-order g * g * → cc partonic mechanism. For the first time we have included in this context next-to-and next-to-next-toleading order mechanisms for the differential distributions. We have proposed a new scheme for calculating the charm quark cross section including in addition the 2 → 3 and 2 → 4 higher-order contributions at the tree-level. The calculations of the g * g * → gcc and g * g * → ggcc mechanisms have been done also in the framework of the k T -factorization, with off-shell initial state partons, for two different unintegrated gluon densities from the literature -Kimber-Martin-Ryskin and recent Parton Branching uPDFs. We have proposed special conditions in order to avoid the problem of double counting when calculating the higher-order corrections at the tree-level.
We have made a detailed comparison of the results for charm production obtained in the standard 2 → 2 calculations with the KMR uPDF and those from the proposed 2 → 2 + 3 + 4 scheme with the higher-order contributions taken into account at the treelevel. Both approaches were found to lead to very similar results. This conclusion applies exclusively for the KMR uPDF model. The analogous analysis have been done also for the PB-NLO-set1 uPDFs. In the latter case, the 2 → 2 calculations leads to a significant underestimation of the charm cross section at the LHC. Within this model of the gluon uPDF the experimental data can be described only in the 2 → 2 + 3 + 4 scheme, with higher-order contributions taken into account at the level of hard-matrix elements. This observations may be also valid for other models of the uPDFs from the literature, including different CCFM-fits, that do not allow for extra hard emissions encoded in their evolution.
Several differential distributions, including correlations observables, for open charm mesons for the LHCb experiment have been analyzed. Within the proposed 2 → 2 + 3 + 4 calculational scheme a good quality description of the data have been obtained for both the KMR and the Parton-Branching unintegrated gluon densities.