Transverse momentum dependent parton densities in a proton from the generalized DAS approach

We use the Bessel-inspired behavior of parton densities at small Bjorken $x$ values, obtained in the case of the flat initial conditions for DGLAP evolution equations in the double scaling QCD approximation (DAS), to evaluate the transverse momentum dependent (TMD, or unintegrated) quark and gluon distribution functions in a proton. The calculations are performed analytically using the Kimber-Martin-Ryskin (KMR) prescription with different implementation of kinematical constraint, reflecting the angular and strong ordering conditions. The relations between the differential and integral formulation of the KMR approach is discussed. Several phenomenological applications of the proposed TMD parton densities to the LHC processes are given.


Introduction
A theoretical description of a number of high energy processes at hadron colliders which proceed with large momentum transfer and containing multiple hard scales can be obtained with transverse momentum dependent (TMD), or unintegrated, parton (quark and/or gluon) density functions in a proton [1]. These quantities depend on the fraction x of the proton longitudinal momentum carried by a parton, the two-dimensional parton transverse momentum k 2 T and hard scale µ 2 of the hard process and encode nonperturbative information on proton structure, including transverse momentum and polarization degrees of freedom. They are related to the physical cross sections and other observables, measured in the collider experiments, via TMD factorization theorems in Quantum Chromodynamics (QCD). The latter provide the necessary framework to separate hard partonic physics, described with a perturbative QCD expansion, from soft hadronic physics. At present, there are number of factorization approaches which incorporate the transverse momentum dependence in the parton distributions: for example, the Collins-Soper-Sterman approach [2] (or TMD factorization), designed for semi-inclusive processes with a finite and non-zero ratio between the hard scale µ 2 and total energy s, and high-energy factorization [3] (or k T -factorization [4]) approach, valid in the limit of a fixed hard scale and high energy.
In the high-energy factorization, the TMD gluon density satisfies the Balitsky-Fadin-Kuraev-Lipatov (BFKL) [5] or Ciafaloni-Catani-Fiorani-Marchesini (CCFM) [6] evolution equations, which resum large terms proportional to α n s ln n s ∼ α n s ln n 1/x, important at high energies s (or, equivalently, at small x). Thus, one can effectively take into account higher-order radiative corrections to the production cross sections (namely, a part of NLO + NNLO +... terms corresponding to the initial-state real gluon emissions). The CCFM equation takes into account additional terms proportional to α n s ln n 1/(1 − x) and therefore is valid at both small and large x. The number of phenomenological applications of the high-energy factorization and CCFM evolution is known in the literature (see, for example, [7][8][9][10][11][12][13][14][15][16][17] and references therein).
There are also other approaches determining the TMD gluon and quark density functions in a proton. So, one can evaluate them using the schemes based on the conventional Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) [18] equations, namely the Parton Branching (PB) approach [19,20] and Kimber-Martin-Ryskin (KMR) prescription [21]. Former gives numerical iterative solution of the DGLAP evolution equations for collinear and TMD parton density functions upon using a concept of resolvable and non-resolvable branchings and by applying Sudakov formalism to describe the parton evolution from one scale to another without resolvable branching. The splitting kinematics at each branching vertex is described by the DGLAP equations and angular ordering condition for parton emissions can be applied instead of usual DGLAP ordering in virtuality. The latter is a formalism invented for constructing the TMD parton distributions from well-known conventional (collinear) parton density functions (PDFs) under the key assumption that the transverse momentum dependence of the parton distributions enters only at the last evolution step. The KMR procedure is believed to take into account effectively the major part of next-to-leading logarithmic (NLL) terms α s (α s ln µ 2 ) n−1 compared to the leading logarithmic approximation (LLA), where terms proportional to α n s ln n µ 2 are taken into account. The KMR approach is currently explored [22] at next-to-leading order (NLO) and commonly used in the phenomenological applications (see, for example, [8,[10][11][12][13][14][15][16] and references therein), where the standard proton PDFs (as obtained, for example, by the NNPDF [23] or CTEQ [24] Collaborations) were taken as an input numerically. The relation between the PB and KMR scenarios was discussed [25] and the connection between the PB and CCFM approaches was established very recently [26].
The KMR formalism is used in the present paper for analytical calculations of the TMD quark and gluon distributions in a proton. The calculations are based on the expressions [27][28][29] for conventional PDFs obtained in the generalized double asymptotic scaling (DAS) approximation [28][29][30]. The latter is connected to the asymptotic behaviour of the DGLAP evolution discovered many years ago [31]. As it was shown, flat initial conditions for DGLAP equations, applied in the generalized DAS scheme, lead to the Bessel-like behaviour for the proton PDFs at small x. Using the results [27][28][29], we derive the analytical expressions for the TMD quark and gluon densities at leading order (LO) and present them in a quite compact form. We implement different treatment of the kinematical constraint involved in the KMR prescription (namely, angular and strong ordering conditions) and discuss the relations between the differential and integral formulation of the KMR procedure pointed out recently in [32]. Finally, we present some phenomenological applications of obtained TMD parton densities to hard LHC processes, sensitive to the quark and gluon content of the proton. To be precise, using the k Tfactorization QCD approach, we consider the inclusive production of b-jets and bb-dijets at √ s = 7 TeV, inclusive production of the Higgs bosons (in the diphoton decay mode) at √ s = 13 TeV and charm and beauty contributions to the deep inelastic proton structure function The outline of our paper is following. In Section 2 we briefly describe our theoretical input. The calculations are explained in detail in Section 3, where we also present the modifications of PDFs and TMDs at low Q 2 -values and outside of the standard small x range. Section 4 present numerical results and discussions. Section 5 contains our conclusions.

Theoretical framework
Since the KMR approach is based on the standard PDFs, here we present a review of small x behaviour of parton densities. Moreover, we introduce the basic formulas of the KMR approach itself.

PDFs and proton structure function
The fairly reasonable agreement between HERA data [33][34][35][36][37] and the results of NLO perturbative QCD evaluations is observed for Q 2 ≥ 2 GeV 2 (see reviews [38,39] and references therein). Therefore, it can be concluded that pQCD is capable of describing the evolution of proton structure function F 2 (x, Q 2 ) and its derivatives down to very low Q 2 values.
It was pointed out [40] that the HERA small-x data can be well interpreted in terms of the so-called doubled asymptotic scaling (DAS) phenomenon related to the asymptotic behaviour of the DGLAP evolution discovered many years ago [31]. The study [40] was extended [28][29][30] to include the finite parts of anomalous dimensions (ADs) of Wilson operators and Wilson coefficients 1 . This led to predictions [27][28][29] of the small-x asymptotic form of PDFs in the framework of the DGLAP dynamics, which were obtained starting at some Q 2 0 with the flat function where f a are PDFs multiplied by x, a = q or g and A a are unknown parameters to be determined from the data. We refer to the approach of [28][29][30] as generalized DAS approximation. In this approach the flat initial conditions (1) determine the basic role of the AD singular parts as in the standard DAS case, whereas the contributions coming from AD finite parts and Wilson coefficients can be considered as corrections which are, however, important for achieving better agreement with experimental data.
Hereafter we consider for simplicity only the LO approximation. The structure function F 2 (x, Q 2 ) and PDFs have the following form (see [28,29]) where e = f i=1 e 2 i /f is an average charge squared with f being a number of active (massless) quark flavors. The small-x asymptotic expressions for sea quark and gluon densities f a (x, µ 2 ) can be written as follows (both the LO and NLO results and their derivation can be found [27,28]): where for the color SU (N c ) group, I ν (σ) andĨ ν (σ) (ν = 0, 1) are the combinations of the modified Bessel functions (at s ≥ 0, i.e. µ 2 ≥ Q 2 0 ) and the 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 and β 0 = 11 − (2/3)f 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 [27,41]; they were obtained 2 for α s (M Z ) = 0.1168.
It is convenient to show the following expressions: 2 In the future, by using (2) and (3) and results of [42] we plan to perform the combined fits to the H1 and ZEUS experimental data [37] and [43] for the DIS structure function F 2 (x, Q 2 ) and its charm part F c 2 (x, Q 2 ), respectively.

Kimber-Martin-Ryskin approach
The expressions (3) -(6) can be used as an input for the KMR procedure [21], giving us the possibility to calculate the TMD parton density functions in a proton. Let k ⊥ ≡ k, then the TMD parton distributions in the differential f (d) a (x, k 2 , µ 2 ) and integral f (i) a (x, k 2 , µ 2 ) formulation of the KMR approach can be written as where D a (x, µ 2 ) are the conventional PDFs, f a (x, µ 2 ) = xD a (x, µ 2 ), which obey the DGLAP equations (see (2.1) in [32]): with the splitting functions The Sudakov form factor T a (µ 2 , k 2 ) has the following form (see (2.4) in [32]):

TMDs from differential formulation of KMR approach
Now we can use (8) to find the results for TMD parton densities without derivatives. Derivation of T a (µ 2 , k 2 ) is as follows and derivations of conventional PDFs are as follows So, the result for the TMD parton densities reads The expressions for conventional PDFs and Sudakov form factors are given in the Sections 2 and 3.1, respectively. Since the value ofd + is negative and the factord + /ρ a is large at low x, the TMDs f With increasing x, the results for the latter can be negative.

TMDs from integral formulation of KMR approach
Now we can use (9) to find the results for TMDs without derivatives. After some algebra we have where Using the relations from (7), we can obtain: where

Cut-off parameter ∆
For the phenomenological applications, the cut-off parameter ∆ usually has one of two basic forms: that reflects the two cases: ∆ 1 is in the strong ordering, ∆ 2 is in the angular ordering (see [32]). In all above cases, except the results for T a (µ 2 , k 2 ), we can simply replace the parameter ∆ by ∆ 1 and/or ∆ 2 . For the Sudakov form factors, we note that the parameters ∆ i (with i = 1, 2) contribute to the integrand in (12) and, thus, their momentum dependence changes the results in (14). To perform the correct evaluation of the integral (12), we should recalculate the p 2 integration in (12). So, we have The analytic evaluation of T (i) 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 (i) a (µ 2 , k 2 ).

Comparison of f
To perform the comparison between the two obtained expressions, it is convenient to rewrite (19) as where Using the relations from (7), we can obtain: It is convenient to compare the results for f or So, the difference is regular at ∆ → 0 and it does not give large contributions. This is clearly illustrated in Fig. 1, where we show the ratio of gluon densities f a (x, k 2 , µ 2 ) calculated with angular ordering condition applied as a function of k 2 T ≡ k 2 at several values of x and µ 2 . As one can see, with increasing k 2 the difference between these two approaches becomes more pronounced.

Infrared modification of the strong coupling
The equations (3) at s < 0 were used in [29], where the higher-twist corrections through twist six were added to find good agreement with the experimental data for the deep inelatic proton structure function F 2 (x, Q 2 ) for Q 2 ≥ 0.5 GeV 2 . However, such application is not so useful here because the case with s < 0 may lead to the negative TMDs and, hence, to the negative cross sections of the physical processes.
To overcome these problems, which emerge at small k 2 values, we investigate an alternative possibility following to [27]; namely, the modification of the strong-coupling constant in the infrared region. Specifically, we consider two modifications, which effectively increase the argument of the strong coupling constant at small µ 2 values, in accordance with [45,46]. In the first case, which is more phenomenological, we introduce a freezing of the strong-coupling constant by changing its argument as µ 2 → µ 2 + M 2 ρ , where M ρ is the ρ meson mass [47]. Thus, in the formulae of Section 3 we introduce the following replacement The second possibility is based on the idea by Shirkov and Solovtsov [48] (see also the recent reviews [49] and the references therein) regarding the analyticity of the strong coupling that leads to an additional power dependence. In this case, the QCD coupling α s (µ 2 ) appearing in the formulae of the previous Sections is to be replaced as Such replacements have been done [27], where we took the normalizations magnitudes A g and A q . As we can see from [27,41], the fits based on the frozen and analytic strongcoupling constants are very similar and describe the F 2 (x, Q 2 ) data in the small-Q 2 range significantly better than the canonical fit.

Beyond small x
In the phenomenological applications (see Section 4) the calculated TMD parton densities will be used to predict the cross sections of several high-energy processes. According to k T -factorization approach [3,4], the theoretical predictions for the cross sections can be obtained by convolution of these TMD parton densities and the 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).
Our TMD parton densities are exactly expressed through the conventional PDFs f a (x, µ 2 ) as it was shown in the equations (26) and (20). Then, the densities f a (x, µ 2 ) listed in Section 2.1 should be extended in the following form [50,51] (see, for example, the recent paper [52], where similar extension has been done in the case of EMC effect from the study of shadowing [53] at low x to antishadowing effect at Note that such form was successfully used in the conventional PDF parametrizations (see [51,54]). The value of β a (0) can be estimated from the quark counting rules [55]: 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, [56][57][58]) and the results for these values may be quite different, because various groups producing PDF sets use different sets of experimental data or take some privilege for their parts. Moreover, the difference can be attributed to various choices of the initial condition µ 2 0 of the µ 2evolution but it should be not so strong because the µ 2 -dependence is double-logarithmic.
It is convenient to assume that similar relations take place just beyond the standard low x range (x ≤ 0.05). Thus, the TMD parton densities can be modified in the form, similar to (33), that leads to In our analysis, the numerical values of β g (0) have been extracted from the fit to the inclusive b-jet production data taken by the CMS [59] and ATLAS [60] Collaborations in pp collisions at √ s = 7 TeV (see Section 4.1 below). We find that best description of the leading b-jet transverse momentum distributions in a whole kinematical region is achieved with β g (0) = 3.03 and β g (0) = 5.77 for "frozen" and analytic strong coupling constant (37) and (38), respectively. We see that the obtained results are close to ones in (34). The TMD gluon densities in a proton obtained with appropriate treatment of the strong coupling and β g (0) are shown in Fig. 1 as a function of transverse momentum k 2 T for different values of proton longitudinal momentum fraction x and hard scale µ 2 . We have used the integral formulation of KMR procedure as given by (36) for illustration. The solid green and yellow curves correspond to the results obtained with "frozen" and analytic coupling constant with angular ordering condition, while corresponding dashed curves represent the results obtained with strong ordering condition. As one can see, the strong ordering condition leads to a steep drop of the gluon densities beyond the scale µ 2 . It contrasts with angular ordering, where the gluon transverse momentum is allowed to be larger than µ 2 (see [32]). We also show here the TMD gluon distributions calculated numerically in the traditional KMR scenario, where the conventional parton densities from standard MMHT'2014 (LO) set [61] were used as an input (red curves on Fig. 1). The results of our analytical calculations are nicely agree with the latter for k 2 T ≥ 10 GeV in wide x region (up to x ≤ 0.05), that demonstrates the applicability of the generalized DAS approximation. At k 2 T < µ 2 0 ∼ 1 GeV 2 the numerically calculated KMR gluon density is modelled to be a flat according to the prescription [21] under strong normalization condition which is often used in the KMR scheme. Such determination, of course, leads to a low k 2 T plateau, clearly seen in Fig. 1. In contrast, our formalism with appropriate modifications of strong coupling as described above results in continuous TMD quark and gluon density functions, well defined in a whole k 2 T region. Below we will consider the phenomenological consequences of our approach.

