Predictive power of transverse-momentum-dependent distributions

We investigate the predictive power of transverse-momentum-dependent (TMD) distributions as a function of the light-cone momentum fraction $x$ and the hard scale $Q$ defined by the process. We apply the saddle point approximation to the unpolarized quark and gluon transverse momentum distributions and evaluate the position of the saddle point as a function of the kinematics. We determine quantitatively that the predictive power for an unpolarized transverse momentum distribution is maximal in the large-$Q$ and small-$x$ region. For cross sections the predictive power of the TMD factorization formalism is generally enhanced by considering the convolution of two distributions, and we explicitly consider the case of $Z$ and $H^0$ boson production. In the kinematic regions where the predictive power is not maximal, the distributions are sensitive to the non-perturbative hadron structure. Thus, these regions are critical for investigating hadron tomography in a three-dimensional momentum space.


I. INTRODUCTION
The theoretical study and experimental exploration of the internal structure of nucleons are of fundamental importance to science [1][2][3]. In the past decades, we have obtained a detailed knowledge of the so-called collinear parton distribution functions (PDFs). These collinear PDFs describe the distribution of partons inside a fast moving nucleon as a function of the nucleon's longitudinal momentum fraction x, and thus provide us with a "one-dimensional" (1D) picture of how partons are distributed inside the nucleons. They are indispensable in the predictions involving high-energy hadrons, such as those at the Large Hadron Collider (LHC), in particular for the inclusive observables with one large momentum transfer, e.g., the total cross section of W/Z and H 0 bosons computed in the collinear factorization formalism [4].
On the other hand, for the observables with more than one observed momentum scale, such as the transverse momentum distribution of W/Z and H 0 bosons when the transverse momentum is so much smaller than the mass of the observed particle (q T Q ∼ M W/Z,H 0 ), a more sophisticated factorization framework, namely the transverse-momentum-dependent (TMD) factorization [5][6][7][8][9], is needed. In such a TMD factorization framework, the observables are written in terms of transverse-momentum-dependent PDFs (TMD PDFs), which are usually just called TMDs for simplicity. The TMDs contain not only the aforementioned longitudinal momentum fraction x, but also the partonic transverse momentum k T with respect to the direction of the parent nucleon. Because of this, the TMDs provide us the rich information on "three-dimensional" (3D) motion of the probed active parton inside the nucleon, often referred to as 3D imaging of the nucleon [1,2,10].
Owing to one of the key defining properties of Quantum Chromodynamics (QCD), the color confinement, we do not see any quarks and gluons in isolation. It is therefore critically important to have a reliable and controllable matching between the properties and dynamics of quarks and gluons participating in high energy collisions and the hadrons observed in the detector, which could be achieved by the QCD factorization. Thus, the investigation of the TMDs and the associated TMD factorization becomes extremely important. On one side, they have a strong interplay with high-energy physics, since the uncertainties of hadronic nature as encoded in the TMDs are among the largest ones that dominate the systematic theoretical uncertainties for the QCD calculations of key observables, which could impact our ability to explore the possible scenarios of Beyond the Standard Model physics. On the other hand, a good knowledge of TMDs is essential to map out the nucleon's 3D partonic structure, namely to understand the confined motion of quarks and gluons inside a bound hadron. This is particularly true in light of the rapid progress towards realizing a US-based Electron-Ion Collider (EIC), a machine aiming at investigating the multidimensional structure of hadrons and nuclei. In both frontiers, one of the most important questions to address would be to understand in which kinematic regions the perturbative QCD-based formalism of TMD factorization is most predictive, from the point of view of a controllable perturbative computation. To address this question, it is important to recognize that the probed transverse momentum (k T ) of the active parton in the hard collisions is not the same as the transverse momentum of the same parton inside a bound hadron, sometimes referred to as the intrinsic k T0 and as shown in Fig. 1 for a generic Drell-Yan type hard collision. With the hard collision and the large momentum transfer, a large amount of parton shower is developed during the collision, making the k T of the probed active parton different from the intrinsic transverse momentum k T0 associated with the confined motion inside the bound hadron. The difference between the k T and k T0 depends on the hard scale of the collision, Q, as well as the phase space available for the shower or the total collision energy √ s. The observed Q and √ s determine the momentum fraction x of the active parton participating in the hard collision. The smaller x is, the larger the phase space available for the shower is. The difference between the k T and k T0 is encoded in the QCD evolution of the TMDs in terms of the TMD factorization. As demonstrated quantitatively in this paper, the QCD evolution of the k T -dependence could be dominated by the logarithmic and perturbatively calculable part of the parton shower, leading to a better predictive power. Moreover, the intrinsically non-perturbative TMDs could be further factorized into the non-perturbative 1D collinear PDFs convoluted with calculable contributions from the parton shower. On the other hand, if the evolution of the k Tdependence is dominated by the non-perturbative dynamics of the parton shower, the TMDs and the corresponding observables will be more sensitive to the non-perturbative physics. The detailed and quantitative study presented in this paper will help us figure out in which regions these non-perturbative contributions play a significant role and where the experimental data are most ideal in order to constrain the non-perturbative component of the TMDs.
TMD factorization and evolution have been extensively studied in the literature [5][6][7][8][9][11][12][13][14], together with the matching to collinear factorization [5,[15][16][17][18][19][20], the generalized universality properties [21][22][23][24][25][26][27][28], and the impact on high-energy physics [29][30][31]. Much of the efforts in TMD phenomenology is devoted to the understanding of the role and the size of the non-perturbative corrections in different kinematic domains [13,[32][33][34][35][36]. The study of the kinematic dependence originates from the work of Parisi and Petronzio [37] and Collins, Soper, and Sterman [5], which focused on the value of the hard scale of the process compared to the infrared scale of QCD (Λ QCD ). More recently, it has been shown at the level of the cross section [17,29,38,39] that also the light-cone momentum fraction x, which is effectively a measure of available phase space for the parton shower, could play an important role in determining the relevance of the non-perturbative corrections. In this article we extend those arguments to the context of the modern TMD factorization formalism [9], linking the predictive power of the TMDs to their double scale evolution, i.e. the ultraviolet and rapidity renormalization scales to be defined below. Our detailed study shows that for TMDs with the large hard-scale Q and the small momentum fraction x, the non-perturbative contribution plays a less important role and thus they have the most predictive power. On the contrary, TMDs with the small hard-scale Q and the large momentum fraction x receive significant non-perturbative contributions, and are better suited for constraining non-perturbative parameters in the TMDs.
The paper is organized as follows. In Sec. II we present the structure of a TMD PDF in the coordinate b T space, which is conjugate to the transverse momentum k T . We separate the small and large b T regions, and derive a functional form that extrapolates the physics from the small to the large b T region. In Sec. III we apply the saddlepoint method to the TMD PDF, and we determine the position of the saddle point as a function of the kinematics studying the structure of the double-scale evolution of the distribution. In Sec. IV we analyze the predictive power of the quark and gluon TMDs and comment on the relevance of the large b T region and its components. In Sec. V we study the transverse momentum distributions for Z boson production and H 0 boson production in pp collisions. They are sensitive to quark and gluon TMDs, respectively. We close the paper in Sec. VI and comment on the advantages presented by the complementary kinematic regions accessed by different experiments and possibilities to learn and control the non-perturbative evolution of TMDs.

II. TMDS FROM SMALL TO LARGE bT REGION
Our main focus in this paper is on the unpolarized TMD PDF for a parton with specific flavor a, which carries the collinear momentum fraction x of the parent hadron and has a transverse momentum k T with respect to the hadron's momentum. On the other hand, µ and ζ are the ultraviolet (UV) renormalization and rapidity regularization scales, respectively. As we will discuss in Sec. V, these TMDs are indispensable in describing e.g., the transverse momentum distribution of a vector boson Z and H 0 boson production in the low transverse momentum region q T M Z,H 0 , and carry rich information on the parton's confined motion in a bound hadron, which is a fundamental emergent property of the QCD dynamics. As shown in Fig. 1, the k T -dependence of the TMD PDF probed at the hard collision is a combination of parton's intrinsic k T0 and the amount of k T generated by the parton shower. Since each radiation from the parton shower could be soft and non-perturbative, and convoluted with additional radiation before and after, it could be advantageous to study the TMDs in their Fourier transformed form in the position or b T -space, defined as [5] When b T is small, much less than 1/Λ QCD , the QCD evolution (or scale dependence) of the TMDs' b T -dependence is perturbatively calculable. Otherwise, the QCD evolution is non-perturbative. Once we understand the TMD PDF in the b T -space, we then Fourier transform it back into the momentum space: The zero-th order Bessel function J 0 emerges from the angular part of the integral and the absence of any dependence on the azimuthal angle of the transverse momentum k T in the unpolarized case. The above Fourier transform would require the information of the F a (x, b 2 T ; µ, ζ) for the entire b T ∈ [0, ∞) region. If the Fourier transform is dominated by the information of TMDs at small b T , we will have a good predictive power for F a (x, k 2 T ; µ, ζ) in all relevant k T region, modulo the knowledge of the standard collinear PDFs, as we demonstrate below. On the other hand, if the Fourier transform is sensitive to the large b T region, F a (x, b 2 T ; µ, ζ) will be sensitive to non-perturbative physics since the evolution kernels for the scale-dependence of the TMDs at large b T are non-perturbative.
Below we first review the behavior and evolution of TMDs in the small-b T region, and we then study how one can extrapolate the TMD to the large-b T region by extending the work of Ref. [38]. By further studying the behavior of the TMDs in the b T -space through a saddle-point approximation, we explore quantitatively in which region the TMDs have the most predictive power.

A. TMDs in the small-bT region
The QCD evolution equations of the TMDs take the following form [5] Here the first three equations are well-known and can be found in the literature, see e.g. in Ref. [9], where K (b T µ, α s (µ)) is called the Collins-Soper evolution kernel, and γ F α s (µ), ζ µ 2 is the anomalous dimension of the TMD PDF. The last equation is obtained from the fact that the differential order in ζ and in µ for F a (x, b 2 T ; µ, ζ) is interchangeable [40], i.e., so long as F a (x, b 2 T ; µ, ζ) are differentiable in both µ and ζ in the kinematic regime that we are interested in. In the perturbative region where 1/b T Λ QCD , one can compute all the evolution kernels in the above evolution equations. For example, for a quark TMD PDF with a = q, we have γ F α s (µ), ζ µ 2 = Γ cusp (α s (µ)) ln where we define µ b = c/b T with c = 2e −γ E and γ E = 0.577 the Euler constant. Γ cusp (α s (µ)) and γ (α s (µ)) are the cusp and non-cusp anomalous dimensions, respectively. They generally have the expansion Γ cusp (α s (µ)) = n=1 Γ n−1 αs 4π n , likewise for the non-cusp. For a quark TMD PDF, one has Γ 0 = 4C F , γ 0 = 6C F , etc. At the same time, we have d (1,1) = 2Γ 0 , d (1,0) = 0, etc. Thus at the first non-trivial order, we have The higher-order expressions, and the expressions for gluon TMD PDF can be found in e.g. Ref. [41]. See also Refs. [42][43][44].
Solving the evolution equation, one can obtain the evolved TMD PDF as where µ 0 and ζ 0 are the initial values for the renormalization scales. Integrating Eq. (7) from µ 2 to ζ, one obtains and thus we have Finally when both µ 0 and ζ 0 are in the perturbative region, the TMD PDF F a at the input scales µ 0 and ζ 0 can be re-factorized onto collinear PDFs f b (x, µ 0 ) via an operator product expansion (OPE) at low b T : In practice, one typically chooses the following input values for µ 0 and ζ 0 , to eliminate the logarithms in the coefficient functions C a/b x, b 2 T , µ 0 , ζ 0 . At the same time, one usually chooses µ and ζ to be associated with the hard scale Q, such as the invariant mass of the lepton pair in the Drell-Yan process, Thus in the usual phenomenology we write the perturbative TMD PDF in Eq. (15) in the following form where F a (x, b 2 T ; µ b , µ 2 b ) will be obtained through Eq. (18).
From the discussion of the previous section, we gain the following information. In the kinematic region where Q is large (Q Λ QCD ) and for the small-b T region, one can rely on the perturbative result in Eq. (21) to obtain the information for the TMD PDF F a (x, b 2 T ; Q, Q 2 ). However, for the large-b T region, non-perturbative physics kicks in and the perturbative result is no longer reliable. Several proposals have been introduced to extrapolate the TMD PDF at small-b T into the large-b T region [5]. In this paper, we follow the spirit of Ref. [38] to keep the TMD PDF at small-b T unchanged while we derive a functional form to extrapolate the perturbative result in the small-b T region to the large-b T region. Such an extrapolation would preserve the predictive power of the perturbative calculations in the small-b T region, which is not affected by the extrapolation at large b T , and at the same time, it would provide a physically motivated functional form for the large-b T region. In other words, we write the TMD PDF in the b T -space as where the parameter b max is the largest value of b T at which the perturbative expression for the TMD PDF is trusted (like the input scale at which the DGLAP evolution starts for the 1D PDFs). We choose a rather conservative value for b max = 0.5 GeV −1 throughout this paper. Accordingly, is just the perturbative expression given in Eq. (21) 1 . Here we use the superscript "OPE" to remind that Eq. (21) is connected with the collinear PDFs through an OPE, see Eq. (18). For b T > b max , instead, the non-perturbative correction factor R NP a tames the behavior of F a when the perturbative calculation is not to be trusted. To maintain the continuity of the TMD PDF at b T = b max , the extrapolation function R NP a should satisfy To derive a functional form for R NP a , we take into account the power correction in the evolution kernel [38]. The Collins-Soper kernel K (b T µ, α s (µ)) has an explicit b T -dependence. When b T > b max , we add a power correction into its evolution equation as follows where γ K is an unknown parameter that characterizes the typical size of the high-twist operator. Such a power correction to the evolution equation is also referred to as "dynamical power correction" in [38] and we will continue to use this terminology. For consistency within the TMD evolution equations Eqs. (6) and (7), one would also have With the modified evolution equations, choosing initial scales for the evolution ζ 0 = µ 2 0 = µ 2 bmax and final scales ζ = µ 2 = Q 2 , one would obtain where the input scale µ bmax = c/b max . Setting b T = b max in the above equation, we would obtain 1 In this paper we do not consider the corrections needed for the proper treatment of the region at an extremely small b T [19,32,37,45,46], which is phenomenologically relevant typically at energies lower than the ones considered in our analyses.
where we have used b max µ bmax = c. By comparing Eqs. (26) and (27), we find with the extrapolation function R NP a given by In order to find a reasonable functional form for R NP a , we now need to figure out the following two factors: For the second factor, we turn to the modified evolution equation (24). To proceed, we integrate µ 2 from µ 2 b to µ 2 bmax and obtain To obtain the second line on the right-hand side, we approximate the µ-dependence of γ K (α s (µ)) ≈ γ(µ 2 ) −α with parameters γ and α [38]. We further define the prefactors on the second line to be parameters g 1 and g 2 . Realizing b T µ b = c, we thus obtain Note that the first line on the-right hand side, K (c, α s (µ b )) − K (c, α s (µ bmax )), depends only on b T and b max through the coupling constant α s , and thus such a term can be combined with g 1 -term on the second line (given its connection to the coupling constant), treating g 1 and α as fitting parameters. For the first factor in Eq. (30), we realize that at the input scale µ bmax , one usually mimics the b T -dependence of the TMD PDF F a (x, b 2 T ; µ bmax , µ 2 bmax ) to have a Gaussian form, see e.g. Refs. [47][48][49][50][51], which describes the intrinsic transverse momentum of the partons. With such an approximation, we thus obtain the ratio of TMD PDF at the input scale µ bmax in Eq. (29) as Combining all the above factors, we obtain the following form for the extrapolation function Such a derivation is motivated by the work presented in Ref. [38]. Our derivation is for an individual TMD PDF, while Ref. [38] is for the Drell-Yan differential cross section. This new derivation is based on modern TMD evolution for the TMD PDF, which makes the derivation more transparent and more straightforward. Our derived extrapolation function R NP a automatically satisfies the normalization condition in Eq. (23), i.e. R NP a = 1 at b T = b max . Besides b max = 0.5 GeV −1 we have chosen beforehand, it consists of four parameters α, g 1 , g 2 , and g 2 . While g 2 controls the size of the dynamical power correction,ḡ 2 mimics the intrinsic transverse momentum, which is also referred to as "intrinsic power correction" in [38]. These two parameters are non-perturbative in nature and generally have to be determined from fits to the experimental data. As we will emphasize below, we require F a (x, b 2 T ; Q, Q 2 ) to be smooth at b T = b max , in order to determine the other two parameters g 1 and α in the extrapolation function R NP a . Specifically, we require the first and second order derivatives of F a (x, b 2 T ; Q, Q 2 ) to be continuous at b T = b max . With these two conditions, g 1 and α can be fixed. Accordingly, the R NP a function acquires an implicit x-dependence (equivalent to a √ s-dependence) through these two parameters.

III. SADDLE POINT APPROXIMATION OF A TMD PDF
Once we have the full b T -dependence of a TMD PDF from the extrapolation method discussed in the previous section, we will be able to compute the TMD PDF in the momentum space through the Fourier transformation: where we have set behavior, the integration on the right-hand side, and thus F a (x, k 2 T ; Q, Q 2 ) in the k T -space, will be mainly controlled by the perturbative physics. On the contrary, if F a (x, b 2 T ; Q, Q 2 ) is very sensitive to the large-b T behavior, the non-perturbative physics will play an important role in the behavior of the TMD PDF F a (x, k 2 T ; Q, Q 2 ) in the momentum space. Understanding the TMD PDF F a (x, k 2 T ; Q, Q 2 ) in the momentum k T -space, i.e., whether it is more dominated by perturbative (small-b T ) or non-perturbative (large-b T ) physics, is very important in order to investigate the predictive power of the TMD PDF and of the TMD differential cross sections, which are based on these TMD PDFs. This is the main goal of this and the next sections.
Following Refs. [5,38], we use the saddle-point method to pinpoint if and how the integration on the right-hand side of Eq. (36) is dominated by the small-b T region. The saddle-point approximation, or the method of steepest descent, is often used to approximate the integral when the integrand has the form of e −c S(b T ) , where c is a constant and S a smooth function of b T . As the negative exponential function is rapidly decreasing, one only needs to look at the contribution from where the exponent is at its minimum. Since the TMD PDF in b T -space follows such a form, see Eqs. (21) and (35), it is natural to apply the saddle-point approximation to analyze the TMD PDF. We mainly concentrate on the case where k T = 0. In such a case, J 0 (k T b T ) = 1 and no oscillations are present. When k T > 0, the Bessel function J 0 (k T b T ) further suppresses the large-b T region of the b T integration, and our analysis will be further improved. At k T = 0, we have and thus the integral is dominated by a saddle point at b sp T , which is determined by [5] In the following, we will study in details the kinematic dependence of the saddle point b sp T , in particular the most relevant x and Q dependence: The approximation relates the integral over b T in Eq. (37) to the evaluation of the integrand at the saddle point b sp T . When the saddle point is small, b sp T 1/Λ QCD , i.e., well in the perturbative region, then one would expect the TMD PDF F a (x, k 2 T ; Q, Q 2 ) to be mainly controlled by the perturbative physics (always modulo the collinear PDFs). On the contrary, if b sp T is large, i.e. b sp T 1/Λ QCD , the large-b T non-perturbative contribution is very important and one has to understand/constrain it well, in order to have a full understanding of the TMD PDF. In other words, we use the information on the saddle point b sp T as an indication of the predictive power of the TMD formalism.

A. Saddle point: general behavior
To start, we first use the perturbative contribution to F a (x, b 2 T ; Q, Q 2 ) to compute the saddle-point b sp T . Plugging the perturbative expression in Eq. (21) into (38), we obtain In general, one can evaluate the saddle point b sp T of the TMD PDF by solving numerically the above equation. This is indeed what we do below when we present the results at next-to-next-to-leading logarithmic (NNLL) accuracy and next-to-next-to-leading order (NNLO) in the strong coupling α s . However, at the leading logarithmic (LL) accuracy where one keeps the leading order (LO) result in Γ cusp and in the coefficient functions C a/b , one can solve the above equation and obtain the following simple results where we have introduced µ b = c/b sp T . The function X (x, µ) quantifies the impact of the DGLAP evolution on the position of the saddle point. Its sign changes according to the value of the light-cone fraction x and determines the x-dependence of the saddle point.
The saddle point for the resummed contribution to the Drell-Yan cross section differential with respect to the transverse momentum of the lepton pair has been discussed in Refs. [5,37]. In that treatment the effect of the x-dependence was neglected. In our treatment, neglecting the x-dependence corresponds to setting X = 0. Accordingly, the solution of Eq. (41) reads: where b 0 is the one-loop coefficient of the QCD beta function [52], and n f is the number of active flavors. The expression for b sp (0) T is analogous to the one presented in Refs. [5,37]. It follows the usual wisdom that the larger the value of Q is, the smaller b sp (0) T is, and thus the perturbative contributions to the observable play a more important role. By including the contribution of X , the solution to Eq. (41) acquires an x-dependence: Note that the right-hand side of Eq. (43) depends on b sp T through µ b , and thus Eq. (41) needs to be solved by iterations. A legitimate choice for the first iteration is to evaluate X at b To understand the behavior of X , as well as for the general numerical investigation, below we rely on the LHAPDF6 library [53] and in particular on the central PDF set from NNPDF30 [54] at NNLO accuracy with α s (M Z ) = 0.118. We also use the APFEL library [55] to calculate the X function. The result is in agreement with applying the finite differences method to the NNPDF30 grid. In Fig. 2, we plot X as a function of x for an up quark (left) and a gluon (right), at different scales µ = 1, 5, 10, 100 GeV, respectively. Apart from the gluon case at µ 1 GeV, the function X is positive for x 0.1 and negative for x 0.1. Thus its effect is to reduce the value of the saddle point b sp T with respect to the solution b sp (0) T for x 0.1 and to increase it for x 0.1. Because of this, for the same Q value but smaller x region, the perturbative contribution (from small-b T region) plays a more important role for the TMD PDF. This means that in general, away from the limiting cases, the TMD PDF is more perturbatively dominated at large Q and small x. On the other hand, the TMD PDF is more dominated by the non-perturbative contribution at small Q and large x. This suggests that even for a moderately large Q, the TMD PDF at large x could become quite sensitive to the non-perturbative contribution, due to the x-dependence of the X function.

B. Saddle point: detailed analysis
After the above qualitative understanding of the kinematic dependence of the saddle point, we now turn to a detailed numerical analysis and concentrate on the x and Q dependence. We first choose representative values of x, and study the Q-dependence of the saddle point b sp T . For the small x region, we choose x = 10 −3 which could be relevant to the LHC and the EIC kinematics. While for large x region, we choose x = 0.3 for our illustration below.
Let us first plot the behavior of the x-independent solution b sp (0) T and the x-dependent one, b sp T , both at LL accuracy as given in Eqs. (42) and (43). In Fig. 3  . This is driven by the contribution of a positive X as discussed in the previous section. Similarly, for the large-x region, the purple curves are above the orange curves, i.e., b sp T > b sp (0) T , again consistent with our analysis above. The parameter b max = 0.5 GeV −1 in principle identifies the perturbative region b T < b max , but, considering that this is an arbitrary choice, we can allow for some degree of tolerance and identify the "extended" perturbative region as b T < 1 GeV −1 . In terms of detailed numerical values, we find from Fig. 3 that for the small-x region, the purple curves for up quark is below 1 GeV −1 , i.e. b sp T < 1 GeV −1 when Q 60 GeV, indicating that the perturbative or small-b T contribution plays a more important role for the up quark TMD PDF in the small-x region. On the other hand, for the large-x region, even when Q 120 GeV, the saddle point is still larger than 1 GeV −1 , suggesting that the non-perturbative or large-b T contribution would still play a significant role for the up quark TMD PDF in the large-x region, even though the Q value is already very large. Similar observations apply to the gluon TMD PDF, in an even better way. Due to the larger color factor (C A vs C F ) in Γ cusp , the Sudakov factor makes the gluon TMD PDF more narrowly concentrate in the small-b T region. For example, for a gluon TMD PDF in both the small and large-x regions, the saddle point b sp T would become smaller than 1 GeV −1 for moderate Q 20 GeV already, suggesting that the non-perturbative contribution plays a less important role in determining the gluon TMD PDF. We also note that the LL solution for the gluon at low x becomes non-smooth in the low-to-intermediate Q region: this is essentially due to the non-smooth behavior of the X function.
Let us analyze the saddle point by including the extrapolation term R NP a (x, b T , Q; b max ) in Eqs. (22) and (35). For that, we evaluate the saddle point of the TMD PDF by directly solving numerically Eq. (38) at NNLL and NNLO. In such a setup, we include Γ 0,1,2 and γ 0,1 in the anomalous dimension, and use two-loop results for the coefficient  functions C a/b , as given in Ref. [41]. As previously discussed, there are four parameters in the extrapolation function R NP a (x, b T , Q; b max ), namely α, g 1 , g 2 , andḡ 2 . In Ref. [38], the following quantity is defined and its value at the scale of the W boson mass is determined to be g 2,W = 0.4 GeV 2 through a fit to the experimental data [38]. From a given value of g 2,W andḡ 2 , the value of g 2 can be determined by inverting the above equation. In the analysis below, we either fix the value of g 2 or g 2,W , or vary around their value. As we have mentioned before, we require F a (x, b 2 T ; Q, Q 2 ) to be smooth at b T = b max so to determine the other parameters in the extrapolation function R NP a . Specifically, we require the first and second order derivatives of F a (x, b 2 T ; Q, Q 2 ) to be continuous at b T = b max . With two conditions, two of the parameters can be fixed and we choose to be α and g 1 . There is a subtlety here that requires some caution. In the context of this analysis, which is focused on the high energy regime, we determine α and g 1 through the continuity of the first and second derivative only if the first derivative in b T = b max is negative (∂F a (x, b 2 T = b 2 max ; Q, Q 2 )/∂b T < 0) and thus F a (x, b 2 T ; Q, Q 2 ) decreases as b T increases to be consistent with the expected physical behavior. On the contrary, which is usually the case at very large x, when such a first derivative is positive, we set α and g 1 to zero. This is one of the possible methods to avoid an unphysical extrapolation in the large b T region. Other more flexible strategies that can guarantee non-zero values for α and g 1 can be introduced in order to describe events at low Q, for example in the context of Semi-Inclusive Deep-Inelastic Scattering at fixed-target energies. We leave such a detailed analysis for future studies.
In Fig. 3, we plot the saddle point b sp T for three different scenarios: (1) (g 2,W ,ḡ 2 ) = (0.4, 0.0) GeV 2 , denoted as blue dots, (2) (g 2,W ,ḡ 2 ) = (0.4, 0.2) GeV 2 , denoted as green dots, (3) (g 2,W ,ḡ 2 ) = (0.6, 0.2) GeV 2 , denoted as red dots. It is evident for the small-x and large-Q region that the numerical values of the saddle points are quite stable for both quarks and gluons, such as at Q = M Z = 91.18 GeV (Z boson) and Q = M H 0 = 125.1 GeV (Higgs boson). This suggests that the non-perturbative contributions are mild in these cases. On the other hand, for the quark TMD PDF in the large-x region, the red dots can be different from the blue/green dots even for very large-Q values, suggesting that the non-perturbative contribution could be quite significant. On the other hand, the situation is quite improved for the gluon TMD PDF at large-x, thanks to the strong Sudakov resummation effect. A certain degree of model dependence is left for the gluon at large x and small Q, which anyway vanishes for the gluon at small x, where the saddle point is almost exclusively in the strict perturbative region b T < b max .
In Fig. 4 Fig. 4 is driven both by the x-dependence of the perturbative part and of the non-perturbative part (R N P a ) of the TMD PDF via g 1 and α. Indeed, when b T > b max , if one sets manually g 1 and α to zero, the x dependence is lost. As previously discussed, this is also what happens at (very) large x in all cases apart for the gluon at Q = M H 0 , when the first derivative of the TMD PDF at b T = b max becomes positive. The x dependence generated by the perturbative contribution is generally monotonic increasing. A confirmation of this trend can be found in the shape of the X function. Thus, the changes in concavity in the large x regions are essentially induced by the treatment of the large b T region and thus model dependent.
The overall trend that we can infer from Fig. 4 is that the saddle point for a gluon lies at lower b T values with respect to the quark case at equal or comparable energy scales, again due to the different color factor in the cusp anomalous dimension. It is instructive to point out that, for physical observables which depend on the convolution in momentum space of two TMD PDFs, such as the transverse momentum differential cross section of W/Z and H 0 boson production, the integrand in the b T -space is more peaked in the low b T region than for the single TMD PDF. Thus at large Q and small x region, the predictive power is then guaranteed (see Fig. 8 in Sec. V and Refs. [17,29,38,39]).
In practice, the plots in Fig. 4 suggest that the transverse momentum distribution of H 0 and Z bosons at small-x (or large center-of-mass energy √ s) would be very well controlled by the perturbative contribution. If we are in the small-x region while at the moderate scale of Υ mass, M Υ , the non-perturbative contribution to the gluon TMD PDF could be mild. This suggests that the transverse momentum distribution of the Υ particle could be very well described by the perturbative physics at the collider energy such as the LHC [17,56], where the gluon-gluon fusion channel dominates the production cross section, but not at lower energies. Finally, for the J/ψ production, which is at a very low mass scale M J/ψ GeV, the non-perturbative contribution would be more important and could be even entangled with the formation of the quarkonium [57,58]. For the quark case, the predictive power is well under control at Q = M Z , as we shall see in Sec. V, whereas the physical observables receive significant non-perturbative corrections for Q 10 GeV.
Overall we can conclude that the kinematic domain in which the predictive power is strongest is the large-Q and small-x region, where the saddle point b sp T for the transverse momentum distribution is comparable to or smaller than 0.5 GeV −1 . We emphasize again that in addition to the value of the hard scale Q, this analysis shows that also the value of the light-cone fraction x contributes to determining how relevant the non-perturbative part of the TMD PDF is. This is essential also to understand which experiments and kinematic configurations can be more useful to investigate the properties of the non-perturbative structure of hadrons and which other experimental configurations are more suited for testing the predictive power of the theory. It is certainly important to keep in mind that the predictive power of any theory always depends also on the precision of the specific observable studied in order to test and falsify the theory itself (see Sec. V A).

IV. RELEVANCE OF NON-PERTURBATIVE CORRECTIONS
Apart from the saddle point of the TMD PDF, it is also useful to directly look at the integrand in b T -space of the TMD PDF at k T = 0 , which is simply The shape of this function is also useful to quantify the relevance of the large-b T part of the TMD PDF. In this section, we will assess the relevance of non-perturbative contributions more quantitatively.
In region I, the integrands are completely determined by the perturbative calculation, see also Eq. (22). Note that by construction this region is not affected at all by the details of the model at large b T . The value of b max = 0.5 GeV −1 is marked with a vertical dashed line in Figs. 5 and 6. Region II is a transition region from the perturbative to non-perturbative region. Since we require the TMD PDF to be smooth at b T = b max , the parameters (α, g 1 ) in the extrapolation function R NP a shape the integrand in this region. Finally region III is dominated by the physics beyond the leading power/twist QCD perturbative calculations and non-perturbative, and the values of the parameters (g 2 ,ḡ 2 ) which quantify the strength of the power corrections would mainly determine the behavior of the integrand. Naturally, if the area under region III is very small, the TMD PDF F a (x, k 2 T ; Q, Q 2 ) in the momentum space will be dominated by the perturbative contribution, up to the knowledge of 1D PDFs as indicated in Eq. (18). On the contrary, if such an area is very large, the TMD PDF in the momentum space will be very sensitive to the non-perturbative contributions.  In Fig. 5 and 6, we fixḡ 2 = 0.2 GeV 2 , and vary g 2,W by a factor of 2 up and down from its best fit value 0.4 GeV 2 . As one can see clearly from Fig. 5, for the small-x and large-Q region (x = 10 −3 and Q = M Z ), the non-perturbative contribution from the large b T region b T 1 GeV −1 is moderate. But, at the same time, we find that in this region changing g 2,W by a factor of 2 leads to minor changes in the integrand, as can be seen from the difference in red and blue curves. This suggests that our derived extrapolation function R NP a is mainly determined by g 1 and α and thus be very good in characterizing the non-perturbative contributions in the large-b T region.
For the case of gluons (Fig. 6), the regions I and II dominantly determine the large-b T behavior of the integrand at small x = 10 −3 , for both values of Q = {M H 0 , M Υ }, while the non-perturbative contribution from the large b T 1 GeV −1 to the integrand becomes very small. At large x = 0.3, instead, the power corrections have a mild impact at the H 0 mass scale and a large impact at the Υ mass scale. Once again, this shows that the value of both the hard scale Q and of the collinear momentum fraction x play an important role in determining the relevance of the large-b T input in a TMD PDF.

A. Impact of power corrections
Let us now study the impact of the power corrections: the dynamical power correction as controlled by g 2 and the intrinsic power correction described byḡ 2 , combined in the parameter g 2 (Q) (see Eq. (44)). To quantify the impact of these power corrections on the normalization of F a (x, k T = 0; Q, Q 2 ), we study the following ratio: where the intrinsic power correction is fixed toḡ 2 = 0.2 GeV 2 , and g 2,W = 0.4 GeV 2 . This ratio R pc allows one to focus on the impact of the power corrections only. In Sec. IV B, instead, we will focus on the role of the overall extrapolation term. We consider g 2 (M W ) = 2g 2,W , g 2 (M W ) = g 2,W , and g 2 (M W ) = g 2,W /2, which correspond, respectively, to the blue, the black, and the red curves in Figs. 5 and 6. In Tab. I and Tab. II we present the values of the R pc ratio for x = 10 −3 and x = 0.3, respectively, choosing three different values of Q. Fixing x, the impact of the power corrections is generally larger at lower Q, which means that the TMD PDF is increasingly affected by the non-perturbative corrections at low energies. Viceversa, at fixed Q the impact of the power corrections is more relevant at larger x, which means that in the large-x region TMD distributions are affected by potentially large non-perturbative effects. At small x (Tab. I), by changing Q from M Z to M Υ , the impact of the power corrections on the quark TMD PDF increases by 4 − 5%, whereas at large x (Tab. II) the increase in the quark case ranges from 12% to 65% for the same change in Q. Keeping the value of x and Q fixed, the power corrections are less relevant for the gluon, since its TMD PDF is peaked at a lower value of b T with respect to the quark case (e.g. compare the Q = M Υ cases in Fig. 5 and Fig. 6), due to the Casimir rescaling in the evolution kernel. At small x (Tab. I), the impact of power corrections on the gluon TMD PDF is very low and is not affected at all by changing Q from M H 0 to M Υ , whereas at large x (Tab. II) the impact is comparable to the quark case. Overall, this is a complementary way to prove that TMDs at large Q and small x regions are perturbatively dominated. The choice  k T = 0 is the simplest case since it implies J 0 (0) = 1 in Eq. (3). This eliminates any oscillation from the Bessel function, and allows a better insight into the physics of the small-b T and large-b T regions. When k T > 0, the Bessel function J 0 (k T b T ) further suppresses the large-b T region of the b T integration.

B. Impact of the complete extrapolation term
Let's introduce a cutoff b c for the upper bound of the b T -space integration in Eq. (3): where F a (x, b 2 T ; Q, Q 2 ) is defined in Eq. (22). To test the influence of the large b T -region on the entire TMD PDF, let's introduce the ratio [38]: The ratio R c represents the fraction of the total integral (b c → +∞) generated by the 0 < b T < b c region. Fig. 7 shows the R c (b c , k T = 0) ratio for an up quark at Q = M Z and Q = M Υ , and for a gluon at Q = M H 0 and Q = M Υ . In each panel the ratios computed with x = 10 −3 and x = 0.3 are compared. Let's consider the valueb c such that R c (b c , k T = 0) = 0.75. The latter is highlighted by a horizontal dashed gray line in Fig. 7. For an up quark at Q = M Z ,b c ∼ 1 GeV −1 at low x, whereas at high x one hasb c ∼ 1.5 GeV −1 . Namely, in order to reproduce 75% of the normalization, a wider portion of the large b T region is needed at large x, where it is thus affected by potentially large non-perturbative corrections. The same trend can be observed in the other three cases too. Comparing with Fig. 5 (a), this also confirms that for an up quark at Q = M Z and x = 10 −3 the dominant part of the large-b T correction is the term proportional to g 1 (which is completely determined by imposing the continuity of the first and second derivatives at b max ), whereas at Q = M Z and x = 0.3 also the dynamical and the intrinsic power corrections play a significant role.
Comparing the panels (a) vs (b) and (c) vs (d) in Fig. 7 one sees that, in general, the ratio saturates faster for gluons than for quarks. This is because the gluon TMD PDF is peaked at lower b T values with respect to the quark distributions (see Figs. 5 and 6) due to the stronger suppression in b T space generated by the Collins-Soper kernel K and by the UV anomalous dimension γ F [9,59,60], as already discussed. Looking at Fig. 7 (a) vs (c) and (b) vs (d) one can see that the effect of lowering the value of the hard scale Q is to increase the sensitivity to the power corrections, both for quarks and gluons and both at low x and large x. Comparing Fig. 7 (b) with Fig. 6 (a) and (b) we can see that for a gluon at Q = M H 0 only the term proportional to g 1 is relevant to build the 75% of the total integral. A similar argument holds for Fig. 7 (d), but comparing with Fig. 6 (c) and (d) we can see that lowering Q the distribution becomes increasingly more sensitive also to the power corrections at large x, on top of the g 1 term. In all cases apart from the gluon at Q = M H 0 and x = 10 −3 , comparing with Figs. 5 and 6 one can see that both the g 1 term and the power corrections (namely the overall non-perturbative functions that extrapolates the low b T behavior into the large b T region) become relevant to determine the 95% of the TMD PDF at k T = 0. Eventually, from Fig. 7 (b) we determine that for a gluon with Q = M H 0 and x = 10 −3 the fully perturbative region determined by b T < b max = 0.5 GeV −1 generates the 90% of the TMD PDF at k T = 0, whereas at x = 0.3 it accounts only for 50% of the distribution. This shows that, in principle, the transverse momentum distribution of a Higgs boson produced in gluon-gluon fusion in hadronic collisions can receive non-negligible non-perturbative corrections when one of the two collinear momentum fractions is very large [59], e.g., for the kinematic region far away from the central-rapidity region at the LHC.

V. CROSS SECTIONS
In order to compute the transverse momentum distribution of a Z boson or a Higgs boson produced in hadronic collisions we need to calculate the convolution of two TMD PDFs in momentum space. This corresponds to multiplying the two TMD distribution in the b T -space.

A. Z-boson
For Z production in pp collisions the cross section differential in the transverse momentum q T and in the rapidity of the produced Z in the low q T M Z region reads [32,51,61]: where we have neglected the large q T corrections O(q T /M Z ) to TMD factorization and the corrections O(Λ QCD /M Z ) to collinear factorization. The factors V q and A q are the vector and axial couplings respectively of the Z boson to the quark. The H Z 0 function is the hard function for Z-production: where σ Z 0 is the leading order term [51] and H is the hard function for Drell-Yan with the lowest order normalization H (0) = 1, which we consider at NNLO [62]. We also adopted the narrow-width approximation, i.e., we neglect contributions for Q = M Z . The value of the branching ratio into leptons is B R (Z → + − ) = 0.033658 [63].
As already mentioned, the net effect of multiplying two TMD PDFs in b T space is that the predictive power for the cross section calculation at a specific value of x and Q is increased with respect to the computation of a single TMD distribution, since the product of two TMDs is peaked at a lower b T with respect to a single TMD PDF.
For example, from Fig. 7 (a) we determined that the term ∝ g 1 in the extrapolation function R NP a , as well as the power corrections play a role in determining the value of the quark TMD PDF at Q = M Z , both at low and high x. Instead, in Fig. 8 (a) we show that we can reproduce the data collected by the CMS experiment at the LHC with √ s = 7 TeV and central rapidity −2.1 < y < 2.1 [64] without including any dynamical or intrinsic power correction in the TMD PDF. The g 1 -term in R NP a is sufficient (and necessary) to capture the behavior of the TMD PDF at large b T needed to describe the experimental data. No fit to the data has been performed to reproduce the experimental data in Fig. 8 (a). In Fig. 8 (b), instead, the normalized integrand of the differential cross section in b T space is displayed for q T = 0. For the rapidity values y = ±2.1 and y = 0 the peak of the integrand lies well in the perturbative region. It is also straightforward to check that the peak of ln (b T × the cross section integrand), which corresponds to the analogue of the saddle point for the TMD PDF discussed in Eq. (40), lies at b T < 0.5 GeV −1 . This result is obtained implementing the OPE on the collinear PDFs at small b T at O(α 2 s ) and working at NNLL accuracy with the Collins-Soper kernel K and the UV-anomalous dimension γ F . The x range spanned by the data in Fig. 8 is ∼ [10 −3 , 10 −1 ] (calculated as (Q/ √ s) e ±y ). We checked that for the data collected at more forward rapidity, where one of the momentum fraction x lies in a large x region (e.g. the one . Thus the formalism is predictive and we can describe the data at low qT just with the perturbative contributions plus the g1-term in the extrapolation function, but without intrinsic or dynamical power corrections. We stress that the theory-data comparison in (a) is not the result of a fit. We evaluated numerically the inclusive cross section σ and we find σ = 12.46 nb.
by the LHCb experiment [65][66][67]), the perturbative contribution plus the g 1 -term alone is not sufficient to correctly describe the data, given also their very high precision. This is consistent with our expectation as the relevance of the non-perturbative contribution increases as x gets larger. Indeed it has been recently shown that the large b T part of the TMD PDF is relevant if one wants to describe the very precise LHC data at forward rapidity, and it is also important to take into consideration its kinematic dependence [33][34][35]. Along these lines, we remind that the predictive power is not an absolute concept, but is always related to the precision of the observable under consideration. For example, there might be extremely precise observables for which the perturbative contributions plus the g 1 -term alone are not sufficient to capture the correct behavior at relatively large b T needed to give an accurate description of the quantity considered, even at large Q and small x. This is the case, for example, of the W boson mass, whose determination is sensitive also to the intrinsic transverse momentum dependence and its flavor decomposition [30,31].
Another interesting information available from Fig. 8 is that the TMD cross section given in Eq. (49) (valid in principle at q T M Z , can accurately describe the data in a range of transverse momenta up to ∼ M Z /3, which is comparable to the values quoted in, e.g., Refs. [35,36]. The determination of the range of applicability of the TMD formalism depends both on the perturbative accuracy of the calculation and also on the separation of the b T regions and on the parametrization of the large b T behavior. The determination of the q T range in which the TMD factorization/approximation describes well the data should be, in principle, combined with the error associated to the TMD factorization [20] and can be a useful piece of information in the context of the matching studies [19,20].

B. Higgs boson
In this section we present the calculation for the transverse momentum differential cross section for Higgs boson production from gluon-gluon fusion in pp collisions at √ s = 13 TeV based on the discussed structure for the TMD PDFs. We calculate the cross section in TMD factorization as [60]: where we have convoluted two gluon TMD PDFs in momentum space. The coefficient σ 0 is the Born-level cross section, C t is the coefficient that integrates out the top quark [60], and H is the hard function for Higgs boson production, with the normalization H (0) =1 in the lowest order. For the analytic expression of these coefficients we refer to Ref. [60]. The resummation of large logarithms in the cross section is done by evaluating each perturbative coefficient at its natural scale, and evolving them up to a common scale by using the respective anomalous dimensions [60].
Since there are not enough precise experimental data yet, we compare our formalism to another evaluation of the same observable performed in the framework of collinear factorization with transverse momentum resummation. Specifically, we compare to the resummed result available from the public code 2 HqT [45,68].
Since in this paper we focus only on the unpolarized TMD PDF, we have decided to omit the contribution of the linearly polarized gluons [69] from Eq. (51). Their role in Higgs boson production has been addressed in Ref. [46,60,[70][71][72] and, more recently, in Ref. [73]. Their contribution to the Higgs transverse momentum distribution is known to be of the order of a few percent, depending on the perturbative order and on the implementation of the non-perturbative corrections [46,60,70]. At the phenomenological level the role of the linearly polarized gluons in hadronic collisions is more relevant in the semi-inclusive production of lighter states, such as (pseudo)-scalar quarkonium production at low transverse momentum [57,[74][75][76][77]. In Fig. 9 (a) we compare Eq. (51) at NNLL and NNLO accuracy in TMD factorization with the calculation from HqT at the same perturbative accuracy. The red curve is the calculation based on the formalism presented in this paper assuming an extrapolation to the large b T region without power corrections (g 2 =ḡ 2 = 0). HqT implements the so-called complex-b T [78] prescription to separate the small and the large b T regions. A Gaussian smearing factor in b T space governed by a single parameter g N P is included to account for the potential non-perturbative effects at large b T . The blue band in Fig. 9 (a) has been obtained by setting g N P = 0 GeV 2 and varying the resummation, renormalization, factorization scales by a factor of 2 around the central value M H 0 [45,68]. Changing the parameters controlling the non-perturbative corrections in both approaches has a small impact. In particular, the predictions obtained within our formalism using g 2,W = 0.4, 0.6 andḡ 2 = 0.2 are identical to the red curve in Fig. 9 (a). This is because the support of the b T -space integrand in Fig. 9 (b) is almost entirely in the perturbative region (b T < 0.5 GeV −1 ). The two calculations are in good agreement and compatible within the uncertainty band, and the differences (especially for q T 20 GeV) could be due to the different methods employed to separate the small and the large b T regions.

VI. SUMMARY AND OUTLOOK
In this paper we have discussed the predictive power of unpolarized transverse momentum dependent parton distribution functions (TMD PDFs, or simply, TMDs), as a function of the light-cone momentum fraction x and of the energy scale Q. Such TMD PDFs are essential ingredients in the modern TMD factorization formalism, which generally describes the observables with more than one momentum transfer, such as hadron production in semi-inclusive deep inelastic scattering (SIDIS), and the transverse momentum distribution of vector boson W/Z and H 0 production in hadronic collisions. We have determined that the predictive power is maximal in the large Q and small x kinematic region, for example for vector boson production at hadron colliders with √ s of the order of the TeV and at central rapidity. In other words, the transverse momentum dependence of the TMDs, probed in this region, is dominated by the leading power and perturbatively calculable contributions from the parton shower in the hard collision, and, therefore, the TMDs in this kinematic region, so as the transverse momentum distributions of the bosons, are well predicted by the TMD factorization formalism. Outside of this region, the non-perturbative contributions (as represented by the dynamical and intrinsic power corrections in our study) become increasingly relevant, according to the kinematics explored (non-central rapidity, low Q, large x). Of course this should not be seen as a problem, rather an advantage for probing the nature of hadron structure.
We emphasized that the transverse momentum k T -dependence of parton (quark or gluon) TMDs probed with two-scale observables, Q q T Λ QCD , in high energy scattering is different from the intrinsic k T -dependence of quarks or gluons inside a bound hadron. The difference between the measured k T -dependence of an active parton participating in the hard collision and the parton's intrinsic motion is a result of the QCD evolution of the TMDs. If the evolution is dominated by the perturbatively calculable kernels at small b T , the observed k T -dependence is effectively generated perturbatively, as pointed out in this paper in the region where Q is large and x is small. Such measured k T -dependence of the TMDs is not sensitive to the details of non-perturbative hadron structure other than that included in the 1D PDFs, while its predictiveness is critically important for understanding the production of Higgs particles and other relevant observables. On the other hand, if the measured k T -dependence and its evolution is dominated by the non-perturbative large b T region, which corresponds to the large x and/or not too large Q regime as pointed out in this paper, experimental data of such observables could provide the much needed information for extracting the non-perturbative k T -dependence of the TMDs so long as the TMD factorization formalism is valid. In particular, together with the recent development in extracting the non-perturbative evolution kernels at large-b T from lattice QCD calculations [79,80], we could perform QCD global analysis of such experimental data to extract the intrinsic parton transverse momentum distributions inside a fast moving hadron to shed some lights on the confined motion of quarks and gluons, the fundamental property of hadron structure.
Hadron production at low transverse momentum from SIDIS in the fixed-target mode is probably the configuration where the predictive power from the perturbative contribution alone is the least, and the most sensitive one to the non-perturbative effects [32,49,50,81,82]. It is also the most challenging one from the point of view of factorization theorems [83][84][85][86], given the fairly low value of Q being a couple of GeVs. Vector boson production at RHIC probes a very interesting kinematic region, namely large Q, which guarantees that the factorization approximations are well under control, and relatively large x ∼ 0.1, where the sensitivity to the non-perturbative effects is larger (see Figs. 5 and 7, where we can see a moderate sensitivity to non-perturbative physics for a quark at Q = M W/Z and large x). This could be an optimal kinematic window to study TMD effects, such as the sign change of the Sivers function [25,87,88]. This also naturally applies to the Drell-Yan measurements at COMPASS [89]. Another potentially interesting experimental configuration in the same large-Q/large-x kinematic region is a fixed-target configuration at the LHC [90,91], where several (un)polarized hadron structure measurements could be performed with very high experimental precision, theoretical control on the factorization approximations, and sizable sensitivity to hadron structure effects. Last but not least, also the future US-based Electron-Ion Collider [2] will provide new insights in the quest for hadron structure and hadronization and, in particular, on the TMD PDFs and fragmentation functions (FFs) [92][93][94][95][96]. According to this analysis, a good configuration to probe the quark structure of hadrons in SIDIS at the future EIC could be √ s ∼ 100 GeV for Q ∼ 10 GeV at central rapidity. At the same Q and rapidity and at higher energies, instead, we would be increasingly sensitive to the perturbative structure of the transverse momentum distributions.
This investigation can be expanded in several different directions, for example including small-x resummation effects, polarization effects, studying fragmentation functions, and confronting the given parametrization of R N P a with experimental data from low to high energies. We leave these for future studies.