Different approaches to TMD Evolution with scale

Many apparently contradictory approaches to TMD factorization and its non-perturbative content exist. This talk evaluated the different methods and proposed tools for resolving the contradictions and experimentally adjudicating the results.


Introduction
In the literature there is a bewildering variety of methods for using transverse-momentum-dependent (TMD) parton densities and the associated factorization properties of cross sections. Taken at face value, many of the methods and their uses appear incompatible or contradictory, especially as regards the non-perturbative contributions. The problems are particularly important when planning new experiments to measure polarization-dependent TMD quantities like the Sivers function, since the nonperturbative part of TMD evolution can notably dilute them as energy is increased.
In this talk, I examined and evaluated some of the different methods. I proposed a systematic approach to test treatments of the non-perturbative contributions from large transverse distances (b T ), both from theoretical and phenomenological view points. Then I proposed systematic modifications to the standard parameterizations of the large-b T behavior that could resolve contradictions, especially as regards the apparently incompatible phenomenology of the function controlling evolution of TMD densities. The methods will pinpoint the experimental conditions needed to give incisive experimental probes of the contradictory theoretical statements.

The need for and existence of non-trivial QCD contributions to TMD cross sections
In this article, I use the Drell-Yan process to illustrate issues that apply to TMD factorization in general. For the transverse-momentum distribution in the Drell-Yan process, the simplest model is the parton model, where the TMD cross section is a convolution of the TMD densities for the annihilating quark and antiquark, and the TMD densities do not evolve. In the parton model, the transverse momentum of the Drell-Yan pair directly probes the intrinsic transverse momentum distribution of the quark and antiquark inside their parent hadrons. a e-mail: jcc8@psu.edu

Experimental view
That the parton model description is inadequate in reality (and hence in QCD) is shown by the data in Fig. 1. The graphs also contain several QCD fits to the data. In plot (a) is shown E dσ/d 3 q from the E605 experiment at relatively low Q = 7-18 GeV and √ s = 38.8 GeV. The width is around 1 GeV.
In plot (b) is shown dσ/dq T from the CDF experiment for Z production at √ s = 1800 GeV. This has a much larger width, around 3 GeV. This value is much larger than for the lower energy data, and it also appears incompatible with any reasonable distribution of purely intrinsic transverse momentum. It indicates substantial evolution effects, a specific effect of QCD and other gauge theories.
There is an apparent dramatic difference between the plots at q T = 0. This is merely an artifact of the normalization of the plotted cross section: Plot (b) has an extra factor of q T , which gives a kinematic zero at the origin; for this plot a sensible measure of the width of the distribution is the position of the peak.
The values of parton x are characterized by the ratio Q/ √ s, which is quite different for the two plots. So interpreting the difference between the widths as being associated with evolution with respect to Q is not totally unambiguous; this is recurrent problem. Actual fits [1,2] use other data as well, and appear to unambiguously manifest that there is Q dependence at fixed x.

Need for evolution from QCD
That QCD requires substantial modifications to the parton model is shown on the theoretical side by examining typical graphs that contribute. In Fig. 2(a) is shown the graphical structure of the amplitude for the Drell-Yan process in the parton model. One quark or antiquark out of each of the high-energy incoming hadrons annihilates to make the Drell-Yan pair; the remaining "spectator" parts of the hadrons continue into the final state unchanged, with a big rapidity gap between them. In the parton model, other contributions are assumed to be power suppressed.    However, in QCD there are many other contributions that are not suppressed, as in the example in Fig. 2(b). First, there are final-state interactions that must exist to neutralize the color of the spectator parts. Also, many further contributions exist: Gluons of any rapidity within the kinematic range set by the incoming hadrons can connect any of the other lines, including all of: the active quarks, Fourier trans. of p|ψ WL ψ|p the spectator parts, and the final-state interaction component. Individual graphs do not give a factorized structure. But at leading power in Q, Ward identities and other methods can be used to convert the sum over graphs to a factorized form. The Ward identities are somewhat unusual, and details can be found in [3,Sec. 11.9]. (Earlier literature is lacking fully explicit formulations and proofs.) One consequence is that the parton densities must be defined with Wilson lines, as in Fig. 3. Effectively the Ward identities convert misattached gluons, that link regions of the graph with opposite rapidities, to attachments to Wilson line operators. Further complications involve potential double counting of contributions from different kinematic regions of internal momenta, which must be suitably compensated, and the presence of a soft factor that in recent formulations is absorbed into a redefinition of the TMD densities.
The actual definition of the parton densities is such that the parton densities have extra scale arguments, and must evolve with energy. QCD thereby substantially violates the prediction of the pure parton model that the shape of transverse-momentum distribution scales with energy. The broadening arises because gluons are emitted roughly uniformly into the available range of rapidity, which increases with energy. This applies to both perturbative and nonperturbative gluons.

