Transverse momentum dependent parton densities in processes with heavy quark generations

To study the heavy quark production processes, we use the transverse momentum dependent (TMD, or unintegrated) gluon distribution function in a proton obtained recently using the Kimber-Martin-Ryskin prescription from the Bessel-inspired behavior of parton densities at small Bjorken $x$ values. We obtained a good agreement of our results with the latest HERA experimental data for reduced cross sections $\sigma^{c\overline{c}}_{\rm red}(x,Q^2)$ and $\sigma^{b\overline{b}}_{\rm red}(x,Q^2)$, and also for deep inelastic structure functions $F_2^c(x,Q^2)$ and $F_2^b(x,Q^2)$ in a wide range of $x$ and $Q^2$ values. Comparisons with the predictions based on the Ciafaloni-Catani-Fiorani-Marchesini evolution equation and with the results of conventional pQCD calculations performed at first three orders of perturbative expansion are presented.


Introduction
Recently, important new data have appeared on the cross sections for the open charm and beauty production in neutral current deep inelastic electron-proton scattering (DIS) obtained [1] by combining the results of H1 and ZEUS Collaborations at HERA. Measurements have shown that heavy flavour production in DIS proceeds predominantly via the photon-gluon fusion process, γg → QQ, where Q is the heavy quark. The cross section therefore depends strongly on the gluon distribution in the proton and heavy quark mass. Moreover, an analysis of the data in the framework of perturbative Quantum Chromodynamics (QCD) has been done [1], where the massive fixed-flavour-number scheme and different implementations of the variable-flavour-number scheme were used.
The theoretical description of the heavy quark production processes can also be performed with transverse momentum dependent (TMD), or non-integrated, functions of the density of partons (quarks and/or gluons) in a proton [2,3]. These quantities, depending on the fraction of the longitudinal momentum x of the proton carried by the parton, the two-dimensional transverse momentum of the parton k 2 T , and the hard scale µ 2 of a complex process, contain nonperturbative information about the proton structure, including the transverse momentum. The TMD factorization theorems provide the necessary basis for separating hard parton physics, described in terms of perturbative QCD, and soft hadron physics. Currently, there are a number of factorization approaches that include the dependence of the parton distribution fnctions (PDFs) on the transverse momentum: for example, the Collins-Soper-Sterman [4] approach, developed for semi-inclusive processes with a finite and nonzero ratio between the rigid scale µ 2 and total energy s, as well as the approach to high-energy factorization [5,6] (or k T -factorization [7]), valid at a fixed limit of the hard scale and at high energies.
With high-energy factorization, the TMD density of gluons satisfies the Balitsky-Fadin-Kuraev-Lipatov (BFKL) [8] or Ciafaloni-Catani-Fiorani-Marchesini (CCFM) [9] evolution equations, which resum the contributions of large logarithms proportional to α n s ln n s ∼ α n s ln n 1/x. These terms are important at high energies s (or, equivalently, at low x values). Thus, high-order radiative corrections can be effectively taken into account in the estimated cross sections (namely, the part of the next-to-leading order (NLO) + of the next-to-leading order (NNLO) + ... terms corresponding to the emission of the original gluons). Phenomenological applications of the high-energy factorization approach augmented by the CCFM are well known in the literature (see, for example, [10][11][12][13][14][15][16][17][18][19][20] and references therein).
In addition to the CCFM equation, there are also other approaches to determining the TMD PDFs in a proton. They can be estimated using schemes based on the usual Dokshitser-Gribov-Lipatov-Altarelli-Parisi (DGLAP) [21] equations, namely the parton branching approach (PB) [22] and the Kimber-Martin-Ryskin (KMR) recipe [23]. The first of them gives a numerical iterative solution of the DGLAP evolution equations for collinear and TMD PDFs using the concept of resolvable and non-resolvable branching and the application of Sudakov's formalism to describe the evolution of a parton from one scale to another without decidable branching. The splitting kinematics at each branching vertex is described by the DGLAP equations and the angular ordering condition for parton emission, which can be used instead of the usual DGLAP ordering by virtuality. The KMR approach is a formalism invented to construct TMD PDFs from well-known traditional (collinear) PDFs under the key assumption that the dependence of parton distributions on transverse momentum comes only in the last stage of evolution. It is believed that the KMR procedure effectively takes into account a main part of the following nextto-leading logarithmic (NLL) terms α s (α s ln µ 2 ) n−1 compared to the leading logarithmic approximation (LLA), where terms are proportional to α n s ln n µ 2 taken into account only. The KMR approach is currently being explored in NLO [24] and is commonly used in phenomenological applications (see, for example, [11,[13][14][15][16][17][18][19] and references therein), where the usual PDFs (for example, the NNPDF [25] or CTEQ [26] ones) accepted numerically as input. The relationship between PB and KMR scenarios has been discussed [27], and the relationship between PB and CCFM approaches has been recently created [28].
The KMR formalism was used in our recent article [29] for analytical calculations of the TMD PDFs. These calculations are based on the [30][31][32] expressions for usual PDFs obtained with the generalized double asymptotic scaling (DAS) approach [31][32][33]. This scaling, related to the asymptotic behavior of DGLAP evolution, was discovered many years ago by [34]. It was shown [30][31][32] that the flat initial conditions for the DGLAP equations used in the generalized DAS scheme lead to the Bessel-like behavior of PDFs at small x values. Using the results [30][31][32], analytical expressions for the TMD quark and gluon densities has be obtained [29] in the leading order (LO) 1 . In Ref. [29], we have implemented various kinematic constraints that exist in the KMR recipe (namely, angular and strict ordering conditions), and investigated the relationship between the differential and integral formulations of the KMR procedure recently mentioned in [36].
In the present paper we continue our investigations and analyze the combined H1 and ZEUS experimental data [1] for the (reduced) charm and beauty cross sections σ cc red and σ bb red and charm and beauty contributions to the proton structure functions (SFs) F 2 (x, Q 2 ) and F L (x, Q 2 ) [37][38][39] and their ratio for different Q 2 values. Studying the earlier data on the charm SF F c 2 in the proton from H1 [40] and ZEUS [41] Collaborations at HERA, obtained for x ∼ 10 −4 , it was found that the charm contribution to the total proton SF F 2 is about 25%, which is significantly larger than that found by the European Muon Collaboration at CERN [42] for large x, where SF F c 2 was only 1% of F 2 . Such a large value of F c 2 required extensive experimental and theoretical studies of heavy quark production processes (see, for example, the data [1] studied in this paper, as well as the experimental data [43] from LHCb Coolaboration at CERN for the prompt charm production in pp collisions). Theoretical studies usually serve to confirm that HERA and LHC data can be described by perturbative charm generation within QCD (see, for example, reviews [44][45][46] and references therein). We also note here that, historically, it was in the study of these processes that the k T -factorization was introduced and tested (see [5][6][7]).
The production of charmed mesons ah hadron colliders is dominated by the gg → cc subprocess, and therefore it provides a sensitive probe of the gluon density at small x. In particular the data [43] of LHCb Collaboration provides an information on the gluon at values of x as small as x ∼ 10 −6 (see [47] and discussions therein). This very small-x region is also crucial for the calculations of signal and backgroud processes for ultra-high energy neutrino astrophysics (see [48] and [49] for calculations of the high energy neutrino crosssection and the prompt atmospheric neutrino flux, respectively). A study of the heavy quark production will be contitued at future lepton-hadron and hadron-hadron colliders, such as LHeC, FCC-eh and FCC-hh, respectively (for a review, see [50,51] and references therein).
To study the process of the heavy quark production, we produce the k T -factorization predictions in two ways, namely, using the framework of DAS approach and CCFM evolution equation. The direct comparison of these predictions is interesting and could be rather useful to discriminate between the different approaches to evaluate the TMD parton (mainly gluon) density in a proton. Then, we calculate the high-energy asymptotics of the heavy quark parts of the SFs F 2 and F L at first three orders of perturbation theory and present the numerical comparison of these higher-order predictions with corresponding results of the k T -factorization calculations.
The outline of our paper is following. In Sections 2 and 3 we briefly describe our theoretical input. Section 4 presents our numerical results for the reduced charm and beauty cross sections and charm and beauty parts of SFs F 2 (x, Q 2 ) and F L (x, Q 2 ) in a wide Q 2 range. Section 5 contains our conclusions. In Appendix A we present the high energy asymptotics of the heavy quark contribution to the SFs F 2 and F L at first three orders of perturbation theory. Appendix B contains the simple approximations of these formulas for the ratio of the heavy quark parts of the SFs F 2 and F L , which could be useful for subsequent applications.

