Transverse momentum dependent distribution functions in the threshold limit

We apply the joint threshold and transverse momentum dependent (TMD) factorization theorem to introduce new threshold-TMD distribution functions, including threshold-TMD parton distribution functions (PDFs) and fragmentation functions (FFs). We apply Soft-Collinear Effective Theory and renormalization group methods to carry out QCD evolution for both threshold-TMD PDFs and FFs. We show the universality of threshold-TMD functions among three standard processes, i.e. the Drell-Yan production in pp collisions, semi-inclusive deep-inelastic scattering and back-to-back two hadron production in e+e− collisions. In the end, we present the numerical predictions for different threshold-TMD functions and also transverse momentum distributions at pp, ep, and e+e− collisions.


Introduction
The parton distribution functions (PDFs) and fragmentation functions (FFs) are two important quantities in particle physics to understand the dynamics of a parton inside a hadron.These functions have been studied extensively both by experiment and theory hadron physics community [1].In the last decade, the hadron physics community proposed a new kind of distribution function called transverse momentum dependent (TMD) distribution functions to extend this studies [2,3].These TMD functions provide information on a parton carrying a certain amount of longitudinal momentum fraction x and transverse momentum k T inside a nucleon, which can be used to probe the quantum correlation between the nucleon spin and active quark or gluon polarization as well as its motion.The measurement of TMD observables provides the leading information on the three-dimensional imaging of a nucleon.
The TMD factorization and resummation framework offer a bridge between TMD functions and observations, which was originally obtained by Collins, Soper and Sterman in [4,5] and has also been derived in Soft-Collinear Effective Theory (SCET) [6][7][8][9][10] based on renormalization group (RG) methods [11][12][13].In the small transverse momentum limit, the differential cross section can be factorized as the product of the hard factor and TMD functions at the leading power.Therefore, one can directly probe TMD functions via different processes, including the Drell-Yan, semi-inclusive deep-inelastic scattering (SIDIS) and back-to-back two hadron production in e + e − collisions.The universality of the nonperturbative parametrization of TMD functions has been investigated extensively in [14][15][16][17][18][19][20][21][22][23][24][25].In the context of the small-x limit, the naive TMD factorization formula at the leadingpower might no longer apply, and numerous studies have addressed TMD functions with considerations for gluon saturation effects [26][27][28][29][30]. Conversely, as x approaches larger values, the underlying factorization theorem of leading-power TMD factorization remains robust; nonetheless, the associated resummation formula begins to incorporate non-negligible logarithmic contributions.Intriguingly, recent lattice QCD simulations [31] demonstrate that in this large x regime, lattice results show inconsistencies with the predictions of various existing non-perturbative parametrization of TMD PDF models [19,20,24,25].Extending our comprehension of TMD functions in this limit is thus of pivotal importance.This study aims to extend our understanding of TMD functions in the large x limit.
In the threshold limit, the phase space for real radiation is restricted, and then the infrared cancellations between real and virtual diagrams leave behind large logarithms ∼ ln(1 − x).Therefore, near the threshold limit, it becomes necessary to take into account these large logarithmic corrections to all orders to have a reliable theoretical prediction.In the Mellin space, these large logarithms are transformed into powers of logarithms of the Mellin variable N .A systematic approach has been proposed to resum these large logarithms to all order [32,33] in the Mellin space and the technique is known as threshold resummation.It is noted that, unlike the TMD resummation, the singular threshold logarithms do not appear explicitly in the physical cross section, since they are always convoluted with PDFs or FFs at the hadron level.After analyzing the dynamical origins of the large corrections in both threshold and TMD resummations, the joint resummation framework of threshold and TMD logarithms was first developed in [34].Such a framework has been applied in various processes [35][36][37][38][39] at hadron colliders.Later on, a factorization and resummation formalism based on SCET + [40,41] was introduced in [42,43], which can be used to perform resummation calculation beyond the next-to-leading logarithmic (NLL) accuracy.
In this paper, we introduce a new type of unpolarized TMD functions, threshold-TMD distribution functions, within the joint threshold and TMD factorization and resummation framework.We apply the crossed threshold resummation method [44] to find a close correspondence between resummation for Drell-Yan, SIDIS and e + e − processes.Therefore, one can obtain the resummation formula for each of the processes using the same procedure.In the joint limit, the cross section is factorized as the product of hard factor and threshold-TMD functions, including threshold-TMD PDFs and FFs, and the structure of the logarithms turns out to be identical.Among these three processes, we have the universality among these unpolarized threshold TMDs (TTMDs).Explicitly, we find This property also appears in the standard TMD factorization and resummation formula, which is very useful in the global fitting of different sets of TMD functions.In principle, such property of universality can also be generalized to the polarized TMD functions, which will be discussed thoroughly in future work.In this paper, we present the numerical results on the unpolarized threshold-TMD functions and also the transverse momentum cross section in three processes.In numerics, we restrict ourselves to resummation at NLL, which captures the main effects in the QCD evolution.The rest of this paper is organized as follows.In section 2, we first review the factorization theorem in the joint threshold and small transverse momentum limit for the Drell-Yan process.Then, we introduce the definition of threshold-TMD PDFs and write down the corresponding QCD evolution function.In this section, we also briefly show the factorization theorem for SIDIS and e + e − processes and give the definition of threshold-TMD FFs.We present the numerical results in section 3, where we first give the numerical predictions of the transverse momentum distribution for threshold-TMD PDFs and FFs, and then discuss the cross section for Drell-Yan, SIDIS and also e + e − .We conclude in section 4. The details of perturbative results of the collinear-soft function are provided in the appendix A. We also collect the fitting parameters of PDFs and FFs used in our numerics in the appendix B.