TMD factorization (modernized Collins-Soper form)
In this section I summarize the formulae of TMD factorization in the form I gave in [3]; detailed proofs were given there. Then then I remark on the location of the nonperturbative information.

TMD factorization
The factorization formula itself for the Drell-Yan cross section is Here, dσ is the hard scattering coefficient, while thẽ f j/H (x, b T ; ζ, µ) are TMD parton densities Fourier transformed into transverse coordinate space. We can set their scale parameters to ζ A = ζ B = Q 2 , µ = Q.

TRANSVERSITY 2014
The evolution equations are (4) Here,K(b T ; µ) is a defined function that controls the evolution of the TMD pdfs and fragmentation functions of light quarks with respect to the ζ parameter.
In the parton model, the integral over all transverse momentum of a TMD parton density is the corresponding integrated, or collinear parton density. Equivalently, when the TMD densities are transformed to transverse coordinate space, the integrated density equals the TMD density at zero transverse separation. In any renormalizable quantum field theory, this result generally needs to be modified. Instead, there is a kind of operator-product expansion (OPE) that expresses the TMD density at small b T in terms of the integrated densities: The coefficients are perturbatively calculable provided that the TMD densities are evolved to scales that avoid large logarithms. The lowest-order value of the coefficients is δ j f δ(x/x − 1), which is the parton model result.

Location of non-perturbative information
The TMD-specific non-perturbative information is at large-b T . Given the existence of the evolution equations, the necessary information is • In the parton densities at large b Tf j/A (x A , b T ; ζ A , µ) at one particular scale. One may choose to label this the "intrinsic transverse momentum" distribution if the scale is low, although this terminology is not entirely accurate.
• In the evolution kernelK(b T ; µ) at large b T . This gives a universal character to the evolution, and can be characterized as giving the effect of "soft glue per unit rapidity".
Predictions for cross sections can only be made with the aid of phenomenological fits for these functions, and/or with the aid of non-perturbative theoretical modeling and calculation. The predictive power of the formalism stems from the universality of these functions: they can be measured from a limited set of data and used to predict cross sections in many other situations, with the aid of evolution and of perturbative calculations of the remaining quantities needed.
The OPE at small-b T also needs the values of the ordinary integrated parton densities. These are obtained from fits to other data than is relevant for TMD factorization. This part of the non-perturbative information is therefore the same as in collinear factorization.