2
k T -dependent Wilson coefficient functions The differential cross section σ QQ (hereafter Q = c, b) can be presented in the simple form: where F Q k (x, Q 2 ) (hereafter k = 2, L) are heavy quark parts of the proton SFs F k (x, Q 2 ), x and y are the usual Bjorken scaling variables. Here we present the basic elements of the relations between SFs F Q 2 (x, Q 2 ) and F Q L (x, Q 2 ) and TMD PDFs. More detail analysis can be found in [52].
In the k T -factorization approach, the SFs F Q k (x, Q 2 ) are driven at small x by gluons and related in the following way to the TMD gluon distribution f g (x, k 2 T , µ 2 ): The functions C k,g (x, Q 2 , m 2 f , k 2 ⊥ , µ 2 ) may be regarded as the structure functions of the off-shell gluons with virtuality k 2 ⊥ (hereafter we call them as Wilson coefficient functions). Following [52], we do not use the Sudakov decomposition, which is sometimes quite convenient in high-energy calculations. Here we only note that the property k 2 = −k 2 ⊥ comes from the fact that the Bjorken parton variable x in the standard and in the Sudakov approaches coincide.
The k ⊥ -dependent Wilson coefficient functions have the following form: where C k,g (x) andC k,g (x) corresponds to application of the Feynman P αβ F polarization tensor and additional tensor of gluon polarization P αβ a : and we omitted the dependence of the coefficient functions on the heavy quark mass m, Q 2 , k 2 ⊥ and hard scale µ 2 . The sum of the polarizations tensors produces so-called BFKL polarizations tensor:P Indeed, as it was shown in [52] P αβ After applying the photon polarizations tensorŝ which can be rewritten as the results for the coefficient functions C k,2 (x) have been calculated in [52]: where a s = α s /(4π) is the strong coupling constant, Θ(x 1 − x) is the Heaviside step function and Following to the representation (7) of gluon polarization tensors, we can represent the basic functions f (i) (1, 2) in the following form: where Here and