Phenomenological applications
We are now in a position to apply the obtained TMD parton densities in a proton to several hard QCD processes studied at hadron colliders. In the present paper we consider the inclusive production of b-jets and Higgs bosons at the LHC conditions and charm and beauty contributions to the deep inelastic proton structure function F 2 (x, Q 2 ) measured in ep collisions at HERA. These processes have been already investigated within the k T -factorization approach and found to be strongly sensitive to the gluon content of the proton. To calculate the total and differential cross sections of b-jets, Higgs boson production and proton structure functions F c 2 (x, Q 2 ) and F b 2 (x, Q 2 ) we strictly follow our previous considerations [11,12,15,16,62]. Everywhere below we have used one-loop formula for the strong coupling constant with n f = 4 active quark flavors and Λ QCD = 143 MeV (that corresponds to α s (m 2 Z ) = 0.1168) for analytically calculated TMD quark and gluon densities as described above and apply n f = 5 with α s (m 2 Z ) = 0.13 for KMR partons evaluated numerically. The latter choice is dictated by the parameter setup employed in the MMHT'2014 (LO) PDFs [61], used here as an input for KMR procedure.

Inclusive b-jet production at the LHC
Following [15,16], our consideration is based on the leading off-shell (depending on the transverse momenta of incoming particles) gluon fusion subprocess g * (k 1 ) + g * (k 2 ) → b(p 1 ) +b(p 2 ), where the four-momenta of all particles are indicated in parentheses. According to the k T -factorization prescription [3,4], the corresponding cross section can be written as where σ * (x 1 , x 2 , k 2 1T , k 2 2T , µ 2 ) is the off-shell partonic cross section and k 2 1T and k 2 2T are the non-zero two-dimensional transverse momenta of incoming partons. The detailed description of the calculation steps (including the evaluation of the off-shell amplitudes) can be found in [15,16]. Here we only specify the essential numerical parameters. So, following [63], we set the b-quark mass m b = 4.78 GeV and, as it often done in pQCD calculations, choose the default renormalization and factorization scales µ R and µ F to be equal to leading b-jet transverse momentum. The calculations were performed using newly developed Monte-Carlo event generator pegasus [64].
The CMS Collaboration has measured the double differential cross section dσ/dp T dy of inclusive b-jet production at √ s = 7 TeV in five b-jet rapidity regions, namely, |y| < 0.5, 0.5 < |y| < 1, 1 < |y| < 1.5, 1.5 < |y| < 2 and 2 < |y| < 2.2 as a function of the leading b-jet transverse momentum [59]. In the ATLAS analysis [60], the inclusive bjet cross section has been measured as a function of transverse momentum p T in the range 20 < p T < 400 GeV and rapidity in the range |y| < 2.1. In addition, the bb-dijet cross section has been measured as a function of the dijet invariant mass M in the range 110 < M < 760 GeV, azimuthal angle difference ∆φ between the two b-jets and angular variable χ = exp |y 1 − y 2 | for jets with p T > 40 GeV in two dijet mass regions.
The results of our calculations are shown in Figs. 3 -5 in comparison with the CMS and ATLAS data [59,60]. The solid green and yellow histograms were obtained with the TMD gluon density as given by (20) -(23) with "frozen" and analytic QCD coupling by fixing both the renormalization and factorization scales at their default values. The red histograms represent the results obtained with the numerically calculated KMR gluon distributions. The shaded bands correspond to scale uncertainties of these predictions. As usual, the latter have been estimated by varying the scales µ R and µ F by a factor of 2 around their default values. We have obtained a good description of the b-jet transverse momentum distributions in each of the rapidity subdivisions, both in normalization and the shape. Our predictions are only tend to slightly underestimate the measured cross sections at very high transverse momenta p T ∼ 200 − 400 GeV, but they agree with the data within the theoretical and experimental uncertainties. The results obtained with numerically calculated KMR gluon density agree with data and analytical TMD gluons at low and moderate transverse momenta, but overestimate both the CMS and ATLAS data at p T > 100 GeV, especially at forward rapidities. We note that here the essentially large x region is probed, so the better description of the data achieved with analytical TMD gluon distributions demonstrates that their large-x extension, as described above in Section 3.7, is rather reasonable.
All the considered TMD gluons show good agreement with the bb-dijet cross sections measured by the ATLAS Collaboration. In particular, the good description of the ∆φ distribution is remarkable, since the latter is known to be a strongly sensitive to the k 2 T shape of the TMD gluon density (see [15,16] and references therein). As it was expected, the χ distribution flattens for large invariant masses M . Note that here an additional acceptance requirement, that restricts the boost of the dijet system to |y boost | = |y 1 + y 2 |/2 < 1.1, has been applied for χ measurements. This requirement significantly reduces [60] the sensitivity to gluon density function at small x and all theoretical predictions for χ distributions are practically coincide. Thus, we conclude that our analytical TMD parton densities given by (20) -(23) does not contradict available LHC data on b-jet production.

