Transverse momentum dependent (TMD) parton distribution functions: status and prospects

We provide a concise overview on transverse momentum dependent (TMD) parton distribution functions, their application to topical issues in high-energy physics phenomenology, and their theoretical connections with QCD resummation, evolution and factorization theorems. We illustrate the use of TMDs via examples of multi-scale problems in hadronic collisions. These include transverse momentum q_T spectra of Higgs and vector bosons for low q_T, and azimuthal correlations in the production of multiple jets associated with heavy bosons at large jet masses. We discuss computational tools for TMDs, and present an application of a new tool, TMDlib, to parton density fits and parameterizations.


I. INTRODUCTION
Experimental information on "3-dimensional imaging" of hadrons, encoded in unintegrated, transverse momentum dependent (TMD) parton density and parton decay functions, comes at present from two main sets of experimental data: deep inelastic scattering (DIS) at high energy, and low-q T Drell-Yan (DY) and semi-inclusive DIS (polarized and unpolarized). In each of these two cases, QCD factorization theorems allow one to relate physical, observable cross sections to TMD parton distributions via perturbatively calculable kernels. These theorems provide the theoretical basis for determining TMD distributions from experimental measurements. They are also essential to formulate and apply methods of perturbative resummation at all orders in the QCD coupling to a large variety of observables in high-energy hadronic collisions. Examples include processes both at the Large Hadron Collider (LHC) and at fixed-target experiments.
This article is based on workshops devoted to these topics held at the University of Antwerp in 2014 1 and provides a concise status report of this field. The purpose of this article is to give a non-technical introduction to the motivations for experimental and theoretical studies of TMDs; to illustrate this with specific examples of application of TMDs to topical issues in high-energy physics phenomenology; to point to future directions of development.
In particular, we examine implications of two sets of QCD factorization theorems based on TMD parton distribution functions: low-q T factorization for heavy particle spectra (including vector bosons, Higgs bosons, heavy flavors) and high-energy factorization. We focus on production processes in hadronic collisions in two limits: i) q T → 0 for fixed invariant mass, and ii) √ s → ∞ for fixed momentum transfer. We illustrate this with examples on transverse momentum spectra and angular correlations at the LHC for Drell-Yan and Higgs boson production and associated multi-jets. We survey computational tools which are being developed to treat the physics of TMDs. We present in particular an application of a new tool, TMDlib, to TMD parton densities based on fits and parameterizations including QCD evolution. The paper is organized as follows. In Sec. II we motivate the use of TMDs. In Sec. III we discuss their role in the physics of large transverse momenta. Sec. IV summarizes experimental prospects and theory developments. Secs. V and VI illustrate the status of fits and parameterizations for TMD parton distributions and of TMD Monte Carlo tools. Final remarks are given in Sec. VII.