Formalisms used
A list of some of the formalisms that have been used in recent years is: Parton model: Here QCD complications, especially TMD evolution, are ignored.
Non-TMD formalisms: These eschew the use of TMD densities in favor of collinear factorization and a resummation of large logarithms in the massless hard scattering. An old example is by Altarelli et al. [4]; a recent one is by Bozzi et al. [5].
Original CSS: Here a non-light-like axial gauge was used to define TMD densities without Wilson lines, and a soft factor appeared in the TMD factorization formula.
Ji-Ma-Yuan [6]: They implemented the CSS method with gauge-invariant TMD densities with non-light-like Wilson lines. They still had a soft factor, and used another parameter ρ beyond the scale parameters of CSS.
New CSS: Here [3] there is a clean up relative to the original CSS version, Wilson lines are mostly light-like, and (square roots of) the soft factor are absorbed into TMD densities, in such a way that rapidity divergences associated with light-like Wilson lines cancel.
Echevarría-Idilbi-Scimemi [8]: This is a SCET-based formalism, but with a different regulator to handle the divergences given by light-like Wilson lines than is used in the CSS and BN formalisms.
Boer [11], Sun-Yuan [12,13]: These authors start from the CSS formalism, but make certain approximations. Sun and Yuan use no non-perturbative function for TMD evolution.
There is disagreement on size of non-perturbative contribution to evolution, i.e., on the form at large b T of the function that CSS callK(b T ); there is even disagreement as to whether this non-perturbative contribution exists.

Parton Model
The factorization formula (1) reduces to the parton model formula when the hard scattering is replaced by its lowestorder approximation, TMD evolution is ignored, and the high-q T correction term is ignored. The parton-model approximation is typically used to fit data at relatively low energies compared with the earlier Drell-Yan fits. At these energies, a particular interest is in fitting polarizationdependent functions like the Sivers and Collins functions, e.g., [14,15]. Typically a Gaussian ansatz is used for the shape of the TMD functions, e.g., [14].
The OPE (5) for the TMD densities at small b T shows that the Gaussian ansatz cannot be exactly correct and that the Gaussian ansatz will fail once large enough transverse momenta are considered. But it evidently allows a good fit to data at low energy.
The neglect of higher-order terms in the hard scattering is reasonable, since α s (Q) is small. It is also reasonable to neglect the high-q T correction when q T is small enough compared with Q. However, in view of the TMD evolution effects definitely seen at high Q, omitting evolution is not correct when a broad enough range of Q is considered.
However, in reality it is found [16,17] that the data indicates that between the energies of the HERMES and COMPASS experiments, TMD evolution appears to exist but is weak. A complication in coming to this conclusion is that in an experiment at fixed energy, x and Q are highly correlated.

Methods without TMD functions
Some authors, e.g., Altarelli et al. [4] and Bozzi et al. [5], eschew completely the use of TMD densities. They use collinear factorization together with a resummation of large logarithms of Q/q T in higher orders of the massless hard scattering coefficient in the collinear factorization framework. If this were fully justified, it would improve predictive power, since the only non-perturbative information used is in the ordinary integrated parton densities.
However, the justification of collinear factorization uses approximations for large Q that are valid only when q T is of order Q or when q T integrated over. The logical foundation fails when q T Q. The errors in collinear factorization relative to the true cross section are suppressed by powers not only of Λ/Q but also of q T /Q. An important symptom of this is that in the leading power "twist-2" collinear factorization, the effects of Boer-Mulders and Sivers functions are missed, whereas at low transverse momentum these functions given leading power effects. See Ref. [18] for a good description of this last issue.
Further, in the resummation formalism, integrals over scale include non-perturbative regions with, e.g., α s (k 2 ) at small k. A proper TMD factorization shows what to in this region.

Original CSS
In the original CSS formalism [19,20], TMD parton densities were defined in a non-gauge-invariant way with use of non-light-like axial gauge; this was used to cut off the rapidity divergences that would appear if the most natural definition, with light-cone gauge, were used. The CSS evolution formula, of the form of (2), gave the dependence of TMD functions on this rapidity cut off. There was a separate soft function in the factorization formula. Furthermore the evolution equations have power-suppressed corrections, which are dropped in phenomenological applications.
CSS recognized that there are non-perturbative effects at large transverse distance b T . To separate these from perturbatively calculable phenomena, they proposed [21] their b * prescription. The combination of TMD factorization, TMD evolution and the definitions of the TMD densities etc determined what kinds of functions to use for parameterization of non-perturbative parts of the cross section.
On the theoretical side, a difficulty with the use of axial gauge to define parton densities is that the singularities in gluon propagators prevent the direct use of the contour deformations that are used in showing that the effects of the Glauber region cancel in the inclusive Drell-Yan cross section. CSS did not present an explicit solution to this problem. Nevertheless the structure of the formula they presented for the solution of the evolution equations remains as an actually implemented method for comparison with data, and agrees with later results.