Inclusive Higgs boson production at the LHC
Our consideration is mainly based on the off-shell amplitude of the gluon-gluon fusion subprocess g * (k 1 ) + g * (k 2 ) → H(p) calculated using the effective Lagrangian [65,66] for the Higgs coupling to gluons and extended recently to the subsequent H → γγ, H → ZZ * → 4l (where l = e or µ) and H → W + W − → e ± µ ∓ νν decays. The details of the calculations are explained in [11,12] and here we strictly follow our previous consideration. Everywhere below, we set the Higgs boson mass m H = 125.1 GeV and its full decay width Γ H = 4.3 MeV. The default values of the renormalization and factorization scales are chosen to be equal to Higgs mass. The cross sections were produced with Monte-Carlo generator pegasus [64].
The latest measurements of the inclusive Higgs boson production (in the diphoton decay mode) were performed by the CMS [67] and ATLAS [68] Collaborations at the LHC energy √ s = 13 TeV. In the CMS analysis, two isolated final state photons originating from the Higgs boson decays are required to have pseudorapidities |η γ | < 2.5, excluding the region 1.4442 < |η γ | < 1.566. Additionally, photons with largest and next-to-largest transverse momentum p γ T (so-called leading and subleading photons) must satisfy the conditions of p γ T /M γγ > 1/3 and p γ T /M γγ > 1/4 respectively, where M γγ is the diphoton pair mass, M γγ > 90 GeV. In the ATLAS measurement [68] both of these decay photons must have pseudorapidities |η γ | < 2.37 (excluding 1.37 < |η γ | < 1.52) with the leading (subleading) photon satisfying p γ T /M γγ > 0.35 and p γ T /M γγ > 0.25, while invariant mass M γγ is required to be 105 < M γγ < 160 GeV. We have implemented experimental setup in our numerical program. The Higgs transverse momentum p T , absolute value of the rapidity y and cosine of photon helicity angle cos θ * (in the Collins-Soper frame) were measured [67,68]. Both p T and y probe the production mechanism and parton distribution functions in a proton, while cos θ * is related to spin-CP nature of the decaying Higgs boson.
The results of our calculations are shown in Figs. 6 and 7 in comparison with latest LHC data. One can see that our predictions with both analytic and "frozen" treatment of QCD coupling reasonably agree with the data for all considered kinematical observables, although some tendency to slightly overestimate the LHC data in the low p T region is observed. This tendency results in a some ovestimation of the rapidity and photon helicity angle distributions, but the predictions are still agree with the data within the experimental and theoretical uncertainties, calculated as it was described above. The scale dependence of our predictions, of course, exceeds the uncertainties of conventional higher-order pQCD calculations (which are about of 10 -11%). However, it could be easily understood because only the tree-level LO hard scaterring amplitudes are involved. The strong drop in the | cos θ * | distribution around | cos θ * | ∼ 0.6 is due to the fiducial requirement on the photon system originating from the scalar Higgs boson decay. The calculations based on the numerically evaluated KMR gluon density agree well with the data at low p T and tend to overshoot them at high transverse momenta. Note that we added to our results contributions from weak boson fusion (W + W − → H and ZZ → H), associated HZ or HW ± production and associated ttH production (grey shaded bands in Figs. 6 and 7). These contributions are essential at high p T and have been calculated in the conventional pQCD approach with the NLO accuracy. We take them from [67,68].
Once again, we can conclude that the analytical expressions for TMD parton densities (20) -(23) does not contradict available LHC data in the probed kinematical region, where µ 2 ∼ m 2 H .

