Disentangling long and short distances in momentum-space TMDs

The extraction of nonperturbative TMD physics is made challenging by prescriptions that shield the Landau pole, which entangle long- and short-distance contributions in momentum space. The use of different prescriptions then makes the comparison of fit results for underlying nonperturbative contributions not meaningful on their own. We propose a model-independent method to restrict momentum-space observables to the perturbative domain. This method is based on a set of integral functionals that act linearly on terms in the conventional position-space operator product expansion (OPE). Artifacts from the truncation of the integral can be systematically pushed to higher powers in ΛQCD/kT. We demonstrate that this method can be used to compute the cumulative integral of TMD PDFs over kT≤kTcut\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {k}_T\le {k}_T^{\mathrm{cut}} $$\end{document} in terms of collinear PDFs, accounting for both radiative corrections and evolution effects. This yields a systematic way of correcting the naive picture where the TMD PDF integrates to a collinear PDF, and for unpolarized quark distributions we find that when renormalization scales are chosen near kTcut\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {k}_T^{\mathrm{cut}} $$\end{document}, such corrections are a percent-level effect. We also show that, when supplemented with experimental data and improved perturbative inputs, our integral functionals will enable model-independent limits to be put on the non-perturbative OPE contributions to the Collins-Soper kernel and intrinsic TMD distributions.


Introduction
Transverse momentum-dependent parton distribution functions (TMD PDFs) that encode the longitudinal and transverse partonic structure of hadrons enter factorization theorems for many physical processes at hadron colliders in a universal way. Their knowledge is essential to the ongoing research program at the Large Hadron Collider (LHC) searching for hints of new physics in precision tests of the Standard Model, see e.g. refs. [1][2][3][4][5][6][7][8][9][10][11][12][13], but TMD PDFs are also of great interest on their own as they provide insight into the strongly coupled dynamics of QCD that bind the proton together [14][15][16][17].
The momentum-space TMD PDF f i (x, ⃗ k T ) encodes the distribution of parton i in an energetic hadron at a longitudinal momentum fraction x and transverse momentum ⃗ k T . It is also often helpful to consider its description in the Fourier-conjugate position space, f i (x, ⃗ b T ). A common challenge in the extraction and interpretation of TMD PDFs is that the most intuitive way to think of them is as a distribution in ⃗ k T space, while the evolution of the TMD PDFs and their nonperturbative structure is mathematically most transparent in ⃗ b T space. An additional complication is that physical TMD cross sections in momentum space, where experimental data is taken, involve a convolution of at least two TMD PDFs f i (x, ⃗ k T ) (or TMD fragmentation or jet functions) due to momentum conservation between the two scattered partons and the color-singlet probe, see ref. [18] for a detailed review of the classic TMD processes and refs. [19][20][21][22][23][24] for examples of more recent applications with jets. This again only turns into a transparent product in position space. Another important momentum-space quantity is the cumulative distribution of a single TMD PDF over transverse momenta | ⃗ k T | ≤ k cut T . For k cut T much larger than the QCD confinement scale Λ QCD , the naive expectation is that this integral should include partons of "all" possible transverse momenta in the hadron, and therefore recover the corresponding longitudinal parton distribution function. Interestingly, the formal theoretical status of this relation is unclear to date, and clarifying it is a major goal of this paper. We review the current status of the relation in section 1.1.
The formal advantage of a position-space description of TMD PDFs is that point by point in b T ≡ | ⃗ b T |, one may formally show that it receives contributions from the running coupling α s (µ) exclusively at energy scales µ ≳ 1/b T . This makes it straightforward to delineate regions in b T where the description of the transverse structure of the hadron is perturbative, b T ≪ 1/Λ QCD , or necessarily nonperturbative, b T ≳ 1/Λ QCD , as illustrated in the left panel of figure 1. The same distinction is much harder to make in momentumspace, because at any given value of ⃗ k T , contributions from all scales µ are present. A perturbative calculation of the position-space TMD PDF f pert i (x, ⃗ b T ) can be described by an operator-product expansion (OPE), whose leading term involves a longitudinal momentum convolution between collinear densities f j (y) and a perturbative kernel C ij that is expanded in the strong coupling α s (µ), schematically f pert i (x, ⃗ b T ) = C ij ⊗ f j (x). (See eq. (1.9) below for a complete formula.) From figure 1 it is clear that the perturbative calculation of the position-space TMD PDF f pert i (x, ⃗ b T ) (dashed blue), becomes ill-behaved in the region 1/b T ∼ Λ QCD near the Landau pole µ ∼ Λ QCD of QCD.
A common way to achieve a description of the TMD PDF f i/h across all values of b T is to evaluate the perturbative part of the TMD PDF f pert i/h at some b * (b T ) instead of b T [25], where b * (b T ) → b max < 1/Λ QCD for b T → 1/Λ QCD . This causes the coupling to saturate at a perturbative scale, and avoids hitting the Landau pole. The method for implementing this cutoff is not unique, and a number of different b * prescriptions have been proposed in the literature [25][26][27][28]. The perturbative part of the TMD PDF is then multiplied with a nonperturbative model function F NP b * (b T ) that behaves as 1+O(Λ 2 QCD b 2 T ) in the perturbative region. Thus in this setup the TMD PDF is written as where µ is the renormalization scale, and ζ is the Collins-Soper scale. 1 This result can then be Fourier transformed to compute the desired quantity in momentum space. This procedure is used, for example, in the global fit analyses for TMDs [38,39] which fit for a number of parameters in their F NP b * model functions. A problem with this setup is that it does not cleanly separate the perturbative and non-perturbative contributions to the TMD. In particular the interpretation of fit results for F NP b * intrinsically depends on the choice of the b * cutoff prescription, as we have indicated by the subscript. Thus, in mixing up perturbative and nonperturbative effects, b * (or similar) prescriptions make it challenging to directly compare nonperturbative model functions obtained in different fits to experimental data.
This conclusion is a bit counterintuitive. Physical intuition about the behavior of QCD and power counting tell us that at sufficiently large transverse momenta | ⃗ k T | ≫ Λ QCD (or after integrating up to k cut T ≫ Λ QCD ), the result for f i (x, ⃗ k T ) should be dominated by the short-distance region b T ≪ 1/Λ QCD . In particular, deviations of O(Λ 2 QCD b 2 T ) from the perturbative result should have the scaling O(Λ 2 QCD /| ⃗ k T | 2 ) expected from dimensional analysis (or a direct momentum space OPE). However, when taking a standard Fourier transform of f pert i (x, ⃗ b T ) over all values of ⃗ b T , a purely perturbative baseline is not even defined due to the presence of the Landau pole or, equivalently, the need for a b * prescription. A main goal of this work is to develop a method for treating momentum-space TMDs where this intuition is restored.
In this paper, we give a set of integral functionals that only require the TMD PDF on an interval 0 ≤ b T ≤ b cut T < 1/Λ QCD as input, where the perturbative description is sufficient, and prove that for large | ⃗ k T | ≫ 1/b cut T these functionals asymptote to the full momentum-space spectrum. Our construction is based on placing a hard cutoff b cut T on the Fourier integral, which makes the contributing regions in b T space completely transparent, as illustrated in the right panel of figure 1. We then show that the result can be systematically improved by including higher-order contact terms on the boundary. Importantly, by pushing the functional to higher orders in 1/(k T b cut T ) in this way, the effect of nonperturbative corrections of O(b 2 T ) in b T space can be isolated, and the correction in momentum space can be shown to be linear in the coefficient. This provides a completely model-independent way to compute the leading effect of nonperturbative dynamics on momentum-space quantities, while the evolution and convolution of TMD PDFs can still comfortably be performed in position space at b T ≤ b cut T in our setup. We note that a hard b T -space cutoff on Fourier integrals has previously been considered in refs. [40,41] as a check on model-specific results, and the setup there amounts to the leading term 1 For definiteness, in this paper we will only consider rapidity renormalization schemes in which the hard factor H(Q 2 , µ) in the factorized TMD cross section in eq. (4.5) is renormalized in MS and independent of additional parameters. This fixes the renormalization scheme for the product of the two TMD PDFs, and fixes the scheme for individual TMDs that depend on one rapidity renormalization parameter ζ, up to rescalings of ζ a,b that leave ζaζ b = Q 4 invariant. This class of equivalent TMDs includes the definition of Collins [18], the EIS scheme [29,30], the CJNR scheme [31,32], and various others, see ref. [33] for further details. Our notation excludes the schemes in refs. [34][35][36], which depend on an additional parameter in the hard factor and TMD PDF. However these schemes can be perturbatively matched to the TMDs considered here, see [37] for a derivation of the all-orders relation. Figure 1. Left: A perturbative calculation of the TMD PDF (dashed blue) blows up at the Landau pole µ ∼ 1/b T ∼ Λ QCD . Right: In our construction in this paper, we introduce a cutoff on the Fourier integral at b cut T , below which the perturbative calculation is still a good approximation of the full TMD PDF (red). We then show that by supplementing the result with higher-order contact terms on the boundary, a purely perturbative baseline for the Fourier transform at transverse momenta | ⃗ k T | ≫ Λ QCD can be computed and systematically supplemented by higher-order nonperturbative corrections. in our construction. We also show that our construction immediately carries over to the cumulative distribution, and allows us to produce model-independent predictions for the transverse-momentum integral of the TMD PDF. As a byproduct of our analysis, we analyze the scaling of artifacts induced by b * prescriptions in momentum space, including analyzing schemes whose artifacts scale with different powers or formally fall faster than any power.
Finally, we propose an application to Drell-Yan transverse-momentum (q T ) spectra measured at the LHC. We are motivated here by the fact that the most recent measurements by CMS and ATLAS [42,43] have reached sub-percent precision, making them far more precise than other data sets sensitive to TMD physics. Model-based global TMD PDF fits typically only consider the subdominant experimental uncertainty when computing goodness-of-fit tests, which carries the risk of overconstraining certain flavor combinations and x regions. At the same time, the LHC data sets are strongly dominated by perturbative physics Λ QCD ≪ q T ≪ Q. For this reason it is desirable to perform a modelindependent fit to LHC data that explicitly incorporates this hierarchy by restricting to the leading nonperturbative effect, which is made possible by our setup. This approach is complementary to model-based fits. As an immediate physics goal, the fit we propose would directly constrain the O(b 2 T ) coefficient in the universal Collins-Soper kernel that governs the evolution of TMD PDFs with energy. Existing model-dependent fits to data [38,39] and extractions on the lattice [44][45][46][47][48][49] both prefer a surprisingly small value for this coefficient, and it would be valuable to confront the extremely precise LHC data with this observation in a model-independent way. This paper is structured as follows: In section 1.1, we review the status of the relation between TMD and longitudinal PDFs, and point out why considering the cumulative distribution for the TMD PDF is helpful. In section 2 we introduce the truncated integral functionals that allow us to compute Fourier transforms in terms of test function data on 0 ≤ b T ≤ b cut T , both for transverse-momentum spectra and cumulative distributions. We also discuss the consequences of our construction for b * prescriptions. In section 3, we apply our setup to the cumulative distribution of a single (unpolarized) TMD PDF and compare the result to the one-dimensional longitudinal PDF. In section 4 we outline the proposed application of our setup to Drell-Yan transverse-momentum spectra at the LHC. We conclude in section 5.