Ji-Ma-Yuan
Ji, Ma and Yuan [6] converted the CSS formalism so that the TMD densities were defined gauge-invariantly, with non-light-like Wilson lines. Their factorization formula still has a separate soft factor, like that of CSS. The way in which they derived factorization entail the use of an extra (dimensionless) ρ parameter in the hard scattering etc, with ρ being large. There are associated large logarithms, and the ρ parameter is in addition to the scale parameters of the CSS formalism. There should have been evolution equation for ρ, but such an equation appears not to have been given.
I know of no fits that actually use this scheme. Fits continued to use the CSS method.

New CSS
In [3], I derived an updated, improved version of the CSS results. On the theoretical side: • Covariant gauge was used throughout, with suitable Wilson lines in gauge-invariant definitions of all the TMD functions.
• Full proofs (at least to all orders of perturbation theory) were given, including a proof of cancellation of the effects of the Glauber region that applies both to collinear and to TMD factorization. (This entails particular directions for the Wilson lines.) • A square root of the soft factor was absorbed into each TMD parton density and fragmentation functions (in a rather unexpected, but unique way).

TRANSVERSITY 2014
• As many Wilson lines were made light-like as possible.
The limits are quite non-trivial to formulate, which is a problem that stymied Ji, Ma and Yuan.
• The evolution equations are strictly homogeneous.
The result is substantially cleaner methods relative to the original CSS work. From a phenomenological viewpoint, the new results should be regarded as being at most a scheme change from the original CSS method, as represented by the solution of the evolution equations.

Becher-Neubert
Becher and Neubert [7] obtained a kind of TMD factorization in the framework of soft-collinear effective theory (SCET) in the Beneke-Smirnov style. The results are intended to be valid for large Q with q T Q, but with a restriction to q T Λ (unlike the CSS framework, which does not have this last restriction). By the restriction to q T Λ, they evade issues of a full TMD formalism and the need for non-perturbative information at large b T . But this also means that their method does not apply in the region of low q T , which is of much experimental interest. Thus their methods also do not include the physics associated with Sivers and Boer-Mulders functions, etc, which at leading power show their characteristic effects primarily in the region of non-perturbative q T .
Furthermore they could not define separate TMD pdfs; only the product of two TMD pdfs was defined and free of divergences. This represents an inadequacy of the Beneke-Smirnov approach.
However, the Becher-Neubert method has given an important tool for NNLO calculations of the coefficient functions in the OPE (5) -see [22,23].