II. WHY TMDS
Transverse momentum dependent parton distributions encode nonperturbative information on hadron structure, including transverse momentum and polarization degrees of freedom, which is essential in the context of QCD factorization theorems for multi-scale, non-inclusive collider observables. A classic example is given by Drell-Yan hadroproduction of electroweak gauge bosons. Fig. 1 [1] shows the differential cross section for Z-boson production in pp collision at the LHC as a function of the Z-boson transverse momentum q T , in the lepton pair's invariant mass range 60 GeV < M < 120 GeV. In the spectrum of Fig. 1 we distinguish the high-q T region, the peak region, and the low-q T region.
In the high-q T region the cross section is expected to be well represented by an evaluation of the partonic Z-boson cross section to finite order in QCD perturbation theory (leading-order (LO), next-to-leadingorder (NLO), and so forth), combined with factorization in terms of ordinary (collinear) parton distribution functions (pdfs). On the other hand, if this theoretical framework is applied to the region of decreasing q T it will not be able to describe the approach to the peak region in Fig. 1 (q T ≈ O(10 GeV)) nor the turn-over region (q T ≈ O(1 GeV)). Rather, the cross section predicted from any finite order of perturbation theory, convoluted with ordinary parton distributions, will diverge as q T decreases. The reason for this is that the physical behavior of the Z-boson spectrum near the peak region and below [2,3] is controlled by multi-parton QCD radiation, which is not well approximated by truncating the QCD perturbation series to any fixed order but rather requires methods to resum arbitrarily many parton emissions, viz., scattering amplitudes with an infinite number of real and virtual insertions of soft gluons.
This can be accomplished in a systematic manner via a generalized form of QCD factorization [4][5][6] which now involves quark distribution functions that, unlike the ordinary ones, explicitly depend on transverse momentum and polarization (TMD pdfs). Such TMD pdfs obey evolution equations [6][7][8] which generalize the ordinary renormalization-group evolution equations of collinear pdfs. These evolution equations, once combined with the TMD factorization of the physical cross section, allow one to resum logarithmically  enhanced contributions in the ratio M/q T to the perturbation series expansions for the physical observables to all higher orders in the QCD coupling. It is only after this generalized factorization analysis -going beyond the collinear factorization -is carried through that the physical behavior of the Z boson spectrum observed in Fig. 1 can be predicted. A second example concerns the rise of proton's structure functions at small longitudinal momentum fractions. Since in pp collisions the product of initial-state longitudinal fractions scales like 1/s at fixed momentum transfer, where s is the squared centre-of-mass energy, as we push forward the high-energy frontier more and more events at small longitudinal fractions contribute to processes probing short-distance physics. Many hard-production cross sections at the LHC receive sizeable contributions from proton's structure functions in this region. As parton longitudinal momenta become small, the fraction of momentum carried by transverse degrees of freedom becomes increasingly important. Fig. 2 shows the proton's gluon density resulting from global fits [9] to hadronic collision data, performed at LO, NLO, NNLO [10][11][12] of perturbation theory, as a function of the longitudinal momentum fraction x for different values of the evolution mass scale Q 2 . In the low-x regime the perturbative higher-order corrections to structure functions are large, and the gluon pdf uncertainty is large. The strong corrections at low x come from multiple radiation of gluons over long intervals in rapidity [13,14], in regions not ordered in the gluon transverse momenta p T , and are present beyond NNLO to all orders of perturbation theory [15,16]. The theoretical framework to resum these unordered multi-gluon emissions is a generalized form of QCD factorization [17,18] in terms of TMD pdfs. Analogously to the Drell-Yan case discussed earlier, the TMD pdfs obey a suitable set of evolution equations [19][20][21], appropriate to this kinematic region. These provide another generalization, valid in the high-energy limit, of the ordinary renormalization-group evolution. The TMD factorization in this case allows one to resum logarithmically enhanced corrections in the ratio √ s/Q to all higher orders in the QCD coupling.
Besides the above examples of Drell-Yan and structure functions, TMD factorization theorems apply to a wide variety of processes at the LHC. In particular, with extensive measurements of Higgs boson production at the LHC Run II, a new set of QCD processes becomes available in which the Higgs boson acts as a color-singlet, pointlike source (in the heavy top limit) which couples to gluons. This is to be contrasted with Drell-Yan and deep-inelastic scattering cases, based on weak and electromagnetic currents providing color-singlet pointlike sources coupled to quarks. This opens up the possibility of a new program of precision QCD measurements in gluon fusion at high mass scales in the LHC high-luminosity runs [22,23].
Analogously to the case of vector bosons in the example of Fig. 1, theoretical predictions for the Higgsboson production differential spectrum over the whole range in transverse momenta accessible at the LHC require generalized QCD factorization, based on initial-state gluon distributions that include polarization and transverse-momentum degrees of freedom. Compared to the vector boson case, however, new features arise which are associated with the role of gluon polarizations in gluon-gluon scattering.
More precisely, in the high-energy limit √ s m H the Higgs boson production from gluon fusion is dominated by a single eikonal gluon polarization [25]. The contribution of this polarization depends on the gluon transverse momentum and can be rewritten in terms of the high-energy projection operator defined in [16]. A complete set of operators for polarization dependent and transverse momentum dependent gluon distributions is given in [26]. In the region of low Higgs-boson transverse momenta, q T m H , the contributions of polarized gluons to the Higgs spectrum have been studied both perturbatively [27][28][29][30][31] and nonperturbatively [24,[32][33][34][35]. An example is shown in Fig. 3 [24], where the unpolarized and linearly polarized gluon distributions contributing to the Higgs boson spectrum at small q T are plotted as a function of transverse momentum. The presence of polarized gluon components (even in unpolarized beams) characterizes gluon fusion processes and has no analogue in the Drell-Yan case. In particular the component in the right hand side plot of Fig. 3 is a gluon TMD distribution with double spin flip (see Table II ahead, top right corner). From the point of view of perturbative power counting, double spin flip effects start to contribute to the Higgs q T spectrum at the NNLO (but may contribute earlier in more complex, less inclusive observables associated with Higgs production). Detailed measurements of Higgs boson final states will allow the QCD dynamics of polarized gluons and their correlations to be explored experimentally for the first time.
For both the Drell-Yan and Higgs cases, in addition to the inclusive spectra an extensive experimental program at the LHC is devoted to the associated production of heavy bosons with jets. The region in which the boson and leading jet are nearly back-to-back presents features comparable to the discussion given above for the low-q T part of the inclusive spectra. For instance, a study of TMD gluon contributions to Higgs + jet final states in which the imbalance between the boson and leading-jet transverse momenta is small The transverse momentum dependence of the unpolarized (left) and linearly polarized (right) gluon distributions [24] contributing to the gluon-fusion Higgs production spectrum. The results are plotted for evolution scale Q = 20 GeV and longitudinal momentum fraction x = 0.01, and for different values of the nonperturbative parameters discussed in [24].The red and green bands around each curve correspond to variations by factor 2 of the resummation scale and rapidity scale in the calculation [24].
is reported in Fig. 4 [36], showing the boson-jet pair's transverse momentum distribution and azimuthal asymmetries.

FIG. 4:
Theoretical predictions [36] for the transverse momentum distribution (left), cos 2φ asymmetry (middle) and cos 4φ asymmetry (right) in Higgs boson + jet production at small transverse momenta qT of the Higgs + jet pair.
Here K ⊥ represents the average of the Higgs and jet transverse momenta, and the shaded blue areas represent the range of the asymmetries as K ⊥ varies from 0 to ∞.
The case of associated boson + jet production when the imbalance between the boson and leading-jet transverse momenta is not small, on the other hand, probes the physics of final states with multiple jets. The role of TMD parton distributions in scenarios with high jet multiplicity is discussed in the next section, and serves to illustrate the connection of TMDs with the kinematic region of large transverse momenta.
An extension of the methods discussed above for Drell-Yan and Higgs production applies to the transverse momentum spectra of heavy flavor pairs, e.g. top quarks. Unlike the case of color-singlet currents coupled to quarks (as in Drell-Yan production) or gluons (as in Higgs boson production in the heavy top limit), heavy-quark pair production constitutes a composite non-pointlike probe, containing color-charged particles in the lowest-order final state and receiving contribution from both quark and gluon TMD channels. Color correlations over long timescales between initial and final states will break factorization in the region of very small transverse-momentum imbalance of the pair [40][41][42][43][44][45][46][47]. Studies of this region and the interplay of perturbative and nonperturbative contributions will help understand quantitatively these effects.
Another area for applications of TMDs concerns single spin asymmetries and azimuthal asymmetries    [39] shows low-energy measurements [37,38] of the Sivers transverse single spin asymmetry along with results of the fit [39]. For hadron's transverse momenta sufficiently small compared to the virtuality scale Q of the deep inelastic (or Drell-Yan) process, spin asymmetries obey TMD factorization formulas of the same kind [6] discussed above for the unpolarized case of low q T Drell-Yan. A combined understanding of current high-energy unpolarized measurements and low-energy spin asymmetry measurements is important for the planning of future polarized collider [51,52] and fixed-target [53,54] experiments. We conclude this section by presenting the full leading-twist set of polarization dependent and transverse momentum dependent parton densities in a spin-1/2 hadron. These are shown in Table I and Table II, for the quark [55,56] and gluon [26,57] cases respectively, including the distributions in unpolarized hadrons (top rows), longitudinally polarized hadrons (middle rows), transversely polarized hadrons (bottom rows). (See [58][59][60][61][62][63][64] for slightly different classifications.) Gauge-invariant operator definitions may be given for each of the TMD distributions in terms of nonlocal operator combinations, in which appropriate Wilson-line gauge links are associated with quark and gluon fields [6,[65][66][67][68][69]. Operator definitions are instrumental in analyzing both factorization and potential sources of factorization breakdown, and in setting up lattice calculations [70][71][72][73] of parton distributions.

III. TMDS AND LARGE TRANSVERSE MOMENTA
Unlike the low-q T Drell-Yan factorization theorem [4][5][6] and its extensions for gluon fusion processes, the high-energy factorization theorem [16][17][18]74] is valid for arbitrarily large momentum transfer. It is based on the high-energy expansion √ s → ∞ and can be applied in the ultraviolet region of high q T . It allows one, for example, to obtain the structure of logarithmic scaling violations in DIS at high energy (see [10][11][12]) and to resum logarithmic corrections of higher order in α s to Higgs and top quark production cross sections  6: Total transverse energy HT distribution in final states with W -boson + n jets at the LHC, for n ≥ 1 (left), n ≥ 2 (center), n ≥ 3 (right). The purple, pink and green bands correspond to the different methods described in [78] to estimate theoretical uncertainties. The experimental data are from [89], with the experimental uncertainty represented by the yellow band.
(see [75][76][77]). In this section we apply this theorem to discuss the role of TMDs in the region of perturbative transverse momenta, in particular in the high-q T part of the Drell-Yan spectrum in Fig. 1.
The basic observation is that the LHC kinematics leads to copious production of final states in which a high-q T vector boson recoils against multiple hard jets. Ref. [78] studies W -boson + n jets final states using TMD high-energy factorization [18]. The motivation for this is twofold: a) kinematical: it has recently been pointed out [79][80][81] that collinearity approximations, once combined with energy-momentum conservation constraints, give rise to longitudinal momentum shifts and sizeable showering corrections in the Monte Carlo algorithms used to simulate multi-jet final states at the LHC; b) dynamical: it has long been known [82][83][84] that, when the picture of multi-jets from finite-order perturbative matrix elements matched with collinear parton showers is pushed to higher and higher energies, new effects arise in jet multiplicity distributions and angular correlations due to soft but finite-angle multigluon radiation. Both these kinematical and dynamical effects can be taken into account by a TMD treatment of QCD parton shower evolution [84].
To achieve this, Ref. [78] uses the exclusive formalism of CCFM evolution equations [82,85,86] implemented in [87]. The TMD pdfs to which evolution is applied are determined from fits to the precision DIS data [88]. By evolving these TMD pdfs up to the scale of W + jets and coupling them with appropriate, perturbatively calculated high-energy matrix elements, one obtains predictions for W -boson + n jets observables. Fig. 6 shows the total transverse energy H T distribution in final states with W -boson + n jets, with n = 1, 2, 3, at the LHC. For comparison the experimental measurements [89] (jet rapidity |η| < 4.4, jet transverse momentum p T > 30 GeV) are plotted. The uncertainty bands on the theoretical predictions are described in [78], and largely reflect uncertainties on TMDs determinations, estimated according to three different approaches corresponding to the three color bands.
The TMD high-energy factorization predicts azimuthal correlations in the W + multi-jet final states. Fig. 7 shows results for the azimuthal correlation between the two leading jets, along with the transverse momentum of the third jet.
Current limitations of the approach described above and ongoing improvements are discussed in [78,90,91], and include in particular the treatment of TMD quark density distributions, and the accuracy of determinations of the gluon density distribution over the whole range of longitudinal momentum fractions x relevant to the LHC kinematics. The results in Figs. 6 and 7 are however encouraging, and sufficiently general, in the context of approaches that aim to go beyond fixed-order perturbation theory and appropriately   FIG. 7: Azimuthal correlation of the two leading jets associated with W -bosons, for pT > 20 GeV (left) and pT > 30 GeV (center), and transverse momentum of the third jet (pT > 30 GeV) (right). The purple, pink and green bands correspond to the different methods described in [78] to estimate theoretical uncertainties. The experimental data are from [89], with the experimental uncertainty represented by the yellow band. take account of nonperturbative effects.
As TMDs describe nonperturbative transverse momentum dynamics in the hadron, they may provide a suitable framework not only for the factorization of the hard process but also to incorporate effects from soft particle production and multi-parton interactions [92,93].
It is worth noting that while for sufficiently inclusive observables in W + jets production calculations based on collinear parton showers matched with finite-order perturbative matrix elements describe measurements at Run I very well, this may not necessarily be the case for observables sensitive to the detailed structure of multi-parton emission [94,95]. For example, Fig. 8 [96] shows ATLAS measurements of the di-jet invariant mass associated with W production, compared with several Monte Carlo calculations. The comparison with the results from the NLO-matched calculation Blackhat + Sherpa [97] suggests that effects beyond NLO + collinear shower may set in for high invariant masses around and above 500 GeV. In this region of masses a similar behavior is observed in the comparison of experimental measurements with the Alpgen [98] Monte Carlo calculation. In Fig. 9 we also plot the di-jet invariant mass distribution from the approach [78]. For the physics program at Run II it is of much interest to examine the region of very large vector boson transverse momenta of order 1 TeV and higher. Fig. 10 [100] shows CMS measurements of the Z-boson p T in events with Z + 1 jet and Z + 2 jets at Run I. At the highest p T one may see dynamics setting in beyond the level of Madgraph [101] and Sherpa [102] multi-leg jet calculations matched with collinear showers, even supplemented with an NNLO k-factor.