Review of TMD PDFs and their cumulative integral
An interesting question is how a TMD PDF f i (x, ⃗ k T ), which encodes information of both the longitudinal momentum fraction x and the transverse momentum ⃗ k T , relates to the corresponding collinear PDF f i (x), which only encodes the longitudinal information. For definiteness, we will consider f i to be the distribution of unpolarized partons within an unpolarized hadron or nucleon in this paper. The extension of this review and of our results to polarized distributions is straightforward, and we return to it in section 5. In a naive interpretation of both distributions as probability densities, one might expect i.e., the integral of the more differential TMD PDF should exactly recover the one-dimensional distribution. Using the standard identity d 2 ⃗ k T e i ⃗ b T · ⃗ k T = (2π) 2 δ( ⃗ b T ), one can equivalently ask about the value of the Fourier transform f i (x, b T ) of the TMD PDF at b T = 0, where ⃗ b T is the position-space variable conjugate to ⃗ k T and b T = | ⃗ b T |. 2 Indeed, it is easy to work out from their operator definitions [50] that eq. (1.3) is formally satisfied at the level of the bare (unrenormalized) TMD and collinear PDF, where dimensional regularization d = 4 − 2ϵ is used to treat ultra-violet (UV) divergences, and ζ denotes the so-called Collins-Soper scale, which appears due to the need to regulate and subtract rapidity divergences in the construction of the TMD PDF. Specifically, the bare position-space TMD PDF operator is defined as a correlator of parton fields that are separated by some lightlike (longitudinal) distance and by ⃗ b T in the transverse direction, and that are connected to each other by a staple-shaped Wilson line extending to lightlike infinity. For b T = 0, only the lightlike separation remains and the Wilson line collapses onto a single straight path between the two fields, which precisely recovers the operator definition of the collinear PDF. (Similarly, an additional soft component in the definition of the TMD PDF collapses to unity for b T → 0, and the dependence on the rapidity regulator drops out.) Eq. (1.4) has been explicitly verified to hold in momentum space at one-loop level in ref. [51] by performing the ⃗ k T integral in 2 − 2ϵ dimensions. However, after renormalization, the naive result eq. (1.2) must be broken, where µ is the MS scale. This is immediately clear from comparing the renormalization group evolution of the TMD and collinear PDF with respect to µ, Importantly, the µ evolution of the TMD PDF is diagonal both in flavor i and momentum fraction x, while for the PDF it sums over all flavors j and involves a convolution in x. Even more strikingly, only the TMD PDF depends on the CS scale ζ through the Collins-Soper kernel γ i ζ , also known as the rapidity anomalous dimension, while the PDF is independent of ζ. More generally, the two-dimensional structure of the renormalization in the TMD case, with the two directions tied together by a closure condition involving the cusp anomalous dimension Γ i cusp , 8) has no analog in the collinear case. Clearly, these observations forbid a simple relation between the TMD and collinear PDF at the renormalized level, and thus break eq. (1.2) in general. Only the renormalized distributions can be extracted from measurements in d = 4 dimensions, so this is unsatisfactory. At distances b T ≪ 1/Λ QCD that are short compared to the QCD confinement scale Λ QCD , another relation between the position-space TMD and collinear PDFs exists. In this regime, TMD PDFs obey an operator product expansion in terms of local operators at b T = 0, i.e., in terms of collinear PDFs [52], Here the matching coefficients C ij carry the dependence on the short-distance physics at the scale µ ∼ 1/b T . They can be computed perturbatively in terms of α s (µ) ≪ 1 and precisely absorb the mismatch in renormalization to the working order in Λ QCD .
Because the total integral of the TMD PDF over all ⃗ k T is a UV quantity that one expects to be dominated by large, perturbative momenta, one may be tempted to evaluate eq. (1.9) at b T = 0 to relate the integral to the collinear PDF. In this case, one option is to take b T → 0 while holding the renormalization scale µ fixed. This expresses the TMD PDF in terms of the PDF at the scale µ and perturbative corrections starting at O[α s (µ)]. However, the perturbative corrections contain logarithms, ln(b T µ), that prohibit setting b T = 0 to compute the total integral. Alternatively, one can evaluate the TMD PDF at µ 0 ∼ 1/b T and use the evolution equations in eqs. (1.6a) and (1.7) to evolve it back to the overall µ. In this case the perturbative corrections at the scale α s (µ 0 ) actually vanish for b T → 0 due to asymptotic freedom. However, this produces a relation with the PDF at asymptotically large scales µ 0 → ∞, which is hard to interpret.
A physically intuitive way of dealing with these issues, which also produces a practically useful result at finite scales, is to instead consider the cumulative TMD PDF in momentum space up to some finite cutoff, i.e., the distribution integrated over all ⃗ k T with | ⃗ k T | ≤ k cut T , where we can think of k cut T as an additional UV cutoff for high-energy modes. This quantity is straightforward to compute in terms of the position-space TMD PDF, Here we inserted the Fourier transform on the first line, performed the trivial integral over the angle of ⃗ b T on the second line, resulting in a zeroth-order Bessel function J 0 of the first kind. On the third line we interchanged the order of integration to perform the k T integral, resulting in the first-order Bessel function. Eq. (1.10) is a typical example of a momentum-space quantity that receives contributions from all b T , including a long-distance contribution from b T ≳ 1/Λ QCD . Unlike the full integral over all k T , it depends on a physical scale k cut T , which can guide the choice of the overall scales µ, ζ that the TMD PDF is evolved to, and the behavior of the cumulative distribution for large k cut T can be studied explicitly. In ref. [53], cumulative distributions of the unpolarized, helicity, and transversity TMD PDFs at µ = k cut T were studied for a specific choice of nonperturbative model at long distances (with perturbative results at next-to-leading order), and found to be consistent with their collinear counterparts. A primary goal of this paper, which originally motivated developing the general formalism we present in section 2, is to compute the cumulative distribution in eq. (1.10) at perturbative k cut T ≫ Λ QCD in a model-independent way in terms of the leading OPE contribution in b T space and quantifiable corrections.

Truncated functionals
Motivated by the physical applications discussed in section 1, in sections 2.1-2.4 we setup a mathematical formalism for systematically obtaining momentum space distributions from an underlying Fourier distribution that obeys an OPE. We then consider the relation to b * prescriptions in section 2.5 and the nature of an asymptotic series which arises in the construction in section 2.6.

Setup and assumptions
We consider an integral functional for an arbitrary function f that satisfies the following two assumptions: Here we take the ≲ sign to mean that there exists a constant K > 0 and an exponent ρ > 1/2 such that the bound |f (b T )| < Kb −ρ T is satisfied for all sufficiently large b T , and similarly for the first derivative with |f ′ (b T )| < Kb −ρ−1 T . In our application, S[f ](k T ) represents the radial momentum spectrum of any position-space function f (b T ) = f (| ⃗ b T |) that respects azimuthal symmetry, see the first two lines of eq. (1.10). The physical interpretation of assumption (2) is that correlation functions of parton fields (and their spatial variations) should vanish far outside the hadron. If f is not oscillating but becomes monotonic at large b T , the second part of assumption (2) concerning f ′ follows from the first. 3 We introduce a cutoff, b cut T , which divides the total functional S[f ] into two parts: In our applications, b cut T < 1/Λ QCD should be chosen such that for S < [f ] perturbation theory can be used to compute f (b T ≤ b cut T ) at short distances, implying that nonperturbative physics only becomes important in the long-distance contribution S > [f ].
We now focus on developing a formalism to obtain S[f ](k T ) for perturbative momenta k T ≫ 1/b cut T ∼ Λ QCD . Our goal is to approximate the total functional S[f ](k T ) using only information from the region b T ≤ b cut T . To the zeroth order, we define In our application, this ensures that all information about the perturbative domain is already included in the zeroth-order functional, including in particular the resummed dependence on any high scales such as the invariant mass of a color-singlet final state or the Collins-Soper scale. 4 In order to account for the leading contribution from S > [f ], we compute the integral Starting with the first exact equality, we used integration by parts and the fact that the surface term for b T → ∞ vanishes by assumption (2) in eq. (2.2). To obtain the last line we used the asymptotic form of the Bessel J 1 function at large x = k T b cut T ≫ 1 to restrict the surface term to the current working order, We further used that assumptions (1) and (2) together imply that f ′ (b T ) is bounded by for all b T > b cut T for some suitable K ′ > 0, such that we can power count away the remaining integral, 5 This scaling again follows from the asymptotic form in eq. (2.7), where the oscillations of the cosine term suppress the value of the integral in eq. (2.8) by one power of 1/(k T b cut T ). With this leading correction from S > [f ] at hand, we can define 4 An alternative would be to construct approximate Fourier functionals in terms of the behavior of the function as bT → 0, which has previously been considered in the applied mathematics literature [54], see the note added at the end of this manuscript. For our use case, this would require including also the resummed perturbative information only as part of an asymptotic series, and would produce impractical results in terms of the strong coupling and collinear PDFs at asymptotic scales µ → ∞ as discussed in section 1.1, so we do not consider this option further. 5 The first inequality is in fact incorrect in general if f ′ changes sign in a correlated fashion with J1.
The correct proof proceeds via a Mellin transformation of f and f ′ , as well as a more explicit evaluation of Bessel integrals in the asymptotic limit, and is given in section 2.2.
which is an improved approximation of S[f ]. Importantly, the leading contribution from long distances to the spectrum at perturbative k T is obtained from the value of the function at the boundary b T = b cut T .