Echevarría-Idilbi-Scimemi
Echevarría, Idilbi and Scimemi [8] also obtained TMD factorization in a SCET framework. Their methods are characterized by the use of strictly light-like Wilson lines, but with a different kind of regulator for the associated rapidity divergences. (I do not think it obeys gauge invariance, which causes considerable difficulty in constructing full proofs. Full proofs of factorization make essential use of Ward identities or some equivalent to combine and cancel non-factorizing terms from individual graphs.) As with the method of [3], they absorb soft factors into the definition of TMD parton densities, but in a simpler way that depends on their methodology. Individual TMD parton densities are defined, unlike the case for Becher and Neubert's approach.
In phenomenological fits, Gaussian parameterizations are used for the TMD parton densities at an initial scale. But a claim is made that non-perturbative information is not needed in their equivalent of CSS'sK function that controls the evolution of the shape of TMD functions. Instead, forK, they use a resummation of perturbation theory. This is applied up to a scale of b T = 4 GeV −1 = 0.8 fm or beyond.
In Fig. 4 is shown an example of their results forK, in various approximations.
from which it is clear that b X is closely related to b ΛQCD , up to the µ i -dependent prop Eq. (18), at LL accuracy G(t) = t, and thus A(µ i ) = 1 and b X = b ΛQCD at that accura accuracy for G(t) and considers the available information on the β-function as illustrat finds that 1 ≤ A(µ i ) ≤ 2 for 1 GeV≤ µ i ≤ 1 TeV. We conclude that the divergence of D of the Landau pole, as claimed before. One can calculate the numerical value of Λ QCD , which for n f = 5 and α s (M Z ) = 0 correspondingly b ΛQCD ≈ 7 GeV −1 . At this point we are clearly within the non-pertur accessed by perturbative calculations and has to be modeled and extracted from expe In Section IV and in Appendix B we show how to derive an expression for D R up to Using Eqs. (16) and (B4) we get the asymptotic expression of D R when X ∼ 1 at NN from which one can obtain (approximately) the values of b where convergence is lost. Fig. 1. Thus we can trust D R up to b c ∼ 4 GeV −1 for µ i = √ 2.4 GeV and b c ∼ 6 G that we have used different numbers of active flavors depending on the scale µ i , n f n f = 5 for µ i = 5 GeV, since we have set the threshold of the bottom mass to m b = the larger the initial scale µ i is the broader the interval of the impact parameter w acceptable, and where b ΛQCD is the maximum achievable value. The two cases shown phenomenological cases, between which one should choose the initial scale in order to TMDs.
A last comment worth mentioning concerns the convergence of D R in the small b r convergence of the resummed D is only spoiled in the region around the Landau pol the small b region D R is completely resummable (see Fig. 1) and this agrees with ot series in this region [29].
Summarizing, the resummation method explained above allows us to implement the range of the impact parameter while for larger values of b one clearly needs a non-p discussion of such contribution is beyond the scope of the current work. Then, we can

Geography of evolution of cross section
The evolution of TMD parton densities in formulated multiplicatively in the space of transverse position. In Fig. 5 is plotted the b T -space integrand corresponding to the two cross section plots in Fig. 1. Up to an overall normalization factor, the integrand plotted is b T times the integrand in the TMD factorization formula (1) when µ = Q. To get the cross section, this integrand is to be multiplied by the Bessel function J 0 (q T b T ) and integrated over b T from zero to infinity. In general, in going from low to high Q, the peak region of the integrand migrates to ever-smaller values of b T . We now examine the plots with a black solid line and a purple dot-dashed line. These correspond to fits made to the same data by Konychev and Nadolsky [2] with the same theoretical conditions except that b max = 1.5 GeV −1 = 0.3 fm and b max = 0.5 GeV −1 = 0.1 fm, respectively, for the two lines. At lower energies, in graph (a), the two plots do not differ greatly. At high energy, in graph (b), the two lines match even more closely up to about b T = 0.8 GeV −1 , and then they diverge strikingly, so that the line corresponding to the smaller value of b max is a factor of about two below the other line at the right-hand edge of the graph. Although this is a large difference, it occurs in a region where the integrand is small, so that the large difference has little effect on the actual cross section. The calculation of the cross section is dominated by much smaller values of b T , which are in a perturbative region.
In both cases, the non-perturbative part ofK(b T ) was parameterized by a quadratic function of b T , but the coefficient is substantially larger for the fit with the small value of b max = 0.5 GeV −1 . The plot illustrates a general phenomenon. Although the integral to get the cross section needs an integral over all b T , up to ∞, there is little sensitivity at large Q to the detailed properties of the integrand at large b T , and hence little sensitivity to the nonperturbative dependence at large b T .