Factorization and resummation formalism
This section presents the factorization theorems for processes, characterized by the leading order partonic reaction q q → γ * or one of its crossed versions, in the joint threshold and small transverse momentum limits.First, we show the factorization theorem for the Drell-Yan process and present our definition of threshold-TMD PDFs based on RG equations.Then, we generalize our results to SIDIS and e + e − processes and introduce the definition of threshold-TMD FFs, which captures both soft gluons and TMD evolution effects for the inclusive hadron production.

Theory formalism in Drell-Yan
For simplicity, we consider the Drell-Yan process mediated by a virtual photon with timelike momentum q µ , where X denotes the undetected hadronic particles in the final state.The standard TMD factorization theorem reads where q T , Q and Y denote the transverse momentum, invariant mass and rapidity of finalstate lepton pairs respectively.The differential cross section is factorized in terms of hard factor H(Q, µ) and TMD PDFs f TMD q/h i (x i , b T , µ, ζ) from two colliding beams where µ and ζ denote the factorization and Collins-Soper scale respectively, and the soft function has been subtracted as in the standard redefinition procedure [45,46].Here q runs over quarks and anti-quarks participating in the hard scattering, and e q denotes the corresponding charges.The light-cone momentum fractions of the hadron h 1,2 carried by the partons are denoted by x 1,2 , which can be expressed as In addition, for simplicity, we define the Born cross section as with the number of color N c = 3 and the fine structure constant α em , so that the leadingorder hard function H DY (Q, µ) is normalized to unity.The perturbative calculation of TMD PDFs using the operator-product expansion method shows that large logarithmic terms will appear in the limit x → 1 [47][48][49][50][51][52][53][54][55][56][57][58].Therefore, one needs to consider the factorization of the cross section in the joint threshold and TMD limit.Following [44] and performing the Mellin transformation with respect to the threshold variable τ DY , we obtain where from the first line to the second line the Jacobian for converting between (Q 2 , Y ) and (x 1 , x 2 ) can be easily worked out using the definitions of the kinematic variables in (2.3).Therefore, the differential cross section can be re-expressed as where we have applied the inverse Mellin transformation C N • • • to obtain the cross section in the momentum space.In the moment space, TMD PDF is defined as It is worth emphasizing that such transformation is not necessary for the standard TMD factorization theorem, but we keep it for the following joint threshold and TMD factorization.Now we briefly discuss the factorization formula in the joint threshold and TMD limit in SCET formalism.More detailed discussions on the NLL resummation formula can be found in [35,36], and the factorization formula within SCET is given in [42,43].The threshold effects embody the behavior of the cross sections at τ → 1, since near the machine threshold s ≈ Q 2 , the colliding energy is just sufficient enough to produce the final-state lepton pair with the invariant mass Q.Therefore, one immediately realizes that the threshold limit in momentum space τ → 1 corresponds to N → ∞ in Mellin space.This means the above factorization is incomplete in the sense that the threshold kinematics are ignored.Generally, there is a standard operator expansion step to relate TMD PDFs to collinear PDFs at small b T [5,45,59].Explictly, we have where C ij (z, b T , µ, ζ) are perturbative matching coefficients.It is known that to arbitrary order there are threshold logarithms [ln n (1 − z)/(1 − z)] + in C ij as z → 1, which must be resummed to achieve precise results.To incorporate both the TMD and threshold effects, we employ the re-factorization technique presented in [42,43] for the TMD PDF f TMD i/h where we introduce the Collins-Soper scale ζ N ≡ ζ/ N 2 with N = N e γ E and the Euler constant γ E in the threshold limit.The rapidity scale dependence ν is cancelled between the Mellin space unsubtracted collinear-soft function Sunsub ) and the standard TMD soft function S(b T , µ, ν).We refer to Sc (b T , µ, ζ N ) as the genuine collinear-soft function in the joint threshold and TMD limit, which is flavor and spin-independent, but different for quarks and gluons.In the appendix A, we provide the derivation of its oneloop expression in the perturbative region and also discuss its relation to the perturbative matching coefficient of TMD PDFs in threshold limit.Moreover, we derive its expressions up to the three-loop order based on the threshold behaviors of the next-to-next-to-nextto-leading order (NNNLO) perturbative matching coefficients of TMD PDFs [56][57][58].In Eq. (2.9), fi/h (N, µ) is the collinear PDF in Mellin space, which takes the form fi/h (N, µ) (2.10) It is important to emphasize that the aforementioned refactorization formula (2.9) is only proven for small values of b T [42,43].The generalization of this formula to higher-twist distributions remains unclear, and further investigation is required to address this issue.We leave the study of higher-twist distributions for future works.Since both the unsubtracted collinear-soft and TMD soft function depend on the rapidity scale ν, their rapidity scale evolution is governed by the following ν-RG evolution equations [12,60] with γ Sc ν = −γ S ν /2.Upon subtracting the soft function, as previously outlined in Eq. (2.9), we have the genuine collinear-soft function Sc (b T , µ, ζ N ).In particular, we find that Sunsub c (b T , µ, ζ N /ν2 ) depends only on the combination ζ N /ν 2 in the threshold limit1 , which is verified up to the three-loop order in appendix A. This allows us to convert the above two ν-RG evolution equations in Eqs.(2.11) and (2.12) into a similar Collins-Soper evolution equation for genuine collinear-soft function Sc (b T , µ, ζ N ) by following the procedure outlined in [46] where [45,46] or the rapidity anomalous dimension [12,60] and collinear anomaly exponent [11,62] in SCET.Here, from perturbative calculation, κ(b T , µ) can be written as We have verified Eq. (2.13) up to three-loop order based on the NNNLO expression of Sc given in Eq. (A.13).The four-loop perturbative expression of κ(b T , µ) was recently calculated in [63,64], and its nonperturbative analysis can be found in [65,66] and the preliminary nonperturbative numerical simulation in lattice QCD is presented in [67][68][69].
It is noteworthy that the Collins-Soper equation governing Sc as a function of √ ζ N coincides with the equation for TMD PDFs as a function of √ ζ.This correspondence arises from the (ν-)RG invariance analysis elucidated in Ref. [42].Specifically, the rapidity divergence remains the same in the threshold limit, defined by Q ≫ Q(1 − τ ) ≫ q T , where τ is the partonic analog of τ DY .Consequently, the rapidity-RG equation (2.11) for the unsubtracted collinear-soft function is unchanged.It is crucial, however, to highlight that the rapidity scale differs between the unsubtracted collinear-soft function and the unsubtracted TMD PDFs.In essence, all threshold logarithms present in the matching coefficients of TMD PDFs are subsumed into the rapidity logarithms.This has been rigorously confirmed up to three-loop calculations [56][57][58].Therefore, we introduce a modified Collins-Soper scale ζ N in the joint threshold and TMD limit.
After solving the above Collins-Soper evolution equation for the ζ N dependence at the renormalization scale µ = µ b , we have where we can choose T and ζ N,f will be determined from RG consistency [45,46].To obtain its value, we first write down RG equations for the hard function, collinear-soft function and threshold PDFs ) and the corresponding anomalous dimensions are with perturbative coefficients of anomalous dimensions needed at the NLL accuracy 2 as where they are defined by Then the RG consistency Γ h + 2 Γ Sc + 2 Γ fq = 0 implies that (2.24) We observe that the value of ζ N,f is different from the one in the standard TMD factorization formula.In the momentum space the scale ζ N,f can be expressed as We note that its value is reduced from the original Collins-Soper scale ζ f = Q 2 in the standard TMD factorization.This is expected since we have taken into account the threshold effects in the joint factorization formula, and the phase space for the initial collinear radiations is further constrained in the threshold limit.Especially, as , the rapidity evolution effects in (2.14) will be turned off automatically.
All-order resummation formula can be obtained by solving the RG equations in both position and moment spaces and evolving the ingredients from their intrinsic scales to a common scale.The all-order resumed cross section is given as 2 Here, we provide the two-loop cusp anomalous dimensions that are employed in Γ h and Γ fi .However, it is important to note that, for the Collins-Soper kernel κ, we consistently utilize its one-loop results in our NLL resummation calculations.
where we introduce a new type of TMD PDFs, i.e. the threshold-TMD PDFs f TTMD i/h , and then the cross section could be factorized as the product of the hard function and threshold-TMD PDFs.The definition of the threshold-TMD PDF at a scale of Q is where the perturbative evolution kernel S pert involves the contribution from both collinear PDF and collinear-soft functions, which reads Here note that the Γ h term is the same as the standard TMD PDFs since the hard function only receives virtual pQCD corrections and is unchanged in the threshold limit.The Γ fi term may be viewed as a simplified DGLAP equation, where the flavor off-diagonal pieces are non-singular as N → ∞, and thus can be neglected at the leading power formalism.Besides, µ b * and µ F are the intrinsic scales of collinear-soft function and collinear PDFs.We stress that the above perturbative evolution kernel is consistent with that in NLL resummation formula [35,36], and it can also be evaluated beyond the NLL accuracy [42] after including higher order ingredients within the perturbative QCD framework.
To match the perturbative and non-perturbative contributions in (2.26), we first apply the standard b * -prescription [5]  (2.28) In addition to the b * -prescription, we also apply the model in [14,23] to parametrize the non-perturbative contribution at large b T .Explicitly, in (2.26), we define the nonperturbative kernel as We observe that the value of ζ N,f in the threshold limit deviates from its counterpart in standard TMD resummation, as elucidated subsequent to Eq. (2.24).Accordingly, in Eq. (2.29), we opt for the Collins-Soper scale factor Q/ N rather than Q as originally posited in Ref. [14].We must underscore the fact that our adaptation of the model does not account for the complete non-perturbative threshold enhancement effects.Therefore, additional fits in the threshold region are requisite for achieving high-precision results.Lastly, apart from the non-perturbative factor S NP associated with TMDs, it is necessary to incorporate the collinear PDF, encapsulated in the factor fi/h (N, µ F ), defined at the scale µ F .