Full proof and higher-order results for the long-distance contributions
We wish to define an asymptotic series of functionals S (n) [f ] for S[f ] which have successively improving error terms for k T b cut T → ∞. In these functionals, we aim to only use information about the test function from the region b T ≤ b cut T , but systematically account for the contributions from S > [f ] by including higher-order surface terms. In order to prove the form of the surface terms up to and including order n, we need the function f (b T ) to satisfy stronger assumptions: Assumption (2 * ) again follows from the simpler assumption and its derivatives become monotonic power laws (or exponentials) for sufficiently large b T . The reverse is always true, i.e., df follows from assumption (2 * ) for all 0 ≤ k ≤ n by bounding the iterated integral.
Mellin transform. Our strategy is to decompose f (b T ) in terms of basis functions of the form b −N T , with N > ρ a continuous parameter, to make the power counting of individual terms transparent. This amounts to taking the Mellin transform, so that the function f (b T ) can be written using the inverse Mellin transform Here the real part c of the chosen contour satisfies Together with bounded variations of f , eq. (2.13) is sufficient for the existence of the Mellin transform and its inverse, see e.g. theorem 28 in ref. [55]. Both conditions are satisified by assumption (2 * ) for n = 0, 1 and c < ρ. Below we will need c > 1/2, and this is why we required a strict inequality ρ > 1/2. Note that assumption (2 * ) in the same way guarantees the existence of the (inverse) Mellin transform for all derivatives f (k) with k ≤ n − 1. The relation is particularly compact for repeated derivatives with respect to ln b T , (2.14) Bessel integral for a single power. To make contact with the Bessel integral above a cutoff x = k T b cut T ≫ 1, let us consider The integral exists for all complex parameters β with ℜ(β) < −1/2 and 1 F 2 is the hypergeometric function. To obtain the last equality we used the asymptotic expansion of Assembling the long-distance contribution. The above definitions allow us to compute the long-distance contribution to S[f ] in terms of the Mellin transform ϕ of f , On the second line, we plugged in the inverse Mellin transform. On the third line, we changed the order of integration and used the definition of g in eq. (2.15), which is justified because the b T integral exists for each power b 1−N T with ℜ(−N ) = −c < 1/2 along the contour individually. Crucially, on the last line we performed the asymptotic expansion of g(−N, k T b cut T ) and pulled the phase and the additional power suppression of each term out of the Mellin integral. Note that naively, the scaling of the error term would be (b cut T ) −n because the remaining N integral scales as (b cut T ) −c ≤ (b cut T ) −1/2 . However, the coefficient of this additional factor of (b cut T ) −1/2 , which is encoded in ϕ(N ), is unrelated to k T and instead comes from the power law Kb −ρ T bounding the function. We now use that the coefficients c k (−N ) are polynomials in N , and evaluate the remaining Mellin integral in eq. (2.17) using eq. (2.14), which turns powers of N into repeated derivatives on the boundary, This is our key result in this section: At each order n in the power counting, the longdistance contribution can be expressed entirely in terms of the test function and its derivatives evaluated on the boundary. Adding the short-distance contribution We call these functionals of f "truncated functionals", referring both to the truncation of the b T integral and of the asymptotic series. We reiterate that the truncated functional Writing out the derivatives with respect to ln b cut T , we find for the first few coefficients 0 ≤ k ≤ 2, as required for n ≤ 3: The first term reproduces the S (1) [f ] given in eq. (2.9).
In the left panel of figure 2, we test the asymptotic expansion for a simple, but physically relevant toy test function Here the λ 1 term is a proxy for the perturbative Sudakov exponential appearing in TMDs, and the λ 2 term is a proxy for a nonperturbative contribution. The full spectrum S[f toy ] is shown by the dashed black line, while the solid lines display the absolute value of the error left over from subsequent approximations by our truncated functionals S (n) [f toy ]. We clearly see that compared to the full spectrum S[f toy ], the error term of the truncated functional exactly follows the expected scaling given by eq. (2.20) and improves as we go to higher n. T in the full spectrum that were expanded away by making use of the "perturbative" prediction f

Marrying the OPE to the truncated functionals
Consider a function f (b T ) which, apart from satisfying assumptions (1 * ) and (2 * ) in eq. (2.10), can also be written in terms of an asymptotic series expansion, where f (j) (b T ) ∼ Λ j QCD and has at most a logarithmic dependence on b T . In our physics application, the above asymptotic series arises from an operator product expansion (OPE) that expresses the matrix element of a nonlocal operator in terms of local operators suppressed by increasing powers of a low physical scale Λ QCD . We keep the Wilson coefficients of the local operators that carry the non-analytic dependence on b T implicit in f (j) . The OPE may in addition have a finite radius of convergence b OPE T , such that the limit m → ∞ exists for all b T < b OPE T , but we do not require this in the following. In physical terms this means we allow for the violation of quark-hadron duality [57] e.g. by terms ∼ exp −1/(Λ QCD b T ) . For simplicity we will also refer to the first OPE term f (0) (b T ) as the "perturbative" part of f , since in our application it can be evaluated using pertubation theory and known, simpler matrix elements.
Suppose we know the OPE up to order m − 1, which means that the leading error term is ∼ (Λ QCD b T ) m . We now show that in this case, we can approximate the total functional S[f ] using the truncated functionals with only the known OPE terms m−1 j=0 b j T f (j) as input. For simplicity let us first assume that only the leading perturbative term f (0) (b T ) is known, and that the next term is suppressed by two powers of , as is often the case. After the Fourier transform, we have correc- where on the first line we account for the correction from truncating the functional, and on the second line the three types of power corrections account for the correction from truncating the OPE. We see that the third class of power corrections, which is a combination of the neglected (b cut in the surface term and the error term of the functional, is suppressed more strongly than the intrinsic OPE correction (Λ QCD /k T ) 2 for functionals with sufficiently high order n > 5 2 . In particular, there are no corrections of O[(Λ QCD b cut T ) 2 ], and we give a short proof of this statement for the general case below.
For an illustrative example of the effect of truncating the OPE under the functional, see the right panel of figure 2, where we test eq. (2.24) with the toy function in eq. (2.22). We put only the leading perturbative term f (0) toy into the truncated functionals, where We see that increasing the order n of the truncated functionals stops improving the result around the quadratic order, because they become subleading to the intrinsic "nonperturbative" OPE correction −λ 2 b 2 T . This means starting at n = 2, the error term of the truncated functional is "stuck" at a quadratic power law, allowing us to identify the leading OPE correction directly in momentum space. In section 4 we will exploit this to propose a model-independent method to constrain nonperturbative dynamics from momentum-space data.
Now for the more general case, suppose we know the OPE up to order m − 1, i.e.
] to find the leading momentum-space error term from the truncation of the OPE. Note that for T but also satisfies assumptions (1 * ) and (2 * ) in eq. (2.10), e.g. by applying an exponential b * prescription (see section 2.5) to b m T f (m) (b T ) and multiplying the result with a similarly constructed exponential model. Then we have where in the last equality we have used that S[f (m) * ](k T ) is independent of b cut T , and thus by dimensional analysis can only yield the intrinsic OPE correction O[(Λ QCD /k T ) m ]. We see that the second term, which accounts for corrections in (Λ QCD b cut T ) m , is suppressed compared to the first when n > m + 1 2 without having to assume that b cut is individually well defined despite the growth ∼ b j T of the integrand at large b T , since it only involves a compact domain. We could therefore exploit that the truncated functionals act linearly on the terms in the OPE, which enabled us to move the sum out of the functional. Specializing to m = 2 and f (1) = 0 as above, we recover eq.  .27) is the main result of this section. It allows us to compute a momentumspace quantity (a TMD spectrum) at large momenta purely in terms of our perturbative knowledge of the position-space OPE at short distances, without having to place any assumptions on the long-distance behavior. The approximation is fully controlled and can be systematically improved by either pushing the knowledge of the OPE to higher orders m, or by carrying surface terms in the truncated functional to higher orders n. In particular, higher-order terms in the OPE are mapped to momentum space order by order in the power counting thanks to the linearity of the functionals.

Extension to cumulative distributions
In physics applications, we are often interested in a slightly different, but related integral functional, which corresponds to the cumulative distribution over all momenta | ⃗ k T | ≤ k cut T in an azimuthally symmetric position-space function f (b T ), see eq. (1.10).
We may again introduce a perturbative cutoff b cut T and divide the K[f ] into K < [f ] and K > [f ] analogous to eq. (2.4), and apply the same procedure described in section 2.2 to derive a series of functionals These truncated cumulative functionals K (n) [f ] approximate K[f ] up to an error term given by The expressions above are analogous to eqs. (2.19) and (2.20) for the spectrum. We give the first few coefficients needed to construct K (n) [f ] with eq. (2.29) for n ≤ 3: Higher-order functionals can be computed using the coefficients given in appendix A.2.
Restricting the input to the functionals to the first m − 1 orders in the OPE of the test function at small b T , we have which should be compared with eq. (2.27). Once again the linearity of the functional K[f ] was used to pull the sum outside the functional. The left panel of figure 3 shows that the cumulative functionals exhibit the expected power-law scaling of error terms for the toy function in eq. (2.22). In the right panel of figure 3 we again restrict the input of the functionals to the leading term at small b T , toy ], and compute the difference to the full result K[f toy ]. Starting from n = 2, we find that the intrinsic O(b 2 T ) correction results in a power-law behavior ∼ k cut T of this difference, as expected. Note that for the cumulative distributions, n > m − 1 2 rather than m+ 1 2 is already sufficient to expose the quadratic behavior of the intrinsic OPE correction.