The case of on-shell gluons
In the particular case of off-shell initial gluons, when k 2 = 0, we have (see [52] for more details): with and where Here [53]. Indeed, we have which are in agreement with the old results [53].

Kimber-Martin-Ryskin approach
Here we present the main elements of TMD PDFs, based on the KMR prescription in so-called interal formulation (see [36]) and DAS appoach for usual PDFs. More details, including the differential formulation of the KMR prescription, can be found in our previous paper [29].

Sudakov form factors
The Sudakov form factor T a (µ 2 , k 2 ) have the following form (see (2.4) in [36]): When ∆ is a constant, we have where

Conventional PDFs
At LO, the conventional sea quark and gluon densities f a (x, µ 2 ) can be written as follows where for the color SU (N c ) group, I ν (σ) andĨ ν (σ) (ν = 0, 1) are the combinations of modified Bessel functions (at s ≥ 0, i.e. µ 2 ≥ Q 2 0 ) and usual Bessel functions (at s < 0, i.e. µ 2 < Q 2 0 ): where I 0 (σ) =Ĩ 0 (σ) and andd are the singular and regular parts of the anomalous dimensions, β 0 is the first coefficient of the QCD β-function in the MS-scheme. The results for the parameters A a and Q 2 0 can be found in [30,55]; they were obtained for α s (M Z ) = 0.1168.
It is convenient to show the following expressions:

TMDs in KMR approach
Now we can use (26) to find the results for TMDs without derivatives. After some algebra we have We can obtain (see more details in [29]): and 3.4 Other prescriptions 1. For the phenomenological applications, we use the cut-off parameter ∆ in the angular ordering [36] (the case of strong ordering can be found in [29]): In all above cases, except the results for T a (µ 2 , k 2 ), we can simply replace the parameter ∆ by ∆ ang . For the Sudakov form factors, we note that the parameters ∆ contribute to the integrand in (28) and, thus, their momentum dependence changes the results in (29). To perform the correct evaluation of the integral (28), we should recalculate the p 2 integration in (28). So, we have The analytic evaluation of T (ang) a (µ 2 , k 2 ) is a very cumbersome procedure, which will be accomplished in the future. With the purpose of simplifying our analysis, below we use the numerical results for T (ang) a (µ 2 , k 2 ).