Standard fits of TMD evolution give bad low-Q predictions
The standard fits (to Drell-Yan data at Q from 5 GeV to m Z ) use a quadratic form forK,  b T . When the TMD pdfs are evolved backwards, to lower Q, this results in unphysical behavior. To see this, consider the large-b T behavior of the integrand for the cross section, as given in: The x-dependent coefficient is to be obtained from a standard Gaussian fit to data of TMD densities at some initial scale. The coefficient with the ln(Q 2 /Q 2 0 ) factor in the exponent results from applying the CSS equation (2) with a quadratic fit forK(b T ) at large b T .
At low enough Q, the coefficient of b 2 T in the exponent reverses sign, so that the integral diverges at large b T instead of converging. With the BLNY fit, this reversal of sign occurs [13] in a region where there is data and where it is reasonable to apply TMD factorization. This is illustrated in Fig. 6. Even with the KN fit using b max = 1.5 GeV −1 , which gives a smaller coefficient of b 2 T inK, the evolved exponent is well below what is needed to fit HERMES data.
in the nonperturbative form factor, rather than that from the evolution itself. This has to be corrected in order to describe the SIDIS data from the CSS evolution with BLNY and KN parametrizations for the form factors.
On the other hand, for moderate Q 2 variations, we shall be able to understand the Q 2 dependence by directly solving the evolution equation. For example, in the Sudakov resummation formula, Eq. (50), we can, in principle, study the Q 2 dependence by taking the structure functions at lower scale Q 0 as input and calculating the structure function at higher Q by applying the evolution with the direct integral of the kernel from Q 0 to Q. That is the approach we are going to take in comparing SIDIS from HERMES/ COMPASS to Drell-Yan lepton pair production. As was briefly shown in Ref.
[47], this approach works well for the Q 2 range from 2 to 100 GeV 2 and covers SIDIS from HERMES and COMPASS and most of the Drell-Yan processes from the fixed-target experiments. Of course, for extremely high Q such as W=Z boson production at the collider, we have to take into account higher-order corrections and go back to the complete CSS resummation.
In the following, we will show that this evolution approach can describe the transverse momentum distribution in SIDIS and Drell-Yan processes up to Q $ 10 GeV. Since Drell-Yan data can also be understood from the CSS resummation with BLNY (KN) parametrization for the nonperturbative form factors, this provides a natural match between SIDIS and Drell-Yan experiments and helps us understand the TMD evolution in this particular energy range. Once we understand how this works for the unpolarized cross sections, we will extend to the Sivers single-spin asymmetries in these processes.

E. Sun-Yuan approach
In our calculations of the SIDIS from HERMES/ COMPASS, we evolve the cross sections directly from lower to higher scale,W UU ðQ; bÞ ¼ e À W UT ðQ; bÞ ¼ e À F UU ðQ; bÞ ¼ e À F sivers ðQ; bÞ ¼ e ÀS su where the Sudakov form fa The  • This function (or something equivalent) is needed to implement correctly the Q dependence of TMD cross sections.
• Surely b T above about 3 GeV −1 = 0.6 fm is in domain of non-perturbative physics, since we know that the size of the proton is about 1 fm.
• It is difficult to avoid confounding x-dependence with Q-dependence of transverse-momentum distributions.
In measuringK one must be careful to analyze data with different Q at the same value of parton x.
• Fig. 6 strongly suggests that evolution of the shape of TMD parton densities slows down at lower Q compared with what happens in the data fit by BLNY and KN.
• Low Q involves larger (more non-perturbative) b T than high Q.
I propose the following general guidelines for modifying current parameterizations: • One should assume that the KN form (with its b 2 T form) is appropriate only for moderate b T , to fit the higher energy DY data correctly. KN is preferred here over BLNY both because it gives a better fit, and because its value b max = 1.5 GeV −1 = 0.3 fm is not excessively conservative.
• As can be seen from Fig. 5, the data used for the KN and BLNY fits constrainK mostly at b T below about 2 GeV −1 . b * prescription, one has where In Eq. (7),K(b * , µ) is intended to be always perturbative, and all non-perturbative behavior is parameterized in the function g K . To illustrate this, Fig. 7 shows the decomposition ofK with the KN fit. The fitted value of the g K function corrects the cut-off perturbative term,K(b * , µ), and brings the result for the fullK back to its RG-improved perturbative value for b T up to around b T = 2 GeV −1 ; only at higher b T does its curve move away from the diverging pure-PT line. One could therefore argue that the fitting has simply reproduced perturbatively calculable behavior in this extended region, i.e., up to around b T = 2 GeV −1 , perhaps also that the b * method could be improved, and perhaps that b max = 1.5 GeV −1 is still too conservative.