IV. THEORETICAL DEVELOPMENTS AND EXPERIMENTAL PROSPECTS
This section gives a brief overview of ongoing theoretical developments and experimental prospects.
Evolution of TMDs and fits to physical cross sections. The above approaches to low-q T spectra which make use of TMDs currently employ, in practice, either approximate analytic (or semi-analytic) solutions of the evolution equations [6][7][8] or perturbative expansions of the TMDs in terms of collinear pdfs, or a combination of both. A different proposal has been put forward in [151] (TMDlib), based on global fits to experimental data to obtain TMD parton distributions at different evolution scales, and on using these to make predictions for physical quantities. This is similar in spirit (but different in its realization) to what is done in the case of collinear parton distributions. Theoretical predictions for physical cross sections which obey TMD factorization formulas could then be obtained by applying these formulas, using perturbatively calculable coefficients and appropriately evolved TMDs determined from fits to experiment. In this approach, unlike most current implementations of TMD formalisms, the nonperturbative dependence on longitudinal and transverse degrees of freedom is fully coupled, and can be entangled with the dependence on the evolution scale [151]. For phenomenological applications this can be important when for instance comparing theory with experimental measurements over a wide range in x and evolution scales.
Nonlinear evolution of the gluon TMD and Wilson line correlators. The conventional gauge-invariant operator definition of the gluon TMD [7,24,26,34,57,152] is distinct from the Weiszacker-Williams operator definition [66,[153][154][155][156] in terms of Wilson lines often used at x 1 (see also [157][158][159] for discussion of the operator definitions). Correspondingly, these gluon TMDs obey different rapidity evolution equations: in the moderate x ∼ 1 region one has linear double-logarithmic equations, while in the x 1 domain the non-linear single-logarithmic Balitsky-Kovchegov equation applies [160,161]. The relationship between these two regimes is examined in [162,163], where it is clarified that the non-linear small-x evolution transforms into linear rapidity evolution for the conventional gluon TMD. Refs. [164,165] consider applications to diffraction, and Refs. [166][167][168][169] to jets at large rapidity separations. Also, the evaluation of the complex combinations of Wilson lines entering the gluon TMD at small x calls for the development of a dedicated methodology. Essential improvement in the understanding and computation of correlators with Wilson lines can be achieved by the eikonal exponentiation methods [170][171][172], which enable the exact resummation of the diagrams presenting a given correlator as the exponent of series of the so-called web diagrams [173][174][175][176][177].
TMDs and generalized loop space. Renormalization properties of Wilson line correlators control the evolution of TMDs [6,[178][179][180][181]. In particular, the appearance of light-cone, or rapidity, divergences [6,182] in higher-loop corrections to the gauge-invariant correlators calls for a treatment of overlapping divergences, which can be achieved by the introduction of a soft subtraction factor [183][184][185][186][187][188]. The evolution of the gauge-invariant path-dependent TMDs with the light-like cusped Wilson lines can also be associated with the geometric evolution in the generalized space [189][190][191]. The differential shape variations of the underlying contours to the Wilson loops are formulated in terms of the Fréchet derivative [192,193] and the equations of motion in the loop space are dual to the energy and rapidity evolution of the TMDs having the same structure of the Wilson lines [194,195].
Non-universality and Wilson lines. Operator definitions of parton distribution functions in terms of quark and gluon fields involve nonlocal operator combinations. For collinear functions the nonlocality is along the lightcone, for TMDs it is along the lightfront involving also transverse separations. Unavoidably, therefore, additional gluonic fields minimally enter in the Wilson lines that are needed for an unambiguous gauge invariant description. The fact that these Wilson lines depend on the hard process brings in a calculable non-universality [65], which is a generalization of the sign flip between T-odd TMDs in going from SIDIS to DY [49,50]. Other examples where these effects appear are jet+jet or photon+jet final states in hadroproduction [196,197] as compared to Drell-Yan or ZZ production, and Higgs + jet final states as compared to Higgs production into colorless final states. Color entanglement can lead to further sources of non-universality affecting both TMD factorization [44] and collinear factorization [198,199]. A program is ongoing devoted to a careful analysis [67][68][69] of the possible operators that contribute to particular TMD structures followed by the study of their evolution.
Soft particle production and multi-parton interactions. As TMDs encode nonperturbative transverse momentum dynamics in the proton, one may ask whether they are relevant not only for factorization of hard processes but also for the understanding of soft particle production and, in particular, of the multi-parton interactions which are found to be needed at low to moderate transverse momenta for Monte Carlo simulations to describe experimental data on underlying events, particle multiplicities and spectra. Double parton interactions [223,224] including parton's transverse momentum dependence are investigated in [92,[225][226][227][228][229][230][231]. The role of parton's transverse momentum in the interpretation of energy flow measurements is discussed in [93,[232][233][234]. Implications for diffraction are considered in [235][236][237]. TMD effects in multi-parton correlations may be studied in upcoming measurements of charged particle multiplicities and spectra and underlying event at the LHC 13 TeV run.
Experimental prospects have been discussed for identifying TMD effects based on measurements of benchmark cross sections, both at the LHC and at lower energy experiments.
Drell-Yan lepton pair production and Drell-Yan plus jets. As discussed in the previous two sections, both the low-q T part of the spectrum and the high-q T part can be sensitive to TMD effects. Multi-differential measurements are especially important as one can access azimuthal correlations in the lepton + jet final states [78,202] which constitute distinctive TMD predictions. Comparison of Z + jet final states at small transverse momentum imbalance [196,197] with di-boson ZZ final states may shed light on color flow patterns which are eventually responsible for factorization breaking phenomena in hard processes sensitive to very low transverse momentum scales.
Higgs boson production and Higgs boson plus jets. Similar measurements to the Drell-Yan case, including differential cross sections, are relevant for gluon TMDs and QCD studies of polarized gluons and color correlations, once sufficient statistics is reached. Measurements of Higgs versus Drell-Yan at the same invariant mass may be used to reduce the influence of pile-up in the high-luminosity LHC runs [22]. The boson q T spectrum, final-state angular distributions and underlying event observables probe different aspects of the Higgs coupling to gluons [22,238].
Heavy flavor production. Measurements of top quark pair production spectra can provide comparable information to the previous two cases but with additional complexity due to the presence of color charges in the final state. The associated initial-state / final-state color correlations at small q T could be studied to examine factorization-breaking contributions in the region of very small transverse momenta [40][41][42][43][44], provided sufficient resolution can be reached. It will also be interesting to investigate kinematic effects of longitudinal momentum reshuffling in parton showers [79] at top quark scales. Similar studies can be done at lower mass scales with bottom and charm quarks.

V. WORKING WITH TMDS: FITS AND PARAMETERIZATIONS
The polarization dependent and transverse momentum dependent proton's parton densities, in the notation of [55,[67][68][69]259], are given in the Tables I and II in Sec. II. This scheme can be extended to spin-1 targets [260].
Most existing fits and parameterizations of these distributions may be grouped into categories which broadly correspond to three main areas of TMDs applications discussed in Sec. II: • Fits to vector boson q T experimental data in unpolarized Drell-Yan production [114-116, 118, 122, 123, 128, 129, 261-263] based on the low-q T TMD factorization [4,6], in some cases including extension to semi-inclusive DIS data [149,150,264].
• Fits to DIS structure function data [88, based on the high-energy TMD factorization [17,18] or on other approaches (e.g. saturation formalism) to high-energy DIS, in some cases including the precision measurements [286,287] and TMD pdf uncertainties [88,276].
For precision phenomenology it will be essential that results of fits and parameterizations are given in a portable form as a determination of TMD pdfs over a given kinematic range, appropriate to the theoretical method and experimental data used. A first step in this direction has been taken in [151]. The main point is that if results of fits to experimental data are used to provide TMD pdfs at different evolution scales, theoretical predictions for physical cross sections could then be obtained by using these pdfs in factorization formulas (or, eventually, in Monte Carlo event generators implementing these formulas). In [151] a library has been initiated, TMDlib, in order to unify and simplify the access of TMDs, along with a plotting tool, TMDplotter, for easier comparisons. Commonly used pdf sets are implemented in TMDlib, with the goal to provide a library of all available TMDs. In TMDlib pdfs are accessible in an easily callable way within the range of their applicability. The pdfs currently included range from TMD gluon densities obtained from fits to small-x DIS data based on high-energy factorization, to TMD gluon densities from fits based on saturation approaches, to TMD quark densities from parton-model fits to low-energy fixed-target data at large x and small k T . TMD fragmentation functions are not yet implemented, but are foreseen for the future.
[GeV] An example from TMDlib is shown in Fig. 11, plotting the transverse momentum dependence of valence quark distributions, at fixed values of x and renormalization scale p 2 , obtained from the fits [88,307].
In Fig. 12 we show results for gluon distributions [24,88]. The nonperturbative parameters of the distribution represented by the red curve are obtained from the fit [88] to DIS experimental data [286,287] while the distribution represented by the blue curve is not fitted but based on the model [24].

VI. WORKING WITH TMDS: MONTE-CARLO GENERATORS AND TOOLS
Inclusive or semi-inclusive hard cross sections can be calculated by convoluting parton density and decay functions with partonic cross sections. For a detailed description of the exclusive structure of the final states, on the other hand, event generators including parton showers and full hadronization are required.
In the collinear case, cross sections are computed with on-shell initial partons. For many processes, higher order calculations exist, and many of these are implemented in Monte Carlo (MC) simulation tools like Powheg [333,334], Mc@nlo [335], aMc@nlo [336], which combine next-to-leading order partonic calculations with parton showers and hadronization. These simulations all need a reshuffling of kinematic variables, after the parton shower is generated, in order to satisfy energy-momentum conservation, which can lead to significant kinematic shifts in the longitudinal momentum fraction x [79]. This is because transverse momentum is generated by the initial-state parton shower, which is not available when the hard scattering is computed. In certain phase space regions, these longitudinal shifts can affect the accuracy of the calculations significantly. Using TMDs, this kinematic reshuffling can be avoided from the beginning provided the TMDs include transverse momenta generated by perturbative QCD evolution, which in turn can be evaluated according to different approximation schemes such as those in [337][338][339][340], [19][20][21], [82,85,86].
If a Monte Carlo method is used to solve the TMD evolution equation, a further advantage is that the solution of the evolution equation can be directly matched to the simulation of parton showers: the kinematic distributions are the same, whether they come from a solution of the evolution equation or from a simulation of the parton shower [87,341,342].
While a general purpose Monte Carlo at TMD level does not yet exist, examples of such algorithms [343][344][345] have been presented for specific cases. We list a few examples below.
• MC event generators with parton shower and hadronization Cascade [341,[346][347][348] is a full hadron level Monte Carlo event generator using TMDs, originally developed for small x processes in ep, now extended to cover medium and large x and pp processes.
Initial state parton showers are treated according to the CCFM formalism, final state parton shower and hadronization is performed by the Lund package Pythia [349]. Parton polarizations are included according to the high-energy factorization [18]. Proton polarizations are not yet included.
Pythia [349]. With the initial and final state parton showers simulated in Pythia, one may argue that several elements of TMD physics are effectively included. Pythia can be used to mimic spindependent cross sections by reshuffling events (assigning polarization states) [343] according to a given cross-section model. This is especially useful when event topologies are needed (e.g., to simulate the interplay of track correlations with detector performance), or where no explicit physics model is yet available to be employed in dedicated MC generators.
mPythia and mLepto are based on Lepto [350] and Pythia [349] with a modification of the hard process [343] to treat the azimuthal angle of the scattered (light) quark and via momentum conservation of the target remnant according to parameterizations of the Sivers function. While limited to the rather specific case of the Sivers effect it can make use of the hadronization embodied in Jetset [351][352][353].
• MC event generators at parton level with fragmentation functions LxJet (see [354]) is devoted to a calculation of jet cross sections at small x in hadron-hadron collisions. It can be also viewed as an event generator as it allows one to generate unweighted events. It uses high-energy factorization [18].
GMC-Trans (see [343]) is a MC generator, developed by the HERMES Collaboration, applying the parton-model expression of the one-hadron semi-inclusive DIS cross section using several models/parametrization for various leading-twist TMD PDFs and FFs. Pion and charged-kaon production is simulated, both for proton and neutron targets (or combinations thereof) without including nuclear effects. An analytic expression for the semi-inclusive DIS cross section was implemented based on the widely used Gaussian ansatz of the transverse-momentum dependences.
TMDGen (see [343]) is an extended version of GMC-Trans entirely written in C++ focusing mainly on di-hadron production in semi-inclusive DIS. Advances in computation power allowed for other than the Gaussian ansatz of the transverse-momentum dependences by employing numeric integration algorithms. It thus allowed the usage of the spectator model [355] for various TMD PDFs and FFs.
Clas (see [343]) uses a similar approach as GMC-Trans, though restricted to the unpolarized sector and to longitudinal double-spin asymmetries. It uses the fully differential single-hadron DIS cross section to simulate semi-inclusive DIS events. The transverse momentum dependence can be Gaussian, but also light-cone quark-model inspired dependence has been implemented.
• Semi-analytical calculations of semi-inclusive processes Resbos [113,114,118] is a package to calculate analytically resummed distributions of inclusive and semi-inclusive observables. The q T resummation in Resbos and parton showering methods of Monte Carlo event generators are complementary. Both are based on all-order resummation using Sudakov form factors. Resbos allows the user to calculate resummed distributions of the Higgs/vector bosons and their decay products up to NNLL. It follows a prescription for matching the resummed contribution onto the fixed-order result and implements a parameterization of non-perturbative effects at small q T in terms of TMD PDFs.
HqT and DYqT [132,356] are numerical programs which implement the analytical q T resummation formalism [29,[131][132][133] to compute, respectively, the q T spectrum of the Standard Model Higgs and Drell-Yan lepton pair (via vector boson production) in hadronic collisions. The resummed results are matched to the fixed order calculation valid at high q T . The program can be used up to NNLL+NLO, with the resummed part evaluated at NNLL, the fixed order evaluated at NLO (Higgs/vector bosons plus one or two partons) and with the normalization fixed to the total NNLO cross section.
HRes and DYRes [357,358] are numerical programs which extend the calculations in HqT/DYqT by retaining the full kinematics of the Higgs/vector bosons and of its decay products. The programs implement q T resummation up to NNLL+NNLO and allow the user to apply arbitrary cuts on final states and to plot the corresponding distributions in form of bin histograms.

VII. CONCLUSIONS
We studied two sets of examples of multi-scale problems in hadronic collisions which require QCD factorization theorems beyond the collinear approximation and call for the use of TMD parton distributions. In one set of examples, the transverse momentum scale is small compared to the hard process scale; in the other, the transverse momentum is of the order of the hard scale but this is much smaller than the total energy of the scattering. In both cases, factorization theorems in terms of TMD parton distributions are necessary in order both to resum logarithmically-enhanced perturbative corrections to all loops and to properly take into account nonperturbative hadron structure effects.
These multi-scale regimes are relevant to LHC phenomenology. An example is the low-q T region of transverse momentum spectra for vector bosons, Higgs bosons, heavy flavor pairs at the LHC. Another example is the production of multi-jets associated with heavy bosons and heavy flavors at large jet masses. Further examples include any final state produced by events at small longitudinal momentum fraction x, such as final states boosted to high rapidities. Besides the LHC, TMD dynamics is central to spin physics in current low-energy experiments and to the planning of future polarized collider and fixed-target experiments.
As the field moves towards the stage of precision studies, appropriate phenomenological tools will be needed. This includes tools for Monte Carlo event simulations, which require parton shower evolution algorithms and determinations of TMD parton distributions from experimental data. First steps toward a new program of portable and accessible TMD pdfs were illustrated with explicit examples in this report.