2.
As it was shown [32], the results of the fits of the experimental data for SF F 2 are not very well at low Q 2 values. To overcome this problem, following to [30], it is possible to use a modification of the strong-coupling constant in the infrared region. Specifically, usually it is possible to consider two modifications: the "frozen" coupling constant (see, for example, [30,56]) and analytic one [57,58], which effectively increase the argument of the strong coupling constant at small µ 2 values, in accordance with [59,60]. As one can see [30,55,61], the fits based on the "frozen" and analytic strong coupling constants are rather similar and describe the F 2 (x, Q 2 ) data in the small Q 2 range significantly better than the canonical fit. However, as it was shown in [29], the "frozen" coupling constant leads to a better agreements with data sets. Thus, we will use it in our present analysis. So, we introduce a freezing of the strong coupling constant by changing its argument as µ 2 → µ 2 + M 2 ρ , where M ρ is the ρ meson mass [56]. Thus, in the formulae of Section 3 we introduce the following replacement 3. In the phenomenological applications (see Section 4) the calculated TMD parton densities will be used to predict the reduced cross sections σ QQ red and proton SFs F Q k (x, Q 2 ). According to k T -factorization theorem [5,7], the theoretical predictions for these observables can be obtained by convolution (2) of the TMD gluon densities and corresponding off-shell production amplitudes. So, we need the TMD quark and gluon distributions in rather broad range of the x variable, i.e. beyond the standart low x range (x ≤ 0.05).
It was shown (see [29] and discussions therein) that the analytic expressions for TMD parton densities can be modified in the form: that in agreement with a similar modifications of conventional PDFs (see, for example, the recent paper [62], where similar extension has been done in the case of EMC effect from the study of shadowing [63] at low x to antishadowing effect at x ∼ 0.1 − 0.2). The value of β a (0) can be estimated from the quark counting rules [65]: where the symbol v marks the valence part of quark density. Usually the β v (0), β g (0), β q (0) are determined from fits of experimental data (see, for example, [66,67]). In our analysis, we use the numerical values of β g (0) = 3.03 which have been extracted [29] from the fit to the inclusive b-jet production data taken by the CMS [68] and ATLAS [69] Collaborations in pp collisions at √ s = 7 TeV.

Phenomenological applications
We are now in a position to apply the TMD parton densities, obtained in [29] and shown above, for phenomenological applications. In the present paper we consider the reduced charm and beauty cross sections σ cc red and σ bb red and charm and beauty contributions to the deep inelastic proton SFs F 2 (x, Q 2 ), which are directly related with the gluon content of the proton. These observables were measured in ep collisions at HERA with rather good accuracy (see [1] and [37][38][39], respectively.) Additionally, in the evaluations below we will use latest TMD gluon density in a proton, obtained from the numerical solution of the CCFM evolution equation, namely, JH'2013 set 2 one [70]. Our choice is motivated mainly by the fact that the CCFM equation smoothly interpolates between the small-x BFKL gluon dynamics and conventional DGLAP one, as it was mentioned above. The input parameters of starting (initial) gluon distribution implemented into the JH'2013 set 2 were fitted to describe the high-precision DIS data on structure functions F 2 (x, Q 2 ) and F c 2 (x, Q 2 ) at x ≤ 5 · 10 −3 (see [70] for more information). Everywhere below, we set the charm and beauty masses to be equal to m c = 1.65 GeV and m b = 4.78 GeV [72]. We use the one-loop formula for the QCD coupling α s with n f = 4 quark flavours at Λ QCD = 143 MeV (that corresponds to α s (m 2 Z ) = 0.1168) for for analytically calculated TMD gluon density as described above. In the case of CCFM-evolved gluon, we apply the two-loop expression for α s with n f = 4 and Λ QCD = 200 MeV, as it is dictated by the fit [70].