Proton structure functions
The important information on the quark and gluon structure of proton can be also extracted from the data on deep inelastic ep scattering. Its differential cross section can be presented in the simple form: where F 2 (x, Q 2 ) and F L (x, Q 2 ) are the proton transverse and longitudinal structure functions, x and y are the usual Bjorken scaling variables. The charm and beauty contributions to F 2 (x, Q 2 ) are described through perturbative production of charm or beauty quarks and, therefore, directly related with the gluon content of the proton. Our evaluation below is based on the formulas [62] and here we again strictly follow our previous consideration in all aspects. We only note that the charm and beauty masses are set to be equal to m c = 1.65 GeV and m b = 4.78 GeV [63].
Our results are shown in Figs. 8 and 9 in comparison with the latest ZEUS [69] and H1 [70,71] data. We find that the predictions obtained with "frozen" strong coupling are in perfect 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 , both in normalization and shape. These predictions slightly overshoot the ones obtained with analytic treatment of the QCD coupling constant. The difference between these two approaches becomes more clearly pronounced at small x and low Q 2 values, Q 2 ≤ 10 GeV 2 . The predictions based on the analytic QCD coupling still agree with the data on F b 2 (x, Q 2 ) within the theoretical and experimental uncertainties, but clearly underestimate the data on F c 2 (x, Q 2 ), especially at low Q 2 . However, we note that some reasonable variation in charmed quark mass m c = 1.65 ± 0.2 GeV can almost eliminate the visible disagreement (not shown in Figs. 8 and 9). The traditional KMR approach underestimates the HERA data on both F c 2 (x, Q 2 ) and F b 2 (x, Q 2 ) at small x and relatively low Q 2 ≤ 30 − 60 GeV 2 , where the higher-order QCD corrections are known to be important. The difference between all the theoretical predictions becomes negligible with increasing of Q 2 . Thus, we can conclude that the proposed analytical calculations of the TMD parton densities in a proton does not contradict the latest HERA data, although the best description of the latter (with the default parameter set) is achieved with "frozen" QCD coupling constant.