Consequences for b * prescriptions
In the b T -space formalism of TMD PDF modeling, one often introduces a function b * (b T ) to shield the Landau pole in the perturbative computation, as presented in eq. (1.1) in the introduction. Two common b * prescriptions in the literature are the Collins-Soper prescription [25] and the Pavia prescription [26]: where we note that the first subleading terms in the series differ, with b 3 T and b 5 T scaling respectively.
We can formally assess the momentum-space effect of the residual b max dependence by applying our truncated functionals at sufficiently high orders to f (b * (b T )). For example, using S (n) with n ≥ 3 for CS and n ≥ 5 for Pavia, we can show using eq. (2.27) that the leading correction terms in eq. (2.34) carries over to momentum space with the scaling ∼ 1/(b max q T ) n expected from naive power counting, where n = 2 for the Collins-Soper prescription and n = 4 for Pavia. As we push to higher orders of S (n) , we can identify the scaling of the leading corrections in momentum space using a plot similar to figure 2.
To test the impact of different b * choices we consider S[σ] where σ is the leading (double) logarithmic (LL) inclusive q T differential cross section for pp → Z/γ * at the LHC, normalized to the tree-level cross section. We use the LL solution of the running coupling starting from α s (m Z ) = 0.118, which exhibits a Landau pole at Λ LL QCD ≈ 100 MeV. For numerical stability, we keep the scale of the PDFs in the resummed prediction fixed at the high scale µ f = Q = m Z , which is an NLL effect. In figure 4, we show the scaling of the difference between using b A max = 1.87 GeV −1 and b B max = 1.40 GeV −1 in momentum space for different b * models. We also include the value of an individual S[σ] (dashed gray curve) using b * Pavia with b max = 1.87 GeV −1 as reference. All curves are normalized to the Born cross section σ LO . We consider the behavior of b * CS (dotted black) and b * Pavia (dashed orange), as well as an exponential b * prescription (solid blue) defined as (2.35) The benefit of the b * exp prescription is that it does not introduce additional powers beyond those predicted by the OPE. In addition, for b T ≤ ab max any function including the prescription is exactly equal to the input function, f (b * (b T )) = f (b T ). In our numerics we have chosen a = 1/2, which is the largest value allowed such that b * exp is still monotonic. We observe that the scaling of the b max ambiguity of the CS and Pavia prescriptions both stabilize to their corresponding power law (q −2 T and q −4 T , respectively), as expected from our analytic arguments below eq. (2.34). For our newly introduced b * exp prescription we see that the scaling indeed falls off faster than any power for large enough q T , as expected. Unfortunately, we find that the coefficient of the exponentially suppressed ambiguity is rather large, so that the practical usefulness of the prescription, beyond the formal improvement of preserving the OPE coefficients, may be limited.
We note that these conclusions also apply when the b * prescription is only implemented through the choice of boundary scales in resummed perturbation theory (a so-called local b * prescription), rather than also substituting the (global) kinematic dependence of the factorized cross section by b T → b * (b T ), which is needed e.g. in the presence of additional measurements on soft radiation to preserve refactorization relations for the relevant matrix elements [28]. It also applies to the choice of TMD PDF boundary scale µ OPE = b 0 /b T + ∆ with ∆ = 2 GeV at which the OPE of the TMD PDF is performed in ref. [38], and which we expect to translate to linear power corrections O(∆/q T ) in momentum space. In either of these cases, the total b T derivative in eq. (2.34) is replaced by partial derivatives ∂/∂µ X with respect to the relevant scale(s) µ X , but the power-law scaling of the artifacts induced by the prescription still follows from expanding the relevant profile scale function µ X (b T ) for small b T , and will carry over to momentum space in the way described above. We also note that in both the setup of ref. [28] and ref. [38], the modification of the function occurs in a region where left-over fixed-order logarithms of b T µ X may be large, so the µ X dependence is not guaranteed to decrease as the perturbative order increases.

Series stability near the Landau pole
As a final check on our setup, we consider the stability of the truncated functional series when applied to a test function that exhibits an actual Landau pole beyond the support of the integral functionals. This is relevant because in the presence of a pole at b T = b pole T , we generically expect the derivatives of the test function to grow factorially, e.g., As the highest derivative f (n) enters the truncated functional at order n with coefficient (−1) n , see appendix A, this factorial growth feeds through into the result for the functional. We use the same numerical setup as in the previous setup, evaluating a single TMD PDF (or the Drell-Yan cross section) at leading-logarithmic (LL) order with a Landau pole at Λ QCD ≈ 100 MeV, but without applying any b * prescription in this case. We continue to fix the PDF scales at the high scale µ f = √ ζ for the TMD PDF and at µ f = Q for the cross section, for numerical stability.
In figure 5 LL results for the cumulant of a single quark TMD PDF (left) and the Drell-Yan transverse momentum spectrum at the LHC (right) are shown. In both cases, while the overall coefficient of the error term is substantially larger than for the test function in figures 2 and 3, we still observe a stable improvement in the error term when going to higher orders in the functional. This indicates that at least up to the order n = 3 we consider here, the Landau pole does not deteriorate the quality of the asymptotic approximation.