Theory formalism in SIDIS and e + e −
In addition to TMD PDFs, another important set of distributions for probing hadronic three-dimensional structures is the TMD FFs, which can be studied in SIDIS and back-toback two hadron production in e + e − collisions, separately.Similar to the structure of the previous subsection, we first review the factorization in terms of usual threshold variables in SIDIS and e + e − .Then we define the TMD FFs in the threshold limit.As we will see, our formalism has an advantage over the usual pure threshold formalism allowing further discussions.Finally, we present the factorization for SIDIS and e + e − in the joint threshold and TMD limit and a collected version.
For the SIDIS process, we consider a proton p with momentum P µ 1 is probed by a virtual photon with space-like momentum q µ and produces a final-state inclusive hadron h with momentum P µ 2 .Explicitly, we have which probes the short-distance scattering of the electron and a quark inside the proton p by exchanging a virtual photon.The standard unpolarized differential cross section is [23,45,70] where q T is the transverse momentum of the photon in the hardon proton frame, and σ DIS 0 is the leading order electromagnetic scattering cross section given by Besides, H SIDIS denotes the hard function in the SIDIS process, and D TMD h/q is the standard TMD FF.Note that we have included a factor of z 2 into the definition of D TMD h/q (z, b T , µ, ζ).Here, we have employed the usual SIDIS kinematic variables with q = ℓ − ℓ ′ , Q 2 = −q 2 and s = (P 1 + ℓ) 2 .In the Mellin space, the above factorization takes the form where, following [44], we have τ SIDIS ≡ xz that is very similar to the Drell-Yan threshold variable τ DY .For this reason, this is referred to as "crossed threshold variable" in [44].
After the inverse Mellin transformation, the above factorization can be rewritten as In the threshold limit i.e τ SIDIS → 1 in momentum space or N → ∞ in Mellin space, the above factorization is incomplete.One needs to re-factorize the TMD FFs to take into account this threshold effect.We find that the re-factorized formula can be written in Mellin space as (2.36) where Dh/q (N, µ) is the collinear FF in Mellin space.It is not surprising to find that the genuine collinear-soft functions for TMD FFs are the same as the ones in TMD PDFs.Since, in the threshold limit, the longitudinal momentum fractions of TMD FFs z → 1, making the behavior of the TMD FFs similar to that of the TMD PDFs at NLL.Moreover, we stress that it is one of the advantages of our joint threshold and TMD formalism to allow one to define threshold distributions separately.
To obtain the Collins-Soper scale, one exploits the RG consistency by and the corresponding anomalous dimensions are with perturbative coefficients of anomalous dimensions needed at the NLL accuracy as