A possible parameterization
One naive idea is that instead of b 2 T , one uses the following parameterization for g K : This goes to a constant as b T → ∞. There are two parameters in (9). Better parameterizations can be found.

Simple ideas for physics constraints on large b T behavior
Given the evolution equation (2), one can characterizẽ K(b T ) as quantifying the effects of the emission of glue for each extra unit of available rapidity, when the energy of an experiment is increased, at fixed x. So, for extra rapidity range ∆y, let • 1 − c∆y = probability of no relevant emission • c∆y = probability of emitting particle(s) • So another possibility for the non-perturbative part ofK isK Here, I have made an ansatz that the transverse-momentum distribution of non-perturbative particle emission at low transverse momentum is Gaussian, motivated by commonly used parameterizations.
We get yet another parameterization, now with quadratic behavior at small b T , and a non-infinite limit when b T → ∞.
Perhaps an exponential at large b T instead of a Gaussian would be better, given known general behavior of correlation functions at large Euclidean distances, as argued by Schweitzer, Strikman and Weiss [24].

Tool to compare different methods:
The A function In a separate talk, I proposed a tool that can conveniently be used to quantitatively compare different methods for TMD factorization in a scheme-independent way. It will be described in much more detail in a forthcoming paper with Ted Rogers.
The motivation arises as follows: • The shape change of transverse momentum distribution comes only from b T -dependence ofK in the CSS formalism, or from some similar quantity.
• Generally in any TMD factorization scheme, the cross section can be written as a Fourier transformation: • So let us define a scheme-independent function 1 where the second line gives its value in the CSS method.
• QCD predicts that this function is: -independent of Q, x A , x B , -independent of light-quark flavor, -RG invariant, -perturbatively calculable at small b T , -non-perturbative at large b T .
It will be useful to compare the values of A(b T ) that correspond to fits and formula in the different articles on the subject of TMD factorization and evolution. The values of parameters where discrepancies occur can be used as a diagnostic: To show which experimental data will be most incisive in arbitrating the correctness of different treatments, and to diagnose which treatments are in disagreement with QCD and whether the disagreements are significant.

Concluding remarks
• Surely we need non-perturbative contributions to TMD factorization. The values of b T that are important in the Gaussian parameterizations of TMD densities are in a region not far from the proton size. Everybody agrees that some parameterization of the nonperturbative properties of TMD densities is needed to describe data at low enough transverse momentum (and hence at large b T ).
• Therefore one must also understand their evolution in this same non-perturbative region of large b T .
• According to established theorems, evolution of TMD functions is governed by a single universal function,K or some equivalent.
• Extrapolation of earlier DY fits to use them at the values of b T relevant for lower energy SIDIS is incorrect.
• It is essential to use better parameterizations ofK so that at large b T its functional form flattens. The parameterizations should be such that they retain compatibility with the evolution measured in Drell-Yan experiments, where substantially smaller values of b T are important compared those needed for the data from the HERMES and COMPASS experiments.
• Physical and phenomenological arguments were given in support of these assertions.
• It is necessary to redo global fits with better parameterizations, and a clear sense of which data are relevant for which regions of transverse position b T .
• In testing and measuring TMD evolution it is essential to ensure that the data being compared are at fixed x with different Q.
• A large coefficient for the b 2 T term inK (and g K ) at large b T causes substantial dilution of the Sivers asymmetry, etc, at large Q, thereby requiring greater sensitivity in future higher-energy experiments. Getting improved understanding and measurements of the non-perturbative part of TMD evolution is important to planning these future experiments.