Reduced cross sections σ QQ
red and SFs F Q k (x, Q 2 ) Usually the differential cross section of heavy quark production in deep inelastic scattering (1) are represented in terms of reduced cross sections σ QQ red , which are defined as follows: Hence, taking together expressions (1) and (45), the reduced cross section σ QQ red can be easily rewritten through F Q 2 (x, Q 2 ) and F Q L (x, Q 2 ) as where the ratio R Q (x, Q 2 ) in defined as: The evaluation below is based on the formulas (46) and (2) with the coefficient functions as given by (10) - (19). Our numerical results for reduced cross sections σ cc red and σ bb red are shown in Figs. 1 and 2, respectively, in comparison with the latest H1 and ZEUS data [1]. The shaded bands represent the theoretical uncertainties of our calculations. We find that the k Tfactorization predictions obtained using derived analytical expressions for TMD gluon density in a proton are in perfect agreement with the HERA data in a wide region of x and Q 2 within the theoretical and experimental uncertainties, both in normalization and shape. These results tend to slightly overshoot the JH'2013 set 2 predictions in the region of small x and especially at low Q 2 . At larger Q 2 and/or moderate or large x ≥ 10 −2 the CCFM-evolved gluon density tends to overestimate the HERA data, that could be due to the determination of corresponding input parameters at small x only (see [70]). To estimate the scale uncertainties of our calculations, the standard variations (by a factor of 2) in default renormalization and factorization scales, which were set to be equal to µ 2 R = 4m 2 Q + Q 2 and µ 2 F = Q 2 , respectively, were introduced. To show the contribution of the longitudinal structure functions F c L (x, Q 2 ) and F b L (x, Q 2 ), we present also in Figs. 1 and 2 the results for F c 2 (x, Q 2 ) and F b 2 (x, Q 2 ) as dotted curves. The difference between the estimated σ QQ red and F Q 2 (x, Q 2 ) is coming from the contribution of the longitudinal SFs F Q L (x, Q 2 ), as it can been clearly seen from (46). So, our calculations show that these contributions are rather important at low x.
To show the difference between σ QQ red and F Q 2 (x, Q 2 ) more clearly, in Figs. 3 and 4 we present our results for SFs F c 2 (x, Q 2 ) and F b 2 (x, Q 2 ) in comparison with the latest ZEUS [37] and H1 [38,39] data. Our corresponding predictions for the reduced cross sections σ cc red and σ bb red are presented here as dotted curves. One can see again that the results obtained with analytically evaluated TMD gluon density are in good agreement with the latest HERA data for both structure functions F c 2 (x, Q 2 ) and F b 2 (x, Q 2 ) in a wide region of x and Q 2 . The CCFM-evolved gluon JH'2013 set 2 provides a bit worse description of the HERA data, although these results are rather close to the measurements. We find that the discrepancy between two considered approaches tends to be more clearly pronounced at large Q 2 .

