Transverse momentum dependent gluon density from DIS precision data

The combined measurements of proton's structure functions in deeply inelastic scattering at the HERA collider provide high-precision data capable of constraining parton density functions over a wide range of the kinematic variables. We perform fits to these data using transverse momentum dependent QCD factorization and CCFM evolution. The results of the fits to precision measurements are used to make a determination of the nonperturbative transverse momentum dependent gluon density function, including experimental and theoretical uncertainties. We present an application of this density function to vector boson + jet production processes at the LHC.


I. INTRODUCTION
The high-precision combined HERA data [1] for proton's deeply inelastic scattering (DIS) structure functions constrain parton density functions (pdfs) over a wide range of the kinematic variables. These data have been used for determinations of the collinear pdfs and related studies at the LHC [2].
On the other hand, QCD applications to multiple-scale scattering problems and complex final-state observables require in general formulations of factorization [3] which involve transverse-momentum dependent (TMD), or unintegrated, parton density and parton decay functions [4][5][6][7]. TMD pdfs are necessary to describe appropriately nonperturbative physics and to control perturbative large logarithms to higher orders of perturbation theory.
The purpose of this work is to use the combined DIS data on structure functions [1] and charm production [8] for determination of TMD pdfs. A general program for TMD pdfs phenomenology has been proposed in [4]. Our work has a more limited scope than this program as we limit ourselves to considering DIS data in the small-x kinematic region. On the other hand, from the point of view of TMD pdfs, this region is interesting because a welldefined form of TMD factorization holds at high energy [9], which has been applied to sum small-x logarithmic corrections to DIS to all orders in α s at leading and next-to-leading ln x level [10][11][12]. Furthermore, given the high precision of the combined data [1,8], this analysis provides a compelling test of the TMD approach and of the limitations of the logarithmic approximations used at small x. This is to be contrasted with earlier analyses [13,14] based on older and much less precise structure function measurements.
The high-energy factorization [9] expresses the heavy-quark leptoproduction cross section in terms of the TMD gluon density via well-prescribed, calculable perturbative coefficients. This framework is extended to DIS structure functions in [11]. Phenomenological applications of this approach require matching of small-x contributions with contributions from medium and large x [10,[15][16][17][18][19][20]. To do this, in this work we further develop the parton branching Monte Carlo [19] implementation of the CCFM evolution equation [21,22], which we include in the herafitter program [1,23]. The TMD gluon distribution at the initial scale q 0 of the evolution is determined from fits to DIS data, including charm production.
We perform fits to measurements of the F 2 structure function [1] in the range x < 0.005, Q 2 > 5 GeV 2 and to measurements of the charm structure function F (charm) 2 [8] in the range Q 2 > 2.5 GeV 2 . We obtain good fits to F 2 and F (charm) 2 , and we make a determination of the TMD gluon density (as well as of the charm mass m c and of Λ QCD , or the strong coupling α s ) based on these. We find that the best fit to F (charm) 2 gives χ 2 per degree of freedom χ 2 /ndf ≃ 0.63, and the best fit to F 2 gives χ 2 /ndf ≃ 1. 18. The method allows one to assign experimental and theoretical uncertainties to pdfs. We give results for these different kinds of pdf uncertainties. We also carry out an application of the TMD gluon density determined from HERA data fits to LHC physics, by computing predictions for W -boson + jet production in proton-proton collisions.
The paper is organized as follows. In Sec. II we summarize the main elements of the approach based on high-energy factorization and evolution, and discuss a few details on its implementation. In Sec. III we describe the fits to charm and DIS precision data. We discuss F 2 , F (charm) 2 , the TMD gluon density determination and associated uncertainties. In Sec. IV we illustrate the use of TMD gluon density at the LHC. We give conclusions in Sec. V.

II. FACTORIZATION AND EVOLUTION
In the framework of high-energy factorization [9] the deeply inelastic scattering cross section is written as a convolution in both longitudinal and transverse momenta of the TMD parton density function A (x, k t , µ) with off-shell partonic matrix elements, as follows Here x and Q 2 denote the Bjorken variable and photon virtuality, and the DIS cross sections σ j , with j = 2, L, are related to the customary DIS structure functions F 2 and F L by [24] where α is the electromagnetic coupling. The factorization formula (1) allows one to resum logarithmically enhanced x → 0 contributions to all orders in perturbation theory, both in the hard scattering coefficients and in the parton evolution, taking fully into account the dependence on the factorization scale µ and on the factorization scheme [11]. Explicit evaluations have been carried out through next-to-leading logarithmic accuracy [10][11][12] in DIS at x → 0. The physical origin of the logarithmically enhanced x → 0 corrections at higher loops lies in the contribution from regions not ordered in initial-state transverse momenta to the QCD multi-parton matrix elements.
The perturbative higher-loop corrections are generally found to be large at small x. Consider for instance the gluonic hard-scattering coefficient function C g 2 (x, α s , Q 2 /µ 2 ) for the DIS structure function F 2 [24]. Taking Mellin moments with respect to x, the perturbative expansion of C g 2 for N → 0, resummed to all orders in α s via Eq. (1), is given at scale µ 2 = Q 2 in the MS minimal subtraction scheme by [11] where α s = α s C A /π, C A = 3, (3) to next-to-leading-logarithmic higher-loop corrections α 2 s (α s ln x) k−1 . The first two terms in Eq. (4) are the leading-order (LO) [25] and next-to-leading-order (NLO) [26] small-x contributions to C 2 . The next two terms are the three-loop and fourloop small-x contributions. The three-loop coefficient agrees with the complete next-tonext-to-leading-order (NNLO) calculation [27]. The three-loop and four-loop terms are logarithmically enhanced compared to lower orders. Moreover, their numerical coefficients are significantly larger than the one-loop and two-loop ones. Analogous results were obtained in [11] for the coefficient function C L of the longitudinal structure function, and confirmed through three loops by the NNLO calculation [28].
Given these results, there is little theoretical justification for treating small-x DIS by truncating the perturbative expansion to fixed NLO (or NNLO) level. Thus, although phenomenologically successful in giving very good fits to structure function data, fixed-order perturbative approaches are theoretically disfavored, and cannot be expected to describe the physics of the scaling violation in the region of low values of x, where transverse-momentum ordering does not apply.
The motivation of this work is to take a quantitative step toward going beyond fixed-order phenomenology, and confront our results with the high-precision combined data [1,8]. The approach of this work is based on • including the hard-scattering kernelsσ j of Eq. (1), whose k t -dependence, once convoluted with the gluon unintegrated Green's function [9,11], controls the all-order resummation of the higher-loop terms α 2 s (α s ln x) k−1 in the structure functions F 2 and F L ; • including evolution of the transverse momentum dependent gluon density A by combining the resummation of small-x logarithmic contributions [29] with medium-x and large-x contributions to parton splitting [30].
This is done via a parton branching Monte Carlo implementation of the CCFM evolution equation [21,22], which we develop based on [19], and make available within the herafitter program [1,23].
In the remainder of this section we briefly recall the basic elements of this approach and give a few technical details. We start in Subsec. II A by recalling the main features of TMD matrix elements and evolution. In Subsec. II B we include the unintegrated valence quark density according to the method [31]. In Subsec. II C we discuss aspects of the numerical implementation. We give comments on the general approach in Subsec. II D.

A. TMD matrix elements and evolution
The DIS transverse-momentum dependent partonic cross sectionsσ j are evaluated in the second paper of [9] for j = 2 and in [20] for j = L, including the effects of finite quark masses. The small-x resummation for DIS structure functions based on these results is carried out in [11]. Let us describe the off-shell process γ(q) + g(k) → q(p) +q(p ′ ) in terms of lightcone momenta where, for any four-momentum, p µ = (p + , p − , p t ), with p ± = (p 0 ±p 3 )/ √ 2. In the high-energy limit we have In terms of lightcone and transverse momentum components, we have Q 2 = −q 2 t /(1 − y), x = Q 2 /(yS).
In Fig. 1 we plot the partonic kernel σ 2 , normalized to the value for k t = 0, as a function of k 2 t,g = −k 2 t at fixed x and Q 2 , for various values of quark masses. As shown in [9,20], the large-k t tail of the kernels in Fig. 1, once it is folded with unintegrated gluon distributions [21,29,32], is responsible for the logarithmically enhanced higher-loop contributions at high energy. FIG. 1. Transverse momentum dependence of the partonic kernel σ 2 for the off-shell process γ(q) + g(k) → q(p) +q(p ′ ) at different values of Q 2 and quark masses. We set x = 10 −2 .
For evolution of the TMD gluon density, we require that in the limit x → 0 this evolves with the full BFKL anomalous dimension [29]. The CCFM evolution equation [21,22] is an exclusive equation which satisfies this property (see e.g. appendix C of [21], section 7 of the first paper in [22], section III of [33]) and, in addition, includes finite-x contributions to parton splitting. It incorporates soft gluon coherence for any value of x. The evolution equation for the TMD gluon density A(x, k t , p), depending on x, k t and the evolution variable p, reads The first term in the right hand side of Eq. (7) is the contribution of the non-resolvable branchings between the starting scale q 0 and the evolution scale p, and is given by where ∆ is the Sudakov form factor, and A 0 (x, k t , q 0 ) is the starting distribution at scale q 0 . The integral term in the right hand side of Eq. (7) gives the k t -dependent branchings in terms of the Sudakov form factor ∆ and unintegrated splitting function P. The gluons' average momentum does not change with p.
The CCFM evolution equation can be written in a differential form, best suited for the backward evolution approach adopted in the Monte Carlo generator [34,35], as where the splitting variable x ′ is given by and φ is the azimuthal angle of q t . The Sudakov form factor ∆ s is given by with α s = C A α s /π = 3α s /π. The splitting functionP g (z i , q i , k ti ) for branching i is given by [14] P where ∆ ns is the non-Sudakov form factor defined by Quark masses are treated in the fixed flavor number scheme. We include the two-loop running coupling α s , and we apply the kinematic consistency constraint [36,37] in the g → gg splitting function [15], given by [36] The evolution (7) of the TMD density implies that regions of transverse momenta below q 0 can be reached. In this region the branching is performed purely by the Sudakov form factor.

B. Unintegrated valence quark density
Previous determinations of parton distributions based on the CCFM evolution have included only the gluon density [13,14]. In this work we include valence quarks using the method of [31]. We consider the branching evolution equation at the transverse-momentum dependent level according to where p is the evolution scale. The quark splitting function P is given by withᾱ s = C F α s /π. In Eqs. (14),(15) the non-Sudakov form factor is not included, unlike the CCFM kernel given in the appendix B of the first paper in [22], because we only associate this factor to 1/z terms. The term xQ v 0 in Eq. (14) is the contribution of the non-resolvable branchings between starting scale q 0 and evolution scale p, given by where ∆ s is the Sudakov form factor, and the starting distributions at scale q 0 are parameterized, using the CTEQ 6.6 [38] u and d valence quark distributions, as In the numerical calculations that follow we will take σ 2 = q 2 0 /2. For every p we ensure that the flavor sum rule is fulfilled.

C. Numerical implementation
CCFM evolution cannot easily be written in an analytic closed form. For this reason a Monte Carlo method is employed, based on [19]. The Monte Carlo solution is however time-consuming, and cannot be used in a straightforward manner in a fit program. Here we proceed as follows. First a kernelÃ (x ′′ , k t , p) is determined from the Monte Carlo solution of the CCFM evolution equation, and then this is folded with the non-perturbative starting distribution A 0 (x), following the convolution method introduced in [39]: The kernelÃ incorporates all of the dynamics of the evolution, including Sudakov form factors and splitting functions. It is determined on a grid of 50 ⊗ 50 ⊗ 50 bins in x, k t , p. The binning in the grid is logarithmic, except for the longitudinal variable x where we use 40 bins in logarithmic spacing below 0.1, and 10 bins in linear spacing above 0.1. The calculation of the cross section according to Eq. (1) involves a multidimensional Monte Carlo integration which is time consuming and suffers from numerical fluctuations. This cannot be employed directly in a fit procedure involving the calculation of numerical derivatives in the search for the minimum. Instead the following procedure is applied: Here, firstσ(x ′ , Q 2 ) is calculated numerically with a Monte Carlo integration on a grid in x for the values of Q 2 used in the fit. Then the last step in Eq. (19) is performed with a fast numerical gauss integration, which can be used in standard fit procedures.

D. Comparison with other approaches and outlook
The approach described above, which we are going to confront in the next section with high-precision DIS measurements, can be compared with other approaches to DIS data in the literature. On one hand, it can be contrasted with descriptions of data based on the DGLAP equation [30] at fixed perturbative order, e.g. NLO or NNLO (see for instance [2] and references therein). As recalled at the beginning of this section, these descriptions, however successful phenomenologically, have little theoretical justification at small x, due to the structure of the perturbative expansion for x → 0. The need to go beyond fixed-order truncations of perturbation theory leads us to employ a TMD formalism. In particular, in the approach of this paper the transverse momentum dependence of the gluon density arises both from perturbative and from nonperturbative processes. Both the kernel and the initial condition of the evolution equation are k t -dependent. These different sources of k t -dependence have distinct physical effects, for instance on the angular distributions of associated jet final states in DIS, as analyzed in detail in [17]. This feature can be contrasted with TMD approaches focusing on nonperturbative k t -dependence, see e.g. [5].
On the other hand, the approach of this paper can be compared with approaches based on the BFKL equation [29]. Recent fits to DIS data have been performed in this context [40][41][42]. Compared to these works, the main theoretical inputs in the present paper are • the use of TMD matrix elements which can be directly related with the resummation of DIS coefficient functions as in Eq. (4), and • the use of the CCFM evolution for the TMD parton density rather than the BFKL evolution. This results in distinctive properties due to soft gluon coherence of the final states contributing to DIS [15].
A further, distinctive feature of the framework employed in this paper is that the gluon distribution obtained from DIS fits can be directly used to make predictions for final states at the LHC, as we do for example in Sec. IV for W -boson production associated with jets. Our approach relies on perturbative factorization theorems, which classify higher-order corrections according to the logarithmic hierarchy based on high Q 2 and low x. For this reason in the next section we will apply this approach to F 2 structure function measurements in a range Q 2 > Q 2 , x < x, where we choose Q 2 = 5 GeV 2 , x = 5 · 10 −3 . For asymptotically small x one expects the operator product expansion to break down and DIS to become dominated by nonperturbative physics. Thus the low-Q 2 region could require methods beyond the ones applied in this work, see e.g. [43,44]. On the other hand, the evolution approach used in this work may be supplemented with nonlinear corrections [45][46][47] to describe aspects of parton saturation [48][49][50][51]. It will therefore be of interest to investigate the extension of the work presented in this paper to low Q 2 . The inclusion of data at higher Q 2 , relaxing the low-x kinematic cut, will constitute a further development of the TMD formulation. High-x theoretical issues, including TMD quark distributions, are discussed e.g. in [52][53][54]. Analyses of DIS data over the whole available range in x and Q 2 in terms of TMD pdfs will be relevant for the calculational program [55] of off-shell hard cross sections.

III. FITS TO DIS PRECISION DATA
The fit to the HERA structure function measurements is performed by applying the herafitter program [1,23] to determine the parameters of the starting distribution A 0 at the starting scale q 0 . We perform fits by using two possible parameterizations of A 0 : the five-parameter form and the three-parameter form As in Eq. (17), we take σ 2 = q 2 0 /2. The parameters N, B, C, D, E (resp. N, B, C) in Eq. (20) (resp. Eq. (21)) are determined by fitting the high-precision structure function measurements [1] in the range x < 0.005 and Q 2 > 5 GeV 2 , and the charm production measurements [8], which are in the range Q 2 > 2.5 GeV 2 . The results presented here are obtained with the herafitter package by treating the correlated systematic uncertainties separately from the uncorrelated statistical and systematic uncertainties. We fit the charm leptoproduction data [8] and the inclusive structure function data [1] based on high-energy factorization and CCFM evolution as described in Sec. II. In particular, we also include two-loop running coupling, gluon splitting and consistency constraint as in Subsec. II A, and, in addition to the gluon-induced process γ * g * → qq, the contribution from valence quarks is included via γ * q → q as in Subsec. II B by using a CCFM evolution of valence quarks [31].
To obtain a reasonable fit to structure function data, we vary the starting scale q 0 , the QCD scale Λ QCD and the charm quark mass m c . The results in Table I Table I reports the values of χ 2 per degree of freedom for the best fit to the charm structure function F (charm) 2 [8] in the full data range Q 2 > 2.5 GeV 2 , and to the inclusive structure function F 2 [1] in the data range x < 0.005, Q 2 > 5 GeV 2 and for a combination of both, in the cases of the three-parameter fit (21) and five-parameter fit (20). The best-fit χ 2 /ndf is below 1 for the charm structure function, while it is around 1.18 for the inclusive structure function. This is in accord with the expectation that charm production is dominated by the gluon distribution, while the inclusive structure function receives significant contributions from quark channels, for which an improved treatment at unintegrated level is needed.
Also, the gluon density determined from charm measurements turns out to be steeper at small x than that determined from the inclusive structure function. The resulting tension between the two fits may be related to the the fact that the starting scale q 0 is not far from the charm threshold mass 2m c . The χ 2 /ndf for the fit to both F (charm) 2 and F 2 is 1.43 for the three-parameter fit, and is not significantly changed by using the five-parameter form. Fig. 2 shows the description of the charm leptoproduction measurements [8] and inclusive structure function measurements [1], by the individual fits to F (charm) 2 and F 2 and by the combined fit. Plotted are the reduced cross sections defined in [1,8].
In Fig. 3 we show the results of a scan in the charm mass, by plotting the values of χ 2 for the fit to charm and inclusive measurements as a function of the charm quark mass. The minimum is reached for m c = 1.45 GeV. . The values of χ 2 for the fit to DIS high-precision data including charm leptoproduction [8] versus the charm quark mass.

B. Unintegrated TMD gluon density
Here we present two sets of unintegrated pdfs determined from the fits to high-precision DIS measurements described in the previous subsection: JH-2013-set1 is determined from the fit to inclusive F 2 data only; JH-2013-set2 is determined from the fit to both F (charm) 2 and F 2 data. The unintegrated TMD gluon density is shown in Figs. 4 and 5 at different evolution scales, versus the longitudinal momentum fraction x and versus the transverse momentum k t . The results are compared with the older parton distribution set A0 [35], which did not use the precision data and did not include two-loop running coupling, kinematic consistency constraint, nonsingular terms in the gluon splitting function, and unintegrated valence quark density.
We see from the evolution to scale p 2 equal to the Z-boson mass in Fig. 4 that the fit to the DIS high-precision measurements gives significant differences in the TMD gluon density compared to earlier determinations, especially in the region of medium to low k t . For the lower p 2 scale in Fig. 5 the differences are less pronounced but still important especially for small values of x.

C. Uncertainties on the TMD gluon density
In this subsection we consider separately experimental and theoretical uncertainties of the TMD parton distributions.
Experimental uncertainties are obtained within the herafitter package from a variation of the individual parameter uncertainties, following the procedure described in [56] applying ∆χ 2 = 1. These result in 10 to 20 percent gluon uncertainty for medium and large x. The experimental uncertainties on the gluon at small x are small (much smaller than those obtained in standard fits based on integrated pdfs), since only the gluon density is fitted. The   Fig. 7 for lower evolution scale p 2 = 25 GeV 2 .
Next we consider theoretical uncertainties. The first such kind of uncertainty is the dependence on the starting scale q 0 for gluon density evolution. In Figs. 6 and 7 the dotted blue curves show the effect on the gluon distribution from variation in the starting scale q 0 . These uncertainties are small at small x, while they become very large at large x because in this region, since we fit F 2 in the range x < 0.005 and Q 2 > 5 GeV 2 , there is little constraint from data.
We also consider theoretical uncertainties on the TMD gluon density from variation of the factorization scale and renormalization scale. This approach is different from that usually followed in determinations of ordinary, collinear pdfs from fixed-order perturbative treatments [2]. In this case, no uncertainty on the pdfs is considered from scale variation. Only when computing predictions for any specific observable the theoretical uncertainty on the predictions is estimated by scale variation. In our approach we are interested to estimate the uncertainty from varying scales in the theoretical calculation used to determine the pdf. In Figs. 6 and 7 the renormalization scale (blue dashed curves) and in the factorization scale (yellow band) are varied by a factor of 2.

D. Integrated parton distributions
For a cross check with the integrated pdfs we now compute the integral over transverse momenta of the TMD parton distributions. In Fig. 8 we plot the results for gluon and valence quark distributions, obtained from the set JH-2013-set1 of this paper (and also, for comparison, from the gluon in the older set A0 [35]), at two different evolution scales. For comparison we plot the ordinary, integrated distributions obtained from the NLO-DGLAP CTEQ 6.6 [38] fit. We observe good agreement for the integral of valence quark distributions at low scales while differences arise from the different evolution at larger scales. For the gluon case, the difference between the integral of TMD and CTEQ reflects the shuffling of flavor singlet contributions between sea quark and gluon in the two formalisms.
We conclude this section by noting that the CCFM evolution kernel can be approximated to one loop by using collinear ordering [16,17]. This constitutes the DGLAP limit of the evolution equation. If we perform fits to the high-precision F 2 data by using Eqs. (1), (7) in the one-loop approximation mode we find that this approximation is unable to give a good fit based on the TMD gluon only, χ 2 /ndf ∼ 6. We interpret this as a check on the consistency of the physical picture, signaling the need for introducing quark-initiated processes in the collinear framework.

IV. TMD GLUON DENSITY AT THE LHC
The TMD parton distributions determined in Sec. III from fits to the high-precision DIS data can be used to make predictions for hadron-hadron collider processes.
An example is the Drell-Yan vector boson production. We here consider W -boson production in association with jets. This process is important both for standard model physics and for new physics searches at the LHC. In particular it is relevant to studies of parton distribution functions and of Monte Carlo event generators [2], including signals of multi-parton interactions, for which W + 2 jets is a classic channel [57,58].
To compute predictions for W -boson + jets final states [59], we use the CCFM gluon and valence quark distributions determined in the previous section, convoluted with highenergy matrix elements [60,61] with off-shell partons [62,63] for weak boson production. We present results for the inclusive jet multiplicity distribution and leading jet transverse momentum spectrum. The results (obtained with the Rivet -package [64]) are shown in Fig. 9 along with the ATLAS measurements [65].
The solid red curves in Fig. 9 are the predictions from JH-2013-set2, with the blue band corresponding to the pdf uncertainty. Both the jet multiplicity and the transverse momentum are reasonably well described by the predictions. Further discussion on W + jets will be given in [59]. The production of final states with W boson and multiple jets at the LHC receives sizeable contributions from large separations in rapidity between final-state particles. However, the cross sections computed in Fig. 9 are not dominated by very small values of x. As a result, the uncertainty band due to the uncertainties in the TMD pdfs is significant. A comparison with NLO-matched results and corresponding uncertainties is presented in [59]. It is conceivable that combining pp measurements on vector boson production with the DIS measurements analyzed in this paper may help to constrain TMD pdfs especially at medium to large values of x.

V. CONCLUSIONS
In this work we have performed the first determination of the TMD gluon density function from high-precision DIS measurements, including experimental and theoretical uncertainties.   9. Predictions for W + jets production using the unintegrated TMD gluon density JH-2013-set2: (left) inclusive jet multiplicity; (right) leading jet transverse momentum spectrum. The experimental data are from [65]. The yellow band is the experimental uncertainty. The blue band is the theory uncertainty.
We have presented fits, based on QCD high-energy factorization and CCFM evolution, to HERA charm-quark leptoproduction data for the structure function F (charm) 2 [8] in the range Q 2 > 2.5 GeV 2 , and to HERA F 2 structure function data [1] in the range x < 0.005 and Q 2 > 5 GeV 2 . In this approach the charm structure function can be regarded as a physical probe of the unintegrated TMD gluon density. We fit the combined HERA charmquark data [8] over the whole kinematic range of the measurement, and obtain that the best fit gives χ 2 per degree of freedom χ 2 /ndf ≃ 0.63. The inclusive F 2 structure function involves both gluon-density and quark-density channels. We fit the combined HERA F 2 data [1] in the kinematic range x < 0.005, Q 2 > 5 GeV 2 , and obtain that the best fit gives χ 2 /ndf ≃ 1.18. Despite the restricted kinematic range, the great precision of the data provides a highly nontrivial test of the approach. We find a good fit to both charm-quark and inclusive data. Based on this, we make a determination of the TMD gluon density (as well as of the QCD scale Λ QCD and the charm mass m c ) and present new unintegrated pdf sets, JH-2013. As a result of the high-precision data, the JH-2013 distributions differ significantly from earlier sets. We also present experimental and theoretical uncertainties associated with the TMD pdfs. We compute predictions based on the TMD pdfs for Wboson plus jets production at the LHC, and find that the results compare well with the measurements [65] of jet multiplicities and transverse momentum spectra within the pdf uncertainties.
The approach of this work is based on the use of transverse momentum dependent matrix elements and evolution. Both the kernel and the initial condition of the evolution equation are k t -dependent. The transverse momentum dependence of the gluon density arises from both perturbative and nonperturbative processes. The physical picture of DIS scaling violation underlying this approach differs from that of finite-order perturbative QCD fits, e.g. at the NLO level, because it takes into account corrections to the collinear ordering in the initial state evolution to all orders in the QCD coupling α s . On the other hand, it also differs from BFKL evolution because it takes into account, for any x, color coherence associated with soft multi-gluon emission. In this paper we have developed a parton branching Monte Carlo implementation of the CCFM evolution equation and we have included it in the herafitter program [1,23].
The choice of the kinematic range for the F 2 data considered in this paper is motivated by the fact that our approach relies on perturbative factorization theorems, which classify higher-order corrections according to logarithmic hierarchy based on high Q 2 and low x. However, we note that the choice made in this paper is conservative, and the physical picture lends itself to extensions to lower Q 2 and higher x. On one hand, this picture goes beyond DGLAP by including non-collinear emissions to all orders in α s . On the other hand, it goes beyond BFKL by including large-x terms according to the CCFM prescription. Furthermore, the study of W -boson + jets made in this paper suggests that pp measurements of vector boson production at the LHC may be used to extend experimental investigations of TMD parton density functions.