Evolution and operator product expansion of the TMD PDF
The dependence of the TMD PDF on the renormalization scale µ and the CS scale ζ is governed by the evolution equations in eqs. (1.6a) and (1.7). They can be solved in closed form to evolve the TMD PDF from some initial values (µ 0 , ζ 0 ) to (µ, ζ), Here we have chosen the evolution path as (µ 0 , ζ 0 ) → (µ, ζ 0 ) → (µ, ζ). We can in particular use eq. (3.1) to evolve the TMD PDF away from the canonical boundary condition µ 0 , √ ζ 0 ∼ 1/b T at which the TMD PDF is a single-scale quantity and thus free of large logarithms to all orders in perturbation theory and in Λ QCD b T . In explicit numerics we choose the canonical boundary scales as where we include a conventional O(1) numerical coefficient b 0 = 2 exp −γ E ≈ 1.12292 with γ E the Euler-Mascheroni constant. We will also refer to this boundary condition simply as the "intrinsic" TMD PDF. 7 We are interested in values of 1/b T , µ 0 , µ ≫ Λ QCD where the TMD PDF and the rapidity anomalous dimension γ i ζ obey an operator product expansion [18,52,[59][60][61], Most often the terms with odd j vanish. The µ anomalous dimension in eq. (1.6a) is a UV quantity and does not receive nonperturbative corrections, so we simply wrote down the OPE at the overall scale µ, which does not need to be equal to µ 0 . On the other hand, the ζ evolution mixes powers in the b T expansion. Specifically, for the first nonvanishing orders in a TMD PDF with f The structure of γ i ζ at higher orders is relatively well understood. The leading perturbative term γ (0) ζ,i is known to three loops [62][63][64], and also contributes the overall additive µ 7 Other choices of the boundary scale and evolution path are also possible; for examples relevant to recent TMD PDF extractions see e.g. ref. [26], where √ ζ0 is held fixed at a low scale Q0 = 1 GeV, or refs. [27,58], where the boundary scales are set by a recursive procedure taking the (perturbative and nonperturbative) rapidity anomalous dimension as input. We stress that while the result for the rapidity anomalous dimension is unambiguous, different choices of ζ0 amount to moving nonperturbative contributions proportional to the rapidity anomalous dimension in and out of the intrinsic TMD PDF. renormalization of γ i ζ in eq. (1.8). Recently, ref. [61] demonstrated that the O(b 2 T ) contribution can be related to a gluon vacuum condensate. Assuming a tree-level Wilson coefficient and ignoring the running of the condensate as in ref. [61], γ QCD is a single fixed number with units of GeV 2 .
The entries f (j) i of the OPE for a TMD PDF in general are given by multiple convolutions of perturbative Wilson coefficients, which encode the logarithmic b T dependence as b T → 0, and higher-power matrix elements that provide the power suppression by Λ j QCD . The leading term for the unpolarized quark or gluon TMD PDF takes the form are the matching coefficients, and the sum runs over all partons j. The linear correction f (1) i = 0 to the unpolarized TMD PDF vanishes by azimuthal symmetry [59,60].
Not much is known about the OPE of the unpolarized TMD PDF beyond the linear order. A subset of terms that multiply the leading-twist collinear PDF f i (x, µ) was determined to all powers in b T and at tree level in α s in ref. [65]. These correspond to target mass corrections to the TMD PDF because only the hadron (or nucleon) mass M j can provide the power suppression at O(b j T ) in this case; nevertheless, they only constitute a subset of the full OPE at any given order in b T . 8 For the purposes of this paper, we will consider the most general form that the quadratic correction can have based on the RG properties, where Λ (2) i ∼ Λ 2 QCD has units of GeV 2 . To see this, recall that the µ dependence of the quadratic term has to be that of the leading-power one, which motivates writing it as a product involving f (0) . By eq. (3.4), the remaining coefficient Λ Evaluating Λ (2) i at our reference boundary scale √ ζ 0 = b 0 /b T , we are left with a function of the flavor i, x, and b T . As for the rapidity anomalous dimension, we assume leadingorder Wilson coefficients and ignore the running of the higher-twist matrix elements and 8 Interestingly, we find that the Fourier transform of the result in ref. [65] for the unpolarized TMD PDF in their eq. (4.13) actually exists (in the distributional sense), with the cumulative distribution given by In particular, the kT dependence of the TMD PDF in this approximation is cut off at k 2 T ≤ x(1 − x)M 2 , which has a simple interpretation in terms of the remnant of the struck hadron consisting of on-shell states, i.e., (P − k) 2 ≥ 0 with P µ the hadron and k µ the parton momentum. It would be interesting to investigate whether the results of ref. [65] for spin-dependent TMD PDFs can be understood in a similar physical way. soft condensates from Λ QCD to 1/b T . (The running in general would dress the leading scaling ∼ b 2 T with an anomalous exponent, and similarly for Wilson coefficients that start at O[α s (1/b T )].) This leaves us with a function Λ (2) i (x, b T , ζ 0 ) = Λ (2) i (x) of a single variable for each flavor.
Combining both sources of O(b 2 T ) nonperturbative corrections and fixing ζ 0 , we have where we defined L ζ ≡ ln(b 2 T ζ/b 2 0 ) as a shorthand. We also absorbed the leading-power rapidity evolution factor exp 1 2 γ i , which evolves it back to the overall ζ. Importantly, the impact of γ (2) ζ,i on the evolved TMD PDF is linear at this order thanks to the power expansion of the evolution kernel. This is crucial, because our truncated functionals are precisely set up to act linearly on terms in the b T -space OPE, so we will be able to move the coefficient γ (2) ζ,i out of the Fourier transform. This setup will allow us to study the linear impact of the OPE coefficients Λ (2) i and γ (2) ζ,i on momentum-space observables, either as an uncertainty on the transverse-momentum integral of a single TMDPDF in section 3.2 or as a signal template for our proposed fit to data in section 4.

Approximating the cumulative TMD PDF with truncated functionals
The cumulative integral of an azimuthally symmetric TMD PDF over | ⃗ k T | ≤ k cut T , can be computed in terms of the b T -space TMD PDF f i with exactly the K functional we have defined in eq. (2.28). At large k cut T ≫ Λ QCD we can therefore use the formalism described in section 2 to approximate the cumulative distribution of a single TMD PDF in order to eventually study its normalization.
Thanks to the OPE eq. (3.5), we can calculate the resummed leading OPE term f (0) i in b T space in terms of collinear PDFs. We use the numerical implementation of the matching coefficients and anomalous dimensions in SCETlib [66] as also used in ref. [67]. More details on our evolution and resummation setup for f (0) i are given in section 3.3. We evaluate f (0) i at N 3 LL order, which requires the two-loop matching coefficients for the TMD PDF [62,68,69], the three-loop collinear quark anomalous dimension [70], the perturbative Collins-Soper kernel to three loops [63,64], and the four-loop cusp anomalous dimension [71,72], We assume n f = 5 active flavors in the matching coefficients and throughout the TMD evolution. We use the CT18NNLO PDF set [73] as provided by LHAPDF [74] together with four-loop running [75,76] of the coupling starting from α s (m Z ) = 0.118.
Applying our truncated functionals to the leading OPE term, i ], we obtain a purely perturbative baseline for the k T -cumulative TMD PDF. In the following we use n = 3, which we find to be numerically stable and sufficient for our purposes in this section. We set b cut T = b 0 /(1.2 GeV) such that all scales in the resummed TMD PDF stay above 1.2 GeV, i.e., well away from the Landau pole. We use standard adaptive integration Figure 6. The normalized deviation of the transverse-momentum integral (with an upper limit) of the TMD PDF from the collinear PDF as a function of the transverse-momentum cutoff k cut T (left) and of the momentum fraction x (right). The blue band shows the uncertainty from varying the position-space integral cutoff b cut T , and the yellow band shows the uncertainty from varying the quadratic OPE coefficients γ (2) ζ,q and Λ (2) q . We also include the cumulative TMD PDF as predicted by SV and Pavia global fits, both of which show small deviation from the collinear PDF. We approximate the transverse-momentum integral of the TMD PDF with i ], which is purely perturbative. In both figures, the evolution scales are set to be µ = √ ζ = k cut T , and the TMD PDF is for a down quark. See figure 12 in appendix B for the results for other light quark flavors including u,d,ū, and s. tools to evaluate K (0) [f ], and evaluate the required derivatives for the surface terms in a numerically stable fashion using Lanczos' method of differentiation-by-integration [77,78]. We find it convenient to change variables to ln b T in the test function before differentiating at ln b cut T , which both improves the numerical stability and makes it straightforward to push the implementation to high orders using the coefficients in appendix A. When needed, e.g. when comparing to explicit models, we perform Bessel integrals over the full range 0 < b T < ∞ using a double-exponential method for oscillatory integrals [79][80][81] as implemented in SCETlib [66].
We can now compare our approximated cumulant to the collinear PDF in order to quantify the deviation from our naive expectation eq. (1.2). In the left panel of figure 6, we show the normalized deviation of K (3) i ] from the collinear PDF as a function of the transverse-momentum cutoff k cut T while keeping x = 0.01 fixed. In the right panel of figure 6, we show the same deviation from the collinear PDF as a function of the momentum fraction x, while keeping k cut T = 10 GeV fixed. In both cases we set µ = √ ζ = k cut T . We find that the central value of our prediction (solid black) only deviates from the collinear PDF at most at the percent level for k cut T ≥ 5 GeV, and for x = 0.01 is essentially equal to the PDF for k cut T ≥ 10 GeV. For x = 10 −3 (x = 0.1) we find a small positive (negative) deviation.
We consider two sources of nonperturbative uncertainty on our result in this section, deferring an estimate of perturbative uncertainties to section 3.3. The first uncertainty component ∆ cut , shown in blue in figure 6, is estimated by varying b cut T from b 0 /(1.5 GeV) to b 0 /(1 GeV). Since the residual b cut T dependence is intrinsically beyond our working order, it parametrizes our ignorance of the long-distance behavior of the TMD PDF. The second source of nonperturbative uncertainty is the behavior of the TMD PDF at short to intermediate distances ≲ b cut T , where the TMD PDF is governed by the OPE. We can assess the uncertainty related to this region in a fully model-independent way by treating the quadratic OPE coefficients γ (2) ζ,i and Λ (2) i in eq. (3.9) as unknowns and propagating them into the cumulative distribution in momentum-space by injecting the corresponding OPE term into our truncated K (3) functional. Importantly, the action of the functional on each term in brackets in eq. (3.9) is linear, i.e., we have, suppressing the arguments of the TMD PDF, (3.11) The impact of γ (2) ζ,i and Λ (2) i on the predicted transverse-momentum integral of the TMD PDF is therefore also linear, including their a priori unknown sign. While we consider simple estimates of the unknown coefficients here to demonstrate the method, their linear impact makes it completely straightforward to update the uncertainty estimate (or the central value) to any determination of the coefficients from model-based global fits to experimental data or a lattice calculation. It also makes the uncertainty very cheap to evaluate since the three functionals in eq. (3.11) only need to be evaluated once, rather than having to perform a Fourier transform for each choice of model parameters.
For definiteness, we consider a variation of Λ (2) q (x) ∈ [−1 GeV 2 , 1 GeV 2 ] independent of x, which yields an uncertainty ∆ f NP . We stress again that this assumption is straightforward to relax point by point in the right panel of figure 6 by a trivial rescaling. To estimate the uncertainty due to higher OPE corrections to the CS kernel, we adopt a very conservative range of γ (2) ζ,q ∈ [−(0.5 GeV) 2 , (0.5 GeV) 2 ], yielding the second OPE-related uncertainty component ∆ ζ NP ; recent model-based global fits [38,39] and the condensate-based estimate from ref. [61] all point to substantially lower values of |γ (2) ζ,q | ≲ 0.05 GeV 2 ≈ (0.22 GeV) 2 . For readability, both nonperturbative uncertainty components are added in quadrature in figure 6, ∆ ζ NP ⊕ ∆ f NP ≡ (∆ ζ NP ) 2 + (∆ f NP ) 2 1/2 , and are shown in yellow. We find that for x = 0.01 the integral of the TMD PDF is fully compatible with the collinear PDF for all values of k cut T already within the nonperturbative uncertainties, while for general x the nonperturbative uncertainty alone only covers part of the difference. Both uncertainty components quickly decrease as k cut T increases, as expected from their respective power counting.
Finally we also include in figure 6 the cumulative TMD PDF as predicted by the SV [38] and Pavia [39] model-based global TMD fits. We use arTeMiDe [82] to evaluate the SV fit result in b T space, assuming the best-fit N 3 LL model parameters from ref. [38], and use our in-house setup to perform the Bessel transform. In the case of the Pavia global fit, we encounter the difficulty that our treatment of quark mass thresholds is different from the one in the NangaParbat [83] native fit code, where the number n f of active massless flavors is changed discontinuously as a function of b T . This has a very large effect at large b T (about 4% at b 0 /b T ≈ 2 GeV, and up to 10% at b 0 /b T ≈ 1 GeV), which distorts the comparison of the genuine nonperturbative effect. For this reason we have set up an independent implementation of the Pavia fit model, using the central replica of the N 3 LL Pavia fit as provided in ref. [84], and combined it with our implementation of f (0) i in SCETlib evaluated with the b * Pavia prescription of ref. [26], with b min → 0 as in eq. (2.33). We have verified that above the bottom-quark threshold, our result for f (0) is within a permil of the result returned by NangaParbat for b min = 0. 9 Both global fits support our conclusion that the deviation from the naive expectation is within percent level. In particular, they are compatible with our perturbative baseline well within our estimate of the nonperturbative uncertainty, indicating that our model-independent approach is sound.
In figure 12 in appendix B, we show the same plots as in figure 6 for other light quark flavors (u,d,ū, and s). We find that the results are very similar to the discussion for the down quark in this section.