46)
Then Γ h + (Γ Sc + Γ Dq ) + (Γ Sc + Γ fq ) = 0 just implies exactly the same Collins-Soper scale as that in the Drell-Yan process (2.47) Therefore, we define the threshold-TMD FFs in the Mellin moment space at a scale Q as with the perturbative evolution kernel S pert is defined as The nonperturbative factor S D NP is defined as with g D 1 = 0.042 GeV 2 [14,23].The values of Q 0 and g 2 have been given in (2.29).It is noted that the above parametrization is consistent with the model used in [14,23] in the threshold limit z → 1.Also, as mentioned earlier, a new fitting including the threshold effect is required for a higher precision theoretical result.Thus, the all-order resummation formula reads Finally, we move on to the factorization and resummation for the inclusive back-to-back two hadron production in e + e − collision, i.e.
where P µ i is the momentum of the hadron h i , and the threshold variable in this process is defined by with Q 2 = q 2 .As expected, in the joint limit, the cross section is factorized as the product of hard factor and two TMD FFs that describe final-state di-hadron production.Besides, the RG consistence also implies the Collins-Soper scale is the same as that in the Drell-Yan and SIDIS processes.Therefore, we obtain the following all-order resummation formula.
where σ e + e − 0 = 4πα 2 em /Q 2 is the Born cross section and H e + e − is the corresponding hard factor, and the definition of the threshold-TMD FFs at the scale Q is given in (2.48).
From now on, we have obtained the all-order resummation formula for Drell-Yan, SIDIS and e + e − processes, and also find a close correspondence between them.In the joint limit, the cross section is factorized as the product of the hard factor and threshold-TMD functions, including TMD PDFs and FFs.Such a property of universality is significant in the future global fitting analysis for threshold-TMD functions in different processes.To summarize this section, we give the following generic resummation structure for all three processes as where for A = DY; FTTMD , and for A = e + e − ; FTTMD and the Born cross sections σ A 0 and hard function H A are defined accordingly.Here in addition to τ , the Born cross section σ 0 can also depend on additional variables, for example, σ SIDIS 0 depends on the inelasticity y as shown in (2.51).This structure holds in the all-order QCD resummation formula.