Conclusions
We presented the analytical calculations of the transverse momentum dependent parton densities in a proton. These calculations are based on the Bessel-inspired behavior of parton densities at small Bjorken x, obtained 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 applied the leading-order Kimber-Martin-Ryskin approach, which is widely used in the phenomenological applications. We implemented 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. Finally, we demonstrated that the calculated TMD parton distributions does not contradict the LHC data on inclusive b-jet production at √ s = 7 TeV, inclusive Higg boson production (in diphoton decay mode) at √ s = 13 TeV and latest HERA data on the charm and beauty contributions to the deep inelastic proton structure function F 2 (x, Q 2 ) in a wide region of x and Q 2 .
As the next step, we plan to study the longitudinal structure function F L (x, Q 2 ) and also to provide predictions for its heavy quark parts, F c L (x, Q 2 ) and F b L (x, Q 2 ), which will be compare with our other predictions [62]. 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 [27][28][29] for the standard PDFs in the generalized DAS approach. We will accept also the results for the NLO matrix elements (see [25,72] and references and discussions therein).
x = 10 -5 x = 10 -4 x = 10 -3 x = 10 -2 x = 10 -1 Figure 1: The ratio of the TMD gluon densities in a proton obtained using integral and differential formulation of the KMR approach as a function of the gluon transverse momentum k 2 T at different values of longitudinal momentum fraction x and hard scale µ 2 . The angular ordering condition was applied.         The charm contribution to the proton structure function F 2 (x, Q 2 ) as a function of x calculated at different Q 2 . Notation of curves is the same as in Fig. 1. The experimental data are from ZEUS [69] and H1 [70]. Figure 9: The beauty contribution to the proton structure function F 2 (x, Q 2 ) as a function of x calculated at different Q 2 . Notation of curves is the same as in Fig. 1. The experimental data are from ZEUS [69] and H1 [71].