Perturbative uncertainty
Since we rely on the accuracy of the perturbative input f (0) i , it is also important to consider the theoretical uncertainty from the truncation of perturbation theory. We estimate the perturbative uncertainty at different resummation orders by varying the boundary scales entering the evolved TMD PDF. As described below, we adapt the variation procedures developed for the resummed q T spectrum in ref. [67], noting that for our case there are some simplifications compared to the original setup.
The TMD PDF at generic ζ, which may be parametrically separate from ζ 0 ∼ 1/b 2 T , is a multi-scale problem. One way of treating it in resummed perturbation theory is to evaluate it at µ 0 ∼ √ ζ 0 ∼ 1/b T , where it is free of large logarithms, and use eq. (3.1) to evolve it back to the overall (µ, ζ). The uncertainty estimate in that case could be estimated by varying the boundary scales (µ 0 , ζ 0 ).
Alternatively, one can consider a factorization of the TMD PDF at generic ζ into purely collinear and purely soft matrix elements, (3.12) which individually are single-scale problems and feature an additional rapidity renormalization scale ν. In the SCET literature, these matrix elements are known as beam and soft functions. 10 In this setup, the running in ν from the canonical beam scale ν B ∼ √ ζ to the soft scale ν S ∼ 1/b T resums the large logarithms of ζb 2 T . For the central prediction, where all scales are set exactly to their canonical values, this is exactly equivalent to running the TMD PDF as a whole as described below eq. (3.1). The two approaches differ, however, 9 The discrepancy to arTeMiDe in the perturbative region, or when turning off the arTeMiDe nonperturbative model, is somewhat larger, at roughly half a percent. We attribute this to the different choice of boundary scale within the so-called ζ prescription used in arTeMiDe. 10 We use the exponential regulator of ref. [85] that renders the overlap subtraction between the beam and soft function scaleless, see the discussion in Section 2 of ref. [33]. More details on rapidity regulators can be found in Appendix B of ref. [33]. The procedure of ref. [67] is to consider the variations 11 i ]| v i for each of them and then take the maximum envelope as our estimate of the perturbative resummation uncertainty ∆ res .
In figure 7 we show ∆ res for the cumulative TMD PDF calculated at NLL, NNLL, and N 3 LL as a function of k cut T and of x. The parameter choices are the same as in section 3.2. When we calculate the perturbative f (0) i to higher resummation orders, we see a convergence of the ∆ res bands, giving us confidence in the robustness of the uncertainty estimate. The central values at all resummation orders are within a few percent from the collinear PDF for k cut T ≳ 10 GeV, and the cumulative TMD PDF is compatible with the 11 Note that in ref. [67], hybrid profile scales were used to turn off the resummation as qT → Q to implement the matching to power corrections of O(q 2 T /Q 2 ) (also known as the Y term). For a single TMD PDF, which is already one of the objects appearing in the leading-power factorization, these corrections are absent, so we can keep the resummation on everywhere. This amounts to taking x1,2,3 → ∞ in the profile scales in ref. [67]. collinear PDF well within the perturbative uncertainty for all k cut T and x. Notice that the perturbative uncertainty in figure 7 is much larger than the nonperturbative uncertainties shown in figure 6. Figure 13 in appendix B shows the same plots as in figure 7 for other light quark flavors. The results are once again very similar to the discussion for the down quark in this section.

Impact of evolution effects
In the discussions above, we have chosen µ = √ ζ = k cut T , which is the natural scale choice such that the left-hand side of eq. (3.10) only depends on one scale k cut T . In this section, we would like to consider how evolution away from these natural scales impacts our conclusion about the transverse-momentum integral of TMD PDFs.
In figure 8, we show the normalized deviation of the cumulative TMD PDF at k cut T = 10 GeV from the collinear PDF as a function of both the renormalization scale µ and the CS scale ζ. We find that, in the vicinity of the natural scale choice µ = √ ζ = k cut T , our conclusion is robust against evolution effects and the deviation remains within percent level. In particular, along the line of fixed µ = k cut T , we find that the ζ evolution is negligible. This has the intriguing physical interpretation that the ζ evolution of the momentumspace TMD PDF is only a shape effect, i.e., while it changes the shape of the TMD PDF at low k T , the transverse-momentum integral stays the same and is determined by collinear physics.
When we move away from µ = k cut T , sizable deviations from the collinear PDF arise from the large logarithm induced by the cusp anomalous dimension. To further understand this, consider the evolution path (µ a , ζ a ) = (k cut T , k cut T ) → (µ a , ζ b ) → (µ b , ζ b ). Since the ζ evolution at fixed µ a = k cut T is only a shape effect, the first evolution step only changes the transverse-momentum integral of the TMD PDF within a few percent. In the second step, however, the µ evolution factor resums large double logarithms of the form where the ellipses indicate single logarithms and Γ cusp [α s ] = αs 4π Γ i 0 + O(α 2 s ) is the one-loop coefficient of the cusp anomalous dimension. Thus the effect of the µ evolution is stronger for a larger ratio of √ ζ b /µ a at the starting point, which is responsible for the saddle point structure observed in figure 8. This confirms the initial observation in section 1 that an equality between cumulative TMD PDF and collinear PDF cannot hold at all scales due to the different renormalization, which resums double logarithms in the case of the TMD PDF. Nevertheless, a relation in the vicinity of µ, √ ζ ∼ k cut T is not in contradiction with the renormalization, and is in fact observed in figure 8 to also extend to parametrically separate values of ζ at fixed µ = k cut T . We conclude that, even when we take the evolution effects in the vicinity of µ = k cut T into account, our intuitive result that the transverse-momentum integral of TMD PDF is within few percent of the corresponding collinear PDF remains robust. It is also interesting to look at the ζ evolution alone. In figure 9, we show the deviation of the transverse-momentum integral of a TMD PDF from the collinear PDF as a function of √ ζ, while keeping x = 0.01 and µ = k cut T = 10 GeV fixed. In the left (right) panel of figure 9, we consider nonperturbative (perturbative) uncertainties ∆ cut , ∆ ζ NP , and ∆ ζ NP ⊕∆ f NP (∆ res ) as described in section 3.2 (section 3.3). We see that the perturbative uncertainties from resummation still dominate over the nonperturbative uncertainties described in section 3.2, but reduce as the TMD PDF is evolved to very high energies √ ζ ∼ 1 TeV. The ∆ cut uncertainty (light blue in the left panel) likewise reduces at large √ ζ, which is expected as the Sudakov suppression of the long-distance region increases, while the ∆ ζ NP (light green) component comes to dominate the total OPE uncertainty (yellow) both at very small, ζ → 1 GeV, and very large ζ → 1 TeV, i.e., far away from the natural choice ζ ∼ k cut T . Figure 14 and figure 15 in appendix B show the same plots as in figure 8 and figure 9 for other light quark flavors. The results are very similar to the discussion for the down Figure 9. The normalized deviation of the transverse-momentum integral of the TMD PDF from the collinear PDF as a function of the evolution scale √ ζ. We also include the SV and Pavia global fits in the left panel. We approximate the transverse-momentum integral of the TMD PDF with i ], which is purely perturbative. In both figures, we use fixed µ = k cut T = 10 GeV and x = 0.01, and the TMD PDF is for a down quark. See figure 15 in appendix B for the results for other light quark flavors including u,d,ū, and s. quark here. Once again, remarkably good agreement is found for the integral of the TMD PDF and the collinear PDF with the most natural choice for fixing the evolution scales. Figure 10. Impact plots of quadratic "OPE" parameters in the toy function in eq. (2.25) interpreted as nonperturbative uncertainties (left) or fits to data (right). We consider the impact of the "nonperturbative" parameter λ 2 on the toy model in eq. (2.22) with the expansion in eq. (2.25). To produce "data" (solid lines on the right) for coefficients λ 2 > 0, where the original definition has a divergent Fourier transform, we have modified the toy function as described in the text, see eq. (4.2). We use a third-order functional S (3) throughout. The innermost uncertainty bands in light blue or gray, respectively, indicate ∆ cut , which for clarity is only shown for λ 2 = 0 in the right panel.
4 Model-independent constraints on nonperturbative physics from data

Outline of the method
In the previous section we used the action of our truncated functionals on the quadratic terms in the OPE to provide a model-independent estimate of nonperturbative uncertainties in cumulative TMD PDFs, treating the coefficients as unknowns. This approach of course remains valid for transverse-momentum cross sections. In this section we explain how information on nonperturbative contributions to cross sections can be isolated, illustrating the steps with the toy model σ = f toy in eq. (2.25). In the left panel of figure 10 we show the transverse-momentum spectrum S[f toy ], where we vary the quadratic coefficient λ 2 within two different ranges of values, as indicated by the yellow bands. We normalize the impact of λ 2 to a baseline given by the third-order truncated functional applied to the purely "perturbative" component of the function, S (3) [f (0) ]. We see that the relative impact of λ 2 is slightly larger than expected from power counting, e.g., it amounts to ∼ 2% at k T = 10 GeV for λ 2 = 1 GeV, which we attribute to the falloff of the spectrum beyond the Sudakov peak that is present also for this toy function.
There is, however, a second way to make use of the linear impact of the quadratic OPE coefficients in our setup when the transverse-momentum spectrum is known from data. In this case, a fit to data may be performed treating the OPE coefficients as a "signal" that can be discriminated from the purely perturbative "background" in a model-independent way. Importantly, since the action of the functional is linear, the fit template for the signal only needs to be computed once and can then be trivially rescaled in the fit, e.g. for our toy function in eq. (2.25), where λ 2 (the combined "OPE" coefficient at the quadratic order) would be the parameter of interest in the fit. We illustrate this approach in the right panel of figure 10, where the dashed lines are the combined fit template for signal and background for different values of λ 2 . To generate illustrative "data" (solid lines) for λ 2 < 0, where the original toy function in eq. (2.22) has a divergent Fourier transform, we modify it as where f is unchanged. This ensures that the full integral S[f toy ] exists for all λ 2 while maintaining f (2) T for either sign. We see that the fit template predicted by the truncated functionals closely follows the mock data in figure 10 for all values of k T and λ 2 , with the exception of λ 2 = −(1 GeV) 2 where our modification in eq. (4.2) results in a strong modification of f toy at large distances, leading to a larger oscillatory error term.