Numerical results
In this section, we numerically study the threshold effect in TMD PDF and TMD FF using the factorization formula of Drell-Yan, SIDIS, and the back-to-back di-hadron production via the e + e − annihilation process.We apply these extracted TMD functions to provide transverse momentum distributions for these three processes for different experimental kinematics.Throughout our numerical analysis, we use one loop strong coupling constant (α s ).The fine structure constant (α em ) is taken to be 1/137 and the number of active quark flavor n f = 5 in the massless limit.For the TMD and collinear evolution, we choose our initial scale (µ F ) to be 1.3 GeV.As mentioned earlier, threshold factorization has been done in the Mellin space because, in this space, the convolutions become a simple product.After achieving the factorization formula in Mellin space, we need to do the Mellin inversion to achieve results in the x space.During our numerical study, we observe that one needs to modify the Collins-Soper scale in order to achieve results that are consistence throughout the allowed kinematic region of τ .This modified Collins-Soper scale will introduce new poles in the Mellin inversion formula.We discuss these issues in detail in the next subsection and present a prescription for Mellin inversion.

Modified Collins-Soper scale and inverse Mellin transformation
In the NLL resummation formula, the Collins-Soper scale in the non-perturbative factor S NP (2.50) is given by Q/ N .When the moment variable N becomes very large, ζ TTMD f goes down to the nonperturbative scale Q 2 0 , which violates the power counting, Q ≫ Q(1 − τ ) ≫ q T , in the factorization theorem.Therefore, we need to freeze the Collins-Soper scale as which reproduces the resummation formula in the moment variable N as long as Q/ N ≫ Q 0 .Besides, it also allows us to perform a straightforward numerical calculation in the complex-N plane.We leave the investigation of different function forms and corresponding theoretical uncertainties for future studies.
To implement the inverse Mellin transformation, one often rewrites the standard textbook Mellin inversion [71] as an integration over a real variable with adjustable contour: Then contour is characterized by a real number c, which should be the right to the rightmost singularity of φ n , and an arbitrary angle ϕ as shown in Figure 1.The parametrization forms we chose are certain simple combinations of power functions, so they only have poles on the real axis.In this regard, one can always find a suitable c and may try various ϕ to reach higher numerical efficiency.There are different Mellin inversion prescriptions available in the literature [71,72].However, when we modify the Collins-Soper scale term from ζ f to ζ * , we find that this modified term would contribute two new poles with Q ≡ e γ E Q.These two singularities would cause numerical issues even when the contour might not exactly pass the poles for different Q/Q 0 .For simplicity, we chose ϕ = π/2 for our contour, which is safe no matter what the Q value is.Also, we choose c = 1.6 to avoid all the other poles.Numerically, we vary the value of ϕ and c and find the results are stable and do not depend, within errors, on the particular choice of input parameters.The collinear PDF sets used in this analysis are parameterized according to the CT18NNLO prescription, with the parameterized equations specified at the initial scale µ F = 1.3 GeV.