Ratio
Following the results of [52] and using our coefficient functions obtained in Section 2 and TMD gluon density presented in Section 3, now we can produce predictions for the ratio R Q (x, Q 2 ) according to (47). Results of our calculations for R c (x, Q 2 ) are presented in Fig. 5, where we plot this ratio as a function of x in a wide Q 2 range. As earlier, we have applied two TMD gluon densities in a proton discussed above 2 .
Our calculations leads to more or less flat (independent on x) behaviour of R c (x, Q 2 ) with R c ∼ 0.1 at low Q 2 ∼ 5 GeV 2 and R c ∼ 0.3 − 0.4 at higher Q 2 ∼ 200 GeV 2 values.
The results obtained with our TMD gluon and CCFM-evolved one are in a good agreement to each other. The difference between them is visible at very large Q 2 only. Moreover, the obtained predictions are in good agreement with ones [52], which were based in-turn on rather old representations for TMD gluon density (see [73] and more recent [3]). As a next point of our study, we would like to compare the results for the ratio R c (x, Q 2 ), obtained in k T -factorization with the oneR c (x, Q 2 ) (see A4), where the ratiô R c (x, Q 2 ) was obtained in the conventional (collinear) QCD factorization at first three orders of perturbation theory (see Appendix A) represented by the solid, dashed and dotted grey curves in Fig. 5. To evaluate the latter, we have used the LO DAS parton densities presented in Section 3.2. Note that the results for this ratio, obtained using the power-like behaviour x −∆ of the collinear PDFs, can be found in [74]. Our calculations show that the k T -factorization predictions rather close to the results obtained beyond LO of collinear perturbation theory. This is in complete agreement with the usual statement about the property of k T -factorization, which resums the main part of higher order pQCD contributions at small x. Indeed, the LO results obtained in the collinear perturbation theory lead to too small values 3 for the ratioR c (x, Q 2 ). However, the Q 2 -dependence of the ratioR c (x, Q 2 ) in the NLO and NNLO results are noticebly different from the corresponding Q 2 -dependence evaluated with the TMD gluons. In fact, in the k T -factorization approach the ratio R c (x, Q 2 ) grows fastly when Q 2 increased whereas in collinear perturbation theory the ratioR c (x, Q 2 ) grows slowly. Of course, the observed difference between the predicted R c (x, Q 2 ) andR c (x, Q 2 ) ratios at moderate and large Q 2 is unclear and needs an explanation, especially because there are no experimental data for the SF F c L (x, Q 2 ) and, accordingly, for the ratio R c (x, Q 2 ).
Indeed, the k T -factorization with the estimated R c (x, Q 2 ) leads to good agreement between experimental data and theoretical predictions for both reduced cross sections σ cc red and SF F c 2 (x, Q 2 ), as one can see in Figs. 1 and 3. From another side, the experimental data for both σ cc red and F c 2 (x, Q 2 ) are in good agreement with the corresponding theoretical predictions obtained in the framework of collinear approach [1,[37][38][39] (see also section 2.5 in the recent review [46]). However, we would like to note that there is a quite similar situation between the exclusive reduced cross section σ red (x, Q 2 ) and F 2 (x, Q 2 ): experimental data for both of these observables are in good agreement with the corresponding theoretical predictions (see [75] and discussions therein), where the calculated SF F L (x, Q 2 ) is known to be very sensitive to low x resummation (see, for example, Section 9.3 in the recent review [46]). But, apart of F c L (x, Q 2 ), the SF F L (x, Q 2 ) is measured at the HERA (see [76] and references therein). Therefore, it seems that in order to understand the observed difference between the predictions for R c (x, Q 2 ) andR c (x, Q 2 ) ratios at large Q 2 one have to investigate first the SF F L (x, Q 2 ) using the same approaches, that is out of our present consideration. We plan to perform such investigation in forthcoming study and then, having experience in the latter, return to the study its heavy quark parts and, correspondingly, the ratio R c (x, Q 2 ).

Conclusions
We have studied the heavy quark production processes with using the transverse momentum dependent gluon distribution function in a proton obtained recently [29] using the Kimber-Martin-Ryskin prescription from the Bessel-inspired behavior of parton densities at small Bjorken x values. The Bessel-like behavior of parton densities at small Bjorken x was obtained [30][31][32][33] in-turn in the case of the flat initial conditions for DGLAP evolution equations in the double scaling QCD approximation. To construct the TMD parton distributions we implemented [29] the different treatments of kinematical constraint, reflecting the angular and strong ordering conditions and discussed the relations between the differential and integral formulation of the KMR approach. Additionally, we have tested the TMD gluon density in a proton obtained from the numerical solution of the CCFM evolution equation, which smoothly interpolates between the small-x BFKL dynamics and large-x DGLAP ones.
We have considered the (reduced) cross sections σ QQ red (where Q = c, b) and charm and beauty contributions to the deep inelastic proton SFs F 2 (x, Q 2 ) and F L (x, Q 2 ). To show an importanse of the longitudinal structure function F c L (x, Q 2 ) and F b L (x, Q 2 ), we compare the results for σ cc red and σ bb red with the SFs and F c 2 (x, Q 2 ) and F b 2 (x, Q 2 ). We achieved a good agreement between the HERA experimental data for these observables and our theoretical predictions and demostrated the importance of the contributions of F c L (x, Q 2 ) and F b L (x, Q 2 ) at small x. Concerning the ratio of the proton SFs, namely, R c (x, Q 2 ) = F c L (x, Q 2 )/F c 2 (x, Q 2 ), we show that the results of k T -factorization calculations are similar to the ones obtained beyond LO of collinear perturbation theory. This effect is clearly visible for Q 2 ≤ 12 GeV 2 . However, starting with Q 2 ≥ 12 GeV 2 , the k T -factorization leads to larger values for the ratioR c (x, Q 2 ), that needs an additional investigations.
As we discussed already in Section 4.2, as the next step, we plan to study the longitudinal structure function F L (x, Q 2 ) and compare our future results with the previous ones [61,77] and [59,78,79], obtained in the framework of k T -factorization and collinear perturbation theory, respectively. This study is important in itself and will provide some kind of clue for solving the problem of differences in the predictions for the ratio R c (x, Q 2 ) obtained in the framework of k T -factorization approach and conventional (collinear) QCD factorization (see Section 4.2). Moreover, we plan to extend the present analysis beyond the LO approximation. We will obtain the results for the NLO TMD parton densities using the corresponding NLO results [30][31][32] for the standard PDFs in the generalized DAS approach. We will accept also the results for the NLO matrix elements (see [27,80] and references and discussions therein). Such results are seems to be extremely important for future experiments, in particular, for experiments at the Electron-Ion Collider (EIC) and Electron-Ion Collider in China (EiCC) (see [83,84] and discussions and references therein). So, at the EIC, an essential low x (up to x ∼ 10 −4 ) region is expected to be probed, thus providing us with new and precise data for DIS SFs, especially data for longitudinal SF F L (x, Q 2 ). The EiCC could provide a new information on the light and sea quark density in a proton, that is, of course, important to produce and update the theoretical high-order predictions for F L (x, Q 2 ). Moreover, EIC and EiCC measurements could be important to distinguish between the different non-collinear QCD evolution scenarios widely discussed at present (see, for example, review [2]).
where the SFs, obtained in the generalized DAS approach, were marked asF Q k (x, Q 2 ). Here ⊗ is the Mellin convolution: where It can be represented aŝ where f g (x, Q 2 ) is given in (31). In fact, the ratioR Q (x, Q 2 ) depends slowly on nonperturbative input f g (x, Q 2 ), which contribute to the both numerator and denominator of the ratioR Q (x, Q 2 ).
Using the results of k T -factorization and BFKL approach [5,85] (see also [86,87]), the results for the high energy limit of collinear coefficient functions of the heavy quark production process in all orders of perturbation theory were obtained. Thus, using the results [85], below we give formulas for the high energy asymptotics of collinear coefficient functions of the heavy quark production process in the first three orders of the perturbation theory.

LO
Taking the LO Wilson coefficient (20), results (24) and (25) for on-shell coefficient functions and the PDFs considered in the section 3.2, we have for the ratioR Q LO (x, Q 2 ): where the dependence of gluon density f g (x, Q 2 ) should be rather week. In fact, there is no x-dependence at all (see Fig. 5), which is associated with the property of the Mellin convolution (A2) in the low x region (see Appendix B).

NLO
Through NLO, we have The NLO coefficient functions B k,g (x, a) of photon-gluon fusion subprocess are rather lengthy and not published in print; they are only available as computer codes [88]. Following [81,82], it is sufficient to work in the high energy regime, defined by x 1, where they assume the compact form 4 [86,86,87]: being the dilogarithmic function. We would like to note that B k,g (x, a) (see (24) and (25)): So, at the NLO we have for the ratioR Q (x, Q 2 ): where the LO gluon density f g (x, Q 2 ) given by (31) is used because its contribution to the ratioR Q NLO (x, Q 2 ) is strongly suppressed (see Appendix B).

NNLO
Following to the results [87], we have for the coefficient function: where in the high energy regime the coeffcicient B k,g (x, a) has the compact form: where J(a) and I(a) are defined by (A10) and (A11), respectively, and where t is given in (A10) and are the trilogarithmic function Li 3 (x) and Nilsen Polylogarithm S 1,2 (x) (see [89]). The results for K(a) in the form of harmonic Polylogarithms [90] can be found in [87]. So, at the NNLO we have for the ratioR Q (x, Q 2 ): where the LO gluon density f g (x, Q 2 ) given by (31) is used because its contribution to the ratioR Q NNLO (x, Q 2 ) is strongly suppressed (see Appendix B).

Appendix B. Collinear results in DAS approach
The use of the DAS approach for the PDFs makes it possible to significantly simplify the formulas for the relationR Q (x, Q 2 ). We will show this below.
Taking the results (24) and (25) for on-shell coefficient functions and the PDFs considered in the section 3.2, it is easily to obtain the following LO results in the generalized DAS approach (see [81] where M k,g (1, Q 2 , m 2 ) is the first Mellin moment (n = 1) (see (A10)), where the Mellin moments can be defined as where x 2 are given by (A3). In fact, the non-perturbative input f g (x, Q 2 ) does cancel in theR f ratio and we havê in the case, where the moments M k,g (1, Q 2 , m 2 Q ) have no singularities at n → 1.

LO
Taking the integral (B2), which becomes to be (A13) at the LO, we can obtain the results (A9), using (see [81]) the following auxiliary formulas 6 if m = 0 So, at LO the small-x approximation formula (A5) thus readŝ which is x-independent, that is in full agreement with the numerical evaluation of the R Q LO (x, Q 2 ) in (A5).

NLO
At NLO, the coefficient function C g k (x) has the form (A6) with the NLO coefficients B (1) k,g (x, a) given by (A7) and (A8). Its moments M k,g (n, Q 2 , µ 2 ) exhibit the corresponding structure The Mellin transforms of B (1) k,g (x, a) exhibit singularities in the limit n → 1, which lead to modifications in (B1). As was shown [91], the terms involving 1/δ at n = 1+δ → 1 (which correspond to singularities of the Mellin moments M k,g (n) (see (B2)) at n → 1) depend on the exact form of the asymptotic low-x behavior encoded in f g (x, µ 2 ). Using the results for f g (x, µ 2 ) from (31), we obtain the modification is simple modification (see [81] and discussions therein) where ρ g (x, µ) are given by (38).
Using the identity where the ratio has some x-dependence coming from the corresponding x-dependence of δ + in (B12). The x-dependence is in rather good agreement with the numerical results in (A20).

NNLO
At NNLO, the coefficient function C g k (x) has the form (A15) with the NNLO coefficients B (2) k,g (x, a) given in Eqs. (A16) and (A17). Its moment M 2,g (n, Q 2 , µ 2 ) exhibits the corresponding structure k,g (x, a) exhibit singularities in the limit n → 1, which has the form 1 (n − 1) 2 → 1 δ 2 ++ = 1 ρ 2 (x, µ) where all definitions can be found in (32) and (33). So, by analogy with the NLO case, we havẽ M k,g (1, Q 2 , µ 2 ) = e 2 Q a s (µ 2 ) B (0) k,g (1, a) + a s (µ 2 )B (1) k,g (1, a) + a 2 s (µ 2 )B (2) k,g (1, a) . (B18) Using the identity we find the Mellin transform (B2) of (A7) to bẽ where the ratio has some x-dependence coming from the corresponding x-dependence of δ ++ in (B19). The x-dependence is in rather good agreement with the numerical results in (A20). Figure 1: The reduced charm cross sections σ cc red (x, Q 2 ) as a function of x calculated at different Q 2 values. The predictions obtained with analytical TMD gluon density in a proton and CCFM-evolved one are shown by the solid green and yellow curves, respectively. The shaded bands correspond to the scale uncertainties of our calculations. The dashed curves represent the contributions from SF F c 2 (x, Q 2 ), as it is described in the text. The experimental data are from H1 and ZEUS [1]. Figure 2: The reduced beauty cross sections σ bb red (x, Q 2 ) as a function of x calculated at different Q 2 values. Notation of all curves is the same as in Fig. 1. The experimental data are from H1 and ZEUS [1]. Figure 3: The charm contribution to the proton structure function F 2 (x, Q 2 ) as a function of x calculated at different Q 2 values. Notation of all curves is the same as in Fig. 1. The experimental data are from ZEUS [37] and H1 [38]. Figure 4: The beauty contribution to the proton structure function F 2 (x, Q 2 ) as a function of x calculated at different Q 2 values. Notation of all curves is the same as in Fig. 1. The experimental data are from ZEUS [37] and H1 [39].