Proposed application to LHC Drell-Yan data
We next outline an application of our method aimed at constraining nonperturbative corrections using LHC data for the Drell-Yan process, pp → Z/γ * → ℓ + ℓ − , which are some of the most precise data taken at the LHC to date [42,43]. At leading power in q T ≪ Q, but to all orders in Λ QCD /q T , the fiducial cross section differential in the dilepton transverse momentum reads where Φ B = {Q, Y, cos θ} is the Born phase-space for a dilepton pair with total invariant mass Q and total rapidity Y , where the lepton is scattered at a polar angle θ measured in the Collins-Soper frame [86]. (For the Born kinematics considered here, the Collins-Soper frame coincides with the frame boosted by Y along the beam axis from the lab frame.) In eq. (4.3), we have made power corrections in q T /Q explicit, but will suppress these in the following. The acceptance of the fiducial volume evaluated on Born-level kinematics with q T = 0 is denoted by A B (Φ B ). For a typical fiducial volume at the LHC, where Q min,max define an invariant-mass bin (e.g. around the Z resonance) and the lepton transverse momenta p 1 T = p 2 T = sin θ Q/2 and pseudorapidities η 1,2 = Y ∓ ln tan θ/2 are required to be within the fiducial detector volume. (The extension to bins of rapidity as e.g. in ref. [42] is trivial.) In eq. (4.3) we have used our shorthand S[. . . ](q T ) for the integral functional in eq. (2.1) that performs the Fourier transform from position space. It acts on the factorized unpolarized cross section in position space [87], where σ B ij (Q) is the Born cross section for the hard production process ij → Z/γ * that includes the Q-dependent Z lineshape, H(Q 2 , µ) = 1 + O(α s ) is the so-called hard function that encodes virtual corrections to the hard process, 12 and the sum runs over all contributing quark flavors. The two TMD PDFs on the right-hand side of eq. (4.5) are evaluated at with E cm the total center-of-mass energy of the proton-proton collision. Inserting eq. (3.9), it is straightforward to give the factorized unpolarized cross section to O( 12 Here we have neglected singlet contributions to the quark axial and vector form factors starting at O(α 2 s ) and O(α 3 s ), respectively, where the initial-state flavors i, j do not couple directly to the Z boson. These contributions are straightforward to restore for the high-precision background prediction, but can safely be neglected in the recasting step for a measured Λ (2) coefficient that we outline below.
where L Q 2 = ln(b 2 T Q 2 /b 2 0 ). If we try to use eq. (4.7) in eq. (4.3) directly, we face a problem commonly encountered in extractions of (TMD)PDFs, namely that the flavor and x a,b dependence of the model, in this case the Λ (2) i,j (x a,b ), only enters under the sum over flavors and the integral over hard phase space, and enters in different ways for different center-of-mass energies. Resolving this usually requires writing down a model-specific ansatz and performing a global fit to a large data set at once to constrain it. We now demonstrate that thanks to the linear impact of the OPE contribution to the TMD PDF at quadratic order, its x and flavor dependence can be condensed into a single aggregate coefficient Λ (2) per fiducial volume that can be determined by a two-parameter fit with γ (2) ζ,i to the data in that single fiducial volume. We also demonstrate that the measured central value, or limits, for Λ (2) can be straightforwardly recast and checked against specific models by performing a simple weighted tree-level integral over the model function.
To set up some more notation, let us define the total fiducial cross section in b T space, which is possible because we consider the q T -independent Born acceptance A(Φ B ). We then simply have dσ/dq T = S[σ(b T )](q T ). Our goal in the following is to break down the quadratic contribution σ (2) (b T ) to σ(b T ) into two terms σ (2,γ ζ ) + σ (2,Λ) with a single, overall coefficient each encoding the nonperturbative physics. This is straightforward for the contribution from the rapidity anomalous dimension, because it is flavor independent and thus only involves a simple logarithm of Q as an additional weight under the hard phase-space integral over the leading-power cross section.
To cast the contribution from the intrinsic TMD PDFs into the same form, we define where to match eq. (4.7) the average coefficient must be defined by . (4.11) The fact that it is simply a ratio of (weighted) leading-power integrals is due to the linear impact of Λ i,j (x a,b ) on the cross section point by point in phase space. It is always possible to evaluate the right-hand side of eq. (4.11) for a given model to compare it to a measured value of Λ (2) .
However, given the expected sensitivity to Λ (2) at the LHC (discussed below), we anticipate that the right-hand side of eq. (4.11) will only be required at relatively low accuracy to compare (limits on) Λ (2) to models. In this case, we can neglect logarithmic Q dependence across the Q bin in several places to a good approximation, namely in the hard function boundary condition, in its evolution down to the TMD PDF scale µ 0 ∼ b 0 /b T , and in the leading-power rapidity evolution factor. As a consequence, all of these terms cancel in the ratio since we can effectively evaluate them at some reference value of Q within the bin. We can also evaluate the intrinsic leading-power TMD PDF at NLL, leaving only . (4.12) Both the numerator and the denominator are a straightforward tree-level integral of collinear PDFs against the Z line shape and the acceptance, with an additional weight of Λ (2) i,j (x a,b ) in the numerator. The residual dependence on µ 0 can be used as a handle on the theory uncertainty during this recasting step, where In total, our b T -space input for a given fiducial volume integrated over hard phase space reads Applying our truncated functionals at order n = 3 and exploiting their linearity, we have, (4.14) We now illustrate the numerical impact of the two remaining nonperturbative parameters, working at fixed Q = m Z and Y = 0 and applying no fiducial cuts for simplicity, i.e., we consider dΦ B = dQ dY d cos θ The recasting relation in eq. (4.12) in this scenario reduces to a weighted flavor average at We plot the impact of γ where in the first two rows we again interpret the impact of the parameters as a modelindependent nonperturbative uncertainty estimate and compare it to two recent modelbased global fits [38,39] as described in section 3.2. The top row shows our uncertainty estimates as bands around our central result for the Drell-Yan spectrum. The middle row shows the corresponding relative uncertainty. We find that both global fits are roughly compatible with the perturbative baseline already within the truncation uncertainty ∆ cut , i.e., setting both nonperturbative parameters to zero. This result implies that the fits report a vanishing total effect of (quadratic) nonperturbative parameters for central LHC kinematics. The impact of these terms is thus smaller than one would estimate purely from dimensional analysis [61]. Figure 11. Impact plots of quadratic OPE parameters in the LHC Drell-Yan cross section interpreted as nonperturbative uncertainties (top and middle row) or fits to data (bottom). We consider the impact of the quadratic coefficient in the rapidity anomalous dimension (left column) and the intrinsic TMD PDF averaged over hard phase space (right column) on the physical N 3 LL cross section. Our averaging procedure for the intrinsic TMD PDF coefficient Λ (2) is defined in eq. (4.11) and discussed there. We generate mock data (solid lines in bottom row) as described around eq. It would be valuable to check this conclusion with a complementary model-independent method. This is illustrated in the plots in the bottom row of figure 11 to which experimental data would need to be added to determine the favored values. We plot the signal template (dashed lines) for different values of the nonperturbative parameters of interest with γ (2) ζ,q on the bottom-left and Λ (2) on the bottom-right, taking the respective other one to be zero. In lieu of using experimental data, we produce corresponding mock data for the left-hand side of eq. (4.14), shown by solid lines in the bottom row of figure 11. To do so we use a b * Pavia prescription on σ (0) , which as demonstrated in section 2.5 modifies the cross section beyond our working order in eq. (4.14). We then combine this with a nonperturbative model function as in eq. (4.2) that flips sign as needed and recovers the respective quadratic coefficient, where λ = −Λ (2) and k = 0 for the intrinsic TMD PDF, and λ = −γ (2) ζ,q and k = 1 for the rapidity anomalous dimension.
Comparing the dashed and solid lines in figure 11, we find that the proposed fit template approximates the model results well down to q T ∼ 8 GeV before nonlinearities and oscillatory error terms kick in, which suggests q T ≥ 8 GeV as a possible lower bound for the fit window. This is similar in spirit to searches for nonzero d ≥ 6 SMEFT coefficients at the LHC, which also often involve restricting the data set to data below (here: above) a certain cut to ensure that the fit stays in the linear region and that limits remain model independent. We stress again that the signal is linear in the nonperturbative rapidity anomalous dimension after the power expansion, even though it appears in the exponent when working to all orders in b T space. We anticipate that this will simplify breaking the (likely) degeneracy between the rapidity anomalous dimension and the intrinsic TMD PDF by performing the above fit to at least two fiducial volumes since the universal parameter γ (2) ζ,q enters linearly in both.

Overview of theory systematics
From figure 11, we can read off that to be sensitive to |γ (2) ζ,q | ≲ (0.35 GeV) 2 and |Λ (2) | ≲ (0.7 GeV) 2 , we must maintain control over the perturbative baseline at a precision level of well below 2%. For definiteness we will consider a precision goal of 1% on the perturbative baseline in the following and review the effects relevant at this level of precision. The experimental systematics and statistical uncertainties on the normalized Drell-Yan spectrum are negligible at this scale [42,43], so we do not discuss them further. Many items on this list require a dedicated implementation effort, or even are not yet understood theoretically. For this reason, a complete fit to actual LHC data is beyond the scope of this paper.
QCD corrections to σ (0) (Q, Y, b T ). Based on the perturbative uncertainty estimates in ref. [67], the N 3 LL perturbative order we consider here for illustration are at the 3-4% level, and hence not sufficient to reach the 1% goal. The complete three-loop QCD boundary conditions for the leading-power factorized cross section in eq. (4.13) have been computed in the meantime [88][89][90] and enable the resummation at N 3 LL ′ accuracy, which we expect to lower the perturbative QCD uncertainty by a factor of ≈ 2 [91][92][93][94].
Power corrections to TMD factorization. The factorization in eq. (4.3) receives power corrections from several sources. The first are power corrections of O(q T /Q) from the full q T dependence of the acceptance, dubbed "fiducial" power corrections in ref. [67]. They are logarithmically divergent at small q T and thus have to be evaluated in resummed perturbation theory, but this is straightforward since the linear power corrections can be shown [67] to multiply the factorized cross section in eq. (4.7), where the effect of γ (2) ζ,q and Λ i,j (x a , x b ) can be dropped as higher order cross terms. In addition, there are genuine QCD corrections of O(q 2 T /Q 2 ), often referred to as non-singular or Y -term corrections. These can be incorporated by a standard fixed-order matching, implying that the fit window would not be restricted from above by a validity range for TMD factorization like q T /Q ≲ 0.25.
QED effects. Incorporating the effect of QED initial-state radiation into the factorization in eq. (4.5) is well understood [95][96][97]. It was also shown in ref. [67] that QED final-state radiation preserves the form of eq. (4.3) and only leads to O(α em ) corrections to the Born acceptance, as long as the experimental lepton definition is sufficiently inclusive. The key effect of QED is thus to break the factorization of the cross section into a leptonic and hadronic part. Very little is known about these initial-final interference (IFI) contributions in the regime of small transverse momentum. In the vicinity of the Z resonance they can be expected to scale as ∼ α em Γ Z /q T in order to recover the narrow-width approximation at large q T . (Since the leptons distinguish a transverse direction in this case, one can no longer appeal to azimuthal symmetry to push the correction to O(b 2 T ).) Understanding the interplay of IFI effects with the all-order resummation would thus be an important improvement for the fit we are proposing.
Quark mass effects. In our analysis we used n f = 5 massless active flavors throughout. This assumption will break down at small q T , close to the fit window we examined here. Corrections to this limit from massive quarks Q have been calculated [98] that affect both the rapidity anomalous dimension and the TMD PDF boundary condition at O(α 2 s m Q b 2 T ), where m Q = m c , m b is the respective quark mass. (We refer the reader to ref. [98] for a review of earlier work on quark mass effects in the TMD PDF boundary condition.) The importance of treating the quark mass dependence of the CS kernel correctly was also stressed in ref. [59]. As these effects have the same scaling b 2 T as our signal, they must be accounted for in the resummed background prediction in order to not pollute our fit. We remind the reader that in section 3.2, we found large numerical differences to the NangaParbat code used in ref. [39], where the number of massless active flavors was changed dynamically as a function of b T . Like our approach of using a fixed n f = 5, this cannot be a complete description of the region between m c ≲ 1/b T ≲ m b : Both approaches treat all quarks as either completely massless or fully decoupled, and thus they are both missing the exact (nonlogarithmic) dependence on b T m c,b , i.e., an O(1) effect in this region. Nevertheless, taking the two approaches to be the extreme choices, the spread between them can serve as an indicator of the size of these mass effects, and we conclude that it will be critical to correctly account for them in a resummed prediction in the future.

Conclusions
In this work, we developed a formalism to disentangle the long-and short-distance physics in TMD objects in momentum space. This task is challenging because taking the Fourier transform of position-space TMD PDFs mixes up contributions from perturbative shortdistance physics with the models used to treat the Landau pole and to describe the nonperturbative long-distance behavior.
We introduced a position-space cutoff b cut T < 1/Λ QCD in the Fourier transform, such that the momentum-space TMDs at perturbative momenta k T ≫ Λ QCD manifestly depend only on the short-distance physics. We then defined two series of truncated integral functionals, S (n) [f ] for the momentum-space spectrum and K (n) [f ] for the cumulative momentum distribution, which systematically account for the long-distance contributions from the region beyond b cut T . The power corrections can be calculated as contact terms on the boundary of integration, and only depend on the information of the integrand at b cut T . In the short-distance region, where we calculate our integral functionals, TMD PDFs can be perturbatively calculated using an operator product expansion (OPE). Pairing the truncated functionals with the perturbative TMDs computed from the OPE, we constructed a model-independent, systematically improvable perturbative baseline for momentum-space TMD quantities and presented explicit formulas for all sources of uncertainty in the approximation.
As an important application, we showed how to determine the transverse-momentum integral of TMD PDFs by combining the OPE with the K (n) functionals to compute their cumulative distribution up to a transverse momentum cutoff k cut T . This allows us to quantitatively assess the naive expectation that the integral over transverse momentum of a TMD PDF equals the corresponding collinear PDF. We first considered the natural choices for the renormalization and Collins-Soper scales µ = √ ζ = k cut T , performing the renormalization group evolution of the TMD PDF to these scales in b T space at N 3 LL order, and found that the deviation of the integral of the TMD PDF from the collinear PDF is at the percent level for perturbative momenta k cut T ≳ 10 GeV. In particular, the deviation as a function of the momentum fraction x is within ∼ 2% for 10 −3 ≤ x ≤ 10 −2 . These N 3 LL results in particular include the two-loop radiative corrections to the TMD PDF as encoded in the leading OPE, and we find that the cumulative distribution is always compatible with the collinear PDF within our estimate of the perturbative uncertainty. Using our truncated functional construction, we can also assess the nonperturbative uncertainty on our results in a model-independent way, which in particular includes a conservative estimate of the leading nonperturbative contribution from the Collins-Soper kernel. We find that the agreement of the integral of the TMD PDF with the collinear PDF deteriorates for general x if only the nonperturbative uncertainty is considered, motivating both an extension to higher perturbative orders in the future to push down the dominant perturbative uncer-tainty and dedicated model-independent studies of the leading (quadratic) contributions to the Collins-Soper kernel and the intrinsic TMD PDF.
Evolving the TMD PDF to different overall scales µ, ζ, we find that the deviation remains within a few percent in the vicinity of µ = √ ζ = k cut T . Furthermore, the evolution of ζ at fixed µ = k cut T is found to only be a shape effect, meaning that while it changes the shape of the TMD PDF at low k T , the cumulative integral stays the same and is determined by the collinear PDF. We tested our conclusions for the unpolarized TMD PDFs for all light quark flavors.
We also proposed to apply our S (n) functionals to momentum-space spectra for the Drell-Yan process, which have been measured to sub-percent precision at the LHC, in order to develop a method of putting model-independent constraints on the leading nonperturbative effects in TMD PDFs. In our approach, the O(b 2 T ) nonperturbative contributions to the Collins-Soper kernel, γ (2) ζ,q , and the intrinsic TMD PDF, Λ i (2) , are treated as a "signal", which can be discriminated from the perturbative "background" by comparing to experimental data. Since the functionals are linear in the OPE coefficients (the "signal strengths"), we are optimistic that the proposed fit to data will be greatly simplified. The method we propose is model-independent and complementary to model-based global fits. We briefly discussed the major theory systematics in the proposed fit, such as QCD corrections to the perturbative cross section, power corrections to the TMD factorization, QED effects, and quark mass effects, all of which influence the perturbative precision. These systematics are important to consider because our estimates suggest that the background prediction must be pushed to roughly a percent accuracy in order to be able to see the nonperturbative signal for a generic O[(1 GeV) 2 ] size of the coefficients. Our results in this paper lend themselves to several future applications. The truncated functional formalism that we developed is generally applicable to Fourier transforms of spherically symmetric functions f (|⃗ r|) (also known as Hankel transforms) in any number of dimensions, with minimal modifications from the two-dimensional case that we presented. Our results are useful in scenarios where the test function is only known on a radial interval 0 < |⃗ r | < r cut , and provide a systematic expansion of the Fourier transform at large ⃗ k in inverse powers of r cut | ⃗ k| ≫ 1. In practical applications, e.g. in the fit we proposed in section 4, one could consider averaging over several values of the cut parameter (b cut T in our case) to smoothen the residual oscillations, which is equivalent to using a smoother function to switch off the integrand around the central b cut T . The numerical stability of the functionals might also be further improved by moving part of the boundary term under the integral as a subtraction term.
For physics applications to TMDs, we note that a fit template linear in the nonperturbative parameter of interest can also be obtained in modifications of our setup. For example, one could apply the exponential b * prescription introduced in eq. (2.35) to individual terms in the OPE, and multiply the result by a similarly constructed "model function", e.g. 1 − exp −(b T /b max − a) −1 . This yields a b T -space integrand that does not contain spurious power corrections, but does fall off as 1/b T at large b T and thus can be integrated over all 0 ≤ b T < ∞. This may be numerically advantageous due to the ability to use series acceleration for the indefinite integral.
Finally, for the study of individual TMD PDFs, an obvious next step is to study the gluon TMD PDF, which we will address in a separate publication. The gluon TMD PDF is an interesting test case because it features polarization already at leading twist even within an unpolarized hadron [99].
For polarized hadrons, it would be interesting to verify in a model-independent way the observation of ref. [53] that the cumulative helicity (g 1 ) and transversity (h 1 ) TMD PDFs are approximately equal to their respective collinear counter parts for µ = k cut T and large enough k cut T . This extension is straightforward in our framework because the OPE for the helicity and transversity TMD PDFs has the same general form as discussed in sections 1.1 and 3.1, i.e., they both have a leading contribution from (polarized) leading-twist collinear PDFs, so our setup carries over immediately. For TMD PDFs whose OPE receives its first contribution at O(Λ QCD b T ) (at subleading twist), our results imply that the cumulative distribution at large k cut T vanishes as Λ QCD /k cut T . This has interesting consequences e.g. for the Burkardt sum rule for the Sivers function [100], whose possible violation at the level of renormalized quantities can be studied in our numerical setup as a function of k cut T including in particular the effect of (rapidity) evolution. Work on this is also in progress.
Note added: While finalizing this manuscript, we became aware that related techniques for computing asymptotic expansions of integral transforms were developed in applied mathematics several decades ago, see e.g. refs. [54,101,102]. The most closely related results were obtained in ref. [54], where Bessel integrals over an interval 0 ≤ t ≤ a were computed in terms of the test function and its derivatives at t = a, and the result was then used to express full Bessel integrals over 0 ≤ t < ∞ in terms of the behavior of the function near the origin t → 0. Our approach in section 2 differs by instead using (numerical) test function data on a compact interval 0 ≤ b T ≤ b cut T to approximate the full integral, which is better suited for our physics applications. This is the case because the analytic form of derivatives of our test functions as b T → 0 is not readily accessible for TMD PDFs or cross sections once radiative and evolution effects are included.