Threshold-TMD PDFs and FFs
In this subsection, we present the numerical results for threshold-improved TMD PDFs and FFs.To investigate the impact of threshold resummation, we compare the TMD PDFs with and without threshold resummation in Figure 2, which shows the transverse momentum distributions of TMD PDFs for the up (u) quark in the proton.We vary the value of x from 0.1 to 0.8 at 5 GeV.The green curves represent the results obtained without threshold resummation, corresponding to the case where ζ N,f = Q 2 in (2.26).The red and blue curves correspond to the results obtained with threshold resummation using the ζ * scheme (3.1) and naive replacement ζ N,f = Q 2 / N 2 , respectively.In the low x region, all three curves exhibit consistent behavior since the threshold effect is mild.However, as we increase the value of x, the threshold effects become evident.Especially, it is noteworthy that in the naive scheme the theoretical predictions for the distribution rapidly decrease to zero and then turn negative.Consequently, we conclude that the ζ * prescription is necessary to ensure reliable theoretical predictions.
In Figure 3, we present the threshold improved TMD PDF for u and d quarks.We make the use of CT18NNLO [73] parameterized colinear PDF set at a given initial scale µ F = 1.3 GeV.The upper panel is for the u quark and the bottom panel is for the d quark.We present our results for three different choices of hard scales; Q = 2.0, 5.0, 10.0 GeV with three different choices of Bjorken scales x = 0.4, 0.6, 0.8, from left to the right respectively.In the threshold region when the hadronic threshold variable is in the order of 1, the Bjorken scales are also in the order of 1.We present our results for a large Bjorken scale to highlight the threshold effect in this region.The uncertainty bands correspond to the 1-σ variation from CT18NNLO PDFs using the Hessian method.The details of these uncertainties are given in appendix B.
In Figure 4, we present threshold improved TMD fragmentation function of Pion (π + hadron).In this case, we use the parameterized JAM20 fragmentation function [74] at the scale µ F = 1.3 GeV.Similar to Figure 3, the upper panel is for the up-type and bottom panel is for the down-type quark fragmentation functions and z is the fragmentation variable.Here, we present our results for three different values of hard scale (Q) with three different values of fragmentation (z) variables.The uncertainty bands correspond to 1-σ variation of JAM FFs using the replica method.
At this point, we would like to emphasize that the direct comparison between the usual TMD functions and threshold-improved TMD functions is not trivial.In our theoretical formalism, we evolve the TMD evolution factor from the initial scale µ F to the desired final scale.On the other hand, in the standard TMD evolution, one starts from scale µ * b .Besides, both perturbative and nonperturbative parts in the Collins-Soper evolution are also different in these two frameworks.Because of this different evolution scheme, it is not straightforward to compare them and we leave such a study for future investigation.

Predictions for different experiments
Finally, we use previously determined threshold-improved TMD functions (PDFs and FFs) to predict some experimental results.For this, we investigate three different processes namely, DY, SIDIS, and e + e − annihilation process.For the DY process, we consider dilepton transverse momentum distribution from proton-proton collision.We choose the hadronic center of mass energy be √ s = 15.0GeV and Q = 6.0 GeV, which is relevant to the Drell-Yan production at typical Fermilab experiments.The results are presented in the left panel of Figure 5.For the e + e − annihilation process, we consider Pion (π + ) in the final hadronic state.We choose the virtuality of the photon to be Q = 10.58GeV and the threshold variable τ = 0.64 for the BELLE [75] kinematics.The results are presented in the right panel of Figure 5.The uncertainties in all the cross sections are coming from the threshold-TMD functions which are computed from the 1-σ variation of PDF and FF using Hessian and Replica methods respectively.As shown in the previous subsection, the uncertainty from FFs is larger than that from PDFs at a similar value of the light-cone momentum fraction making a wider uncertainty band in the e + e − process than in the Drell-Yan process.
Finally, in Figure 6, we present the results for the SIDIS process.In this case, we make use of four experiments kinematics namely HERMES [76], COMPASS [77], the future Electron-Ion Collider (EIC) [3] at Brookhaven National Laboratory, and the 12 GeV program [78] currently underway at Jefferson Lab (JLab).The upper left panel is for HERMES kinematics of ep collision with hard scale Q 2 = 4 GeV 2 and threshold parameter τ = 0.24 and the upper right panel is for COMPASS kinematics with hard scale Q 2 = 20 GeV 2 and threshold parameter τ = 0.32.In the lower panel, we present the results for EIC (left, Q 2 = 50 GeV 2 and τ = 0.4) and JLab 12 GeV program (right, Q 2 = 3 GeV 2 and τ = 0.48) kinematics.To conclude this section, we notice that although the future EIC will make the most precise SIDIS measurements at small x (down to x ∼ 10 −4 ), the EIC will also increase the precision of the data in the large-x region up to x ∼ 0.5.The JLab 12 GeV program will make precision measurements up to x ∼ 0.6 and smaller values of Q 2 .We expect both EIC and JLab 12 to make important constraints on these threshold TMD PDFs and TMD FFs in the near future.

Conclusion
In this paper, we have introduced the threshold TMD functions namely, PDFs and FFs.To probe this threshold effect, we use three different processes: Drell-Yan, SIDIS and e + e − processes.Considering the extra degrees of freedom in this kinematical region known as collinear-soft degrees of freedom, we re-factorized the formula for these processes and defined threshold-improved TMD distribution functions.We observe that because of the phase space reduction in the threshold region, the Collins-Soper scale is not the same for the usual TMD and threshold-improved TMD distribution functions.We find that ζ N,f = Q 2 / N 2 for threshold TMD whereas ζ f = Q 2 for the usual TMD.During our numerical analysis, we observe that apart from the RG consistency, the Collins-Soper scale has to be chosen in such a way that it is always greater than the non-perturbative scale Q 0 .Therefore, one needs to modify the Collins-Soper scale, and this modified scale introduces a new kind of pole in the integration contour.We provide a Mellin inversion prescription to avoid all kinds of poles in the integration contour.Finally, we provide the numerical predictions for transverse momentum distributions of Drell-Yan, SIDIS and e + e − processes for different experimental kinematics using these threshold-improved TMD functions.
Our theoretical formalism would open a new window into TMD physics.Future experimental analysis and global fitting analysis will certainly help in understanding the non-perturbative mechanism of the TMD functions in the threshold limit and unveiling the three-dimensional picture of a hadron in the large x limit.Our formalism will be important to extract these TMD functions in the threshold region and will be reliable theoretical input to understand the experimental data.In this work, we only present the factorization and resummation formula for the unpolarized cross section, but it can be generalized to the process involving polarized hadron in both initial and final states.The corresponding theoretical predictions on the spin asymmetry in the threshold limit will be explored in future work.
where the one-loop bare soft function reads [60,79] S In addition, we replaced Q 2 / N 2 with the Collins-Soper scale ζ N .Therefore, at one-loop order (A.5) has the form as where rapidity divergences related to the poles in η cancel out during the combination and the remainder divergences in ϵ are absorbed into the renormalization factor Z uv by defining a finite collinear-soft function based on the standard RG methods.Finally, we want to point out that the perturbative results of the collinear-soft function can also be obtained by taking the N → ∞ limit in the TMD PDFs or TMD FFs [52].In the following discussion, we use TMD PDFs as our example.After performing operator product expansion, the x-space unsubtracted TMD PDF can be expressed as where C q←i is the perturbative matching coefficient, and in the moment space at one loop it can be written as where is the one-loop Altarelli-Parisi splitting kernel.In the Mellin N -space the matching coefficient can be written as Then we take the large N limit and can have where the contribution from the flavor off-diagonal splitting kernels is power suppressed in the threshold limit, and so can be neglected in the leading-power factorization formula.
To rigorously validate the factorization theorem in the joint threshold and TMD limit, examining higher-order corrections in the perturbative expansion of collinear-soft functions should be illuminating.Specifically, we employ the threshold asymptotics of the NNNLO expressions for the perturbative matching coefficients C q←q [56][57][58] to ascertain the threeloop collinear-soft function.We represent the expansion coefficients of the unsubtracted collinear-soft functions as Sunsub where the coefficients up to three loop are given by  which allows us to corroborate the RG consistency condition Γ h + 2 Γ Sc + 2 Γ fq = 0 up to the third-loop level.Below, we also provide, for convenience, the anomalous dimensions for the hard function and the threshold PDF [80]:

B Fitting parameters for PDFs and FFs
In this appendix, we will collect parametrization functional forms of PDF and FF at the initial scale µ F , where we consider CT18 PDF sets [73] and JAM20 FF sets used in section 3.In CT18 PDF sets [73], the valence quarks is parameterized as The best-fit values of the parameters a k have been given in [73], and we collect them in Table 1.Besides, we also provide their Hessian uncertainties which have been included in our numerics as the theoretical uncertainty estimation.In order to evaluate the Hessian uncertainties, we apply the following formula The fitting results are also collected in Table 1 as errors of parameters a k .For the pion FFs, we apply JAM20-SIDIS sets, and the corresponding parametrization functional form reads with initial scales are µ F = 1.3 GeV.The best fitting parameters are summarized in Table 2, where we give the error of a k corresponding uncertainties of FFs evaluated by the replica method.

+ b 2 T /b 2 max,
to avoid the Landau pole in the non-perturbative region as b T → ∞, which is b * ≡ b T » 1 with µ b * = b 0 /b * and b max = 1.5 GeV −1 .

Figure 1 .
Figure 1.Mellin inversion contour.We choose ϕ = π/2 and c = 1.6 to adapt various Q in our numerical analysis as shown in the blue dashed line.

Figure 2 .
Figure 2. Transverse momentum distribution of TMD PDFs for the up (u) quark in the proton, obtained for three different Bjorken scales x = 0.1, 0.4, 0.8 at Q = 5 GeV.The green curves correspond to the results obtained without threshold resummation.The red and blue curves represent the results obtained with threshold resummation using the ζ * scheme and without it, respectively.The collinear PDF sets used in this analysis are parameterized according to the CT18NNLO prescription, with the parameterized equations specified at the initial scale µ F = 1.3 GeV.

Figure 3 .
Figure 3. Transverse momentum distribution of threshold-TMD PDFs for u (upper panel) and d (lower panel) quark in the proton for three different hard scales and three different Bjorken scales.Here, we choose parameterized CT18NNLO collinear PDF sets where the parameterized equations are given at initial scale µ F = 1.3 GeV.

Figure 4 .
Figure 4. Transverse momentum distribution of threshold-TMD FFs for the π + from u (upper panel) and d (lower panel) quark.Here, we choose parameterized JAM20 collinear FF sets at a given scale µ F = 1.3 GeV.

Figure 5 .
Figure 5. Cross section distribution for Drell-Yan and e + e − processes.The left panel is for Drell-Yan dilepton production in pp collisions, p + p → γ* → ℓ + + ℓ − + X.The right panel is for back-to-back two hadron production in e + e − collisions, e + + e − → π + + π + + X, where we choose both hadrons to be π + as an example.The e + e − annihilation result is presented for BELLE kinematics.All the parameters for the numerical computations are mentioned in the main section of the numerical results.

Figure 6 .
Figure 6.Cross section distribution for SIDIS process, e + p → e + π + + X, where we choose the final state hadron to be a π + meson as an example.The upper panel is for HERMES (left) and COMPASS (right) kinematics, and the lower panel is for future EIC (left, 5 × 41 beam energies for ep collisions with √ s = 29 GeV) and JLab 12 GeV program (right).
≡ ln ν 2 /ζ N .We can immediately confirm that the above expressions satisfy the standard Collins-Soper equations.Upon factorizing the soft function, as previously outlined in Eq. (A.5), we procure the properly subtracted collinear-soft function Sc (b T , µ, ζ N ).Furthermore, the non-cusp anomalous dimensions are as follows: