# Disentangling Long and Short Distances
in Momentum-Space TMDs

Markus A. Ebert  Johannes K. L. Michel  Iain W. Stewart  and Zhiquan Sun

(January 18, 2022)

###### Abstract

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 $\Lambda_{\rm QCD}/k_{T}$. We demonstrate that this method can be
used to compute the cumulative integral of TMD PDFs over $k_{T}\leq
k_{T}^{\mathrm{cut}}$ 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 $k_{T}^{\mathrm{cut}}$, 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 nonperturbative OPE contributions to the
Collins-Soper kernel and intrinsic TMD distributions.

##  1 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,{\vec{k}}_{T})$ encodes the distribution
of parton $i$ in an energetic hadron at a longitudinal momentum fraction $x$
and transverse momentum ${\vec{k}}_{T}$. It is also often helpful to consider
its description in the Fourier-conjugate position space,
$f_{i}(x,{\vec{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 ${\vec{k}}_{T}$ space, while the evolution of the TMD
PDFs and their nonperturbative structure is mathematically most transparent in
${\vec{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,{\vec{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 $\lvert{\vec{k}}_{T}\rvert\leq k_{T}^{\mathrm{cut}}$. For
$k_{T}^{\mathrm{cut}}$ much larger than the QCD confinement scale
$\Lambda_{\mathrm{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}\equiv\lvert{\vec{b}}_{T}\rvert$, one may formally show that
it receives contributions from the running coupling $\alpha_{s}(\mu)$
exclusively at energy scales $\mu\gtrsim 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}\ll
1/\Lambda_{\mathrm{QCD}}$, or necessarily nonperturbative, $b_{T}\gtrsim
1/\Lambda_{\mathrm{QCD}}$, as illustrated in the left panel of figure -299.
The same distinction is much harder to make in momentum-space, because at any
given value of ${\vec{k}}_{T}$, contributions from all scales $\mu$ are
present. A perturbative calculation of the position-space TMD PDF
$f^{\mathrm{pert}}_{i}(x,{\vec{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 $\alpha_{s}(\mu)$,
schematically $f^{\mathrm{pert}}_{i}(x,{\vec{b}}_{T})=C_{ij}\otimes f_{j}(x)$.
(See eq. (1.9) below for a complete formula.) From figure -299 it is clear
that the perturbative calculation of the position-space TMD PDF
$f^{\mathrm{pert}}_{i}(x,{\vec{b}}_{T})$ (dashed blue), becomes ill-behaved in
the region $1/b_{T}\sim\Lambda_{\mathrm{QCD}}$ near the Landau pole
$\mu\sim\Lambda_{\mathrm{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^{\mathrm{pert}}_{i/h}$ at some $b^{*}(b_{T})$ instead of $b_{T}$ [25],
where $b^{*}(b_{T})\to b_{\mathrm{max}}<1/\Lambda_{\mathrm{QCD}}$ for
$b_{T}\to 1/\Lambda_{\mathrm{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_{b^{*}}^{\mathrm{NP}}(b_{T})$ that behaves as
$1+\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2})$ in the perturbative
region. Thus in this setup the TMD PDF is written as

$\displaystyle f_{i/h}(x,b_{T},\mu,\zeta)$ $\displaystyle=f_{i/h}^{\rm
pert}(x,b^{*}(b_{T}),\mu,\zeta)\>F_{b^{*}}^{\mathrm{NP}}(b_{T})\,,$ (1.1)

where $\mu$ is the renormalization scale, and $\zeta$ is the Collins-Soper
scale.111For definiteness, in this paper we will only consider rapidity
renormalization schemes in which the hard factor $H(Q^{2},\mu)$ in the
factorized TMD cross section in eq. (4.5) is renormalized in
$\overline{\text{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 $\zeta$, up to rescalings of $\zeta_{a,b}$ that leave
$\zeta_{a}\zeta_{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.  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_{b^{*}}^{\mathrm{NP}}$ 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_{b^{*}}^{\mathrm{NP}}$ 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 $\lvert{\vec{k}}_{T}\rvert\gg\Lambda_{\mathrm{QCD}}$ (or
after integrating up to $k_{T}^{\mathrm{cut}}\gg\Lambda_{\mathrm{QCD}}$), the
result for $f_{i}(x,{\vec{k}}_{T})$ should be dominated by the short-distance
region $b_{T}\ll 1/\Lambda_{\mathrm{QCD}}$. In particular, deviations of
$\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2})$ from the perturbative
result should have the scaling
$\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}/\lvert{\vec{k}}_{T}\rvert^{2})$
expected from dimensional analysis (or a direct momentum space OPE). However,
when taking a standard Fourier transform of
$f^{\mathrm{pert}}_{i}(x,{\vec{b}}_{T})$ over all values of ${\vec{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.

Figure -299: Left: A perturbative calculation of the TMD PDF (dashed blue)
blows up at the Landau pole $\mu\sim 1/b_{T}\sim\Lambda_{\mathrm{QCD}}$.
Right: In our construction in this paper, we introduce a cutoff on the Fourier
integral at $b_{T}^{\mathrm{cut}}$, 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
$\lvert{\vec{k}}_{T}\rvert\gg\Lambda_{\mathrm{QCD}}$ can be computed and
systematically supplemented by higher-order nonperturbative corrections.

In this paper, we give a set of integral functionals that only require the TMD
PDF on an interval $0\leq b_{T}\leq
b_{T}^{\mathrm{cut}}<1/\Lambda_{\mathrm{QCD}}$ as input, where the
perturbative description is sufficient, and prove that for large
$\lvert{\vec{k}}_{T}\rvert\gg 1/b_{T}^{\mathrm{cut}}$ these functionals
asymptote to the full momentum-space spectrum. Our construction is based on
placing a hard cutoff $b_{T}^{\mathrm{cut}}$ on the Fourier integral, which
makes the contributing regions in $b_{T}$ space completely transparent, as
illustrated in the right panel of figure -299. 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_{T}^{\mathrm{cut}})$ in this way, the effect of nonperturbative
corrections of $\mathcal{O}(b_{T}^{2})$ 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}\leq b_{T}^{\mathrm{cut}}$ 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 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 $\Lambda_{\mathrm{QCD}}\ll q_{T}\ll Q$. For this reason
it is desirable to perform a model-independent 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 $\mathcal{O}(b_{T}^{2})$ 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\leq b_{T}\leq
b_{T}^{\mathrm{cut}}$, 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.

###  1.1 Review of TMD PDFs and their cumulative integral

An interesting question is how a TMD PDF $f_{i}(x,{\vec{k}}_{T})$, which
encodes information of both the longitudinal momentum fraction $x$ _and_ the
transverse momentum ${\vec{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

$\displaystyle\int\\!\mathrm{d}^{2}{\vec{k}}_{T}\,f_{i}(x,{\vec{k}}_{T})\stackrel{{\scriptstyle?}}{{=}}f_{i}(x)\,,$
(1.2)

i.e., the integral of the more differential TMD PDF should exactly recover the
one-dimensional distribution. Using the standard identity
$\int\\!\mathrm{d}^{2}{\vec{k}}_{T}\,e^{\mathrm{i}{\vec{b}}_{T}\cdot{\vec{k}}_{T}}=(2\pi)^{2}\delta({\vec{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$,

$\displaystyle\int\\!\mathrm{d}^{2}{\vec{k}}_{T}\,f_{i}(x,{\vec{k}}_{T})=f_{i}(x,b_{T}=0)\stackrel{{\scriptstyle?}}{{=}}f_{i}(x)\,,$
(1.3)

where ${\vec{b}}_{T}$ is the position-space variable conjugate to
${\vec{k}}_{T}$ and $b_{T}=\lvert{\vec{b}}_{T}\rvert$.222We will continue to
use the same symbol $f_{i}$ for the unpolarized TMD PDF, its Fourier
transform, and the unpolarized collinear PDF in the following as the meaning
will always be clear from the arguments or the context.

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,

$\displaystyle
f_{i}^{\mathrm{bare}}(x,b_{T}=0,\epsilon,\zeta)=f_{i}^{\mathrm{bare}}(x,\epsilon)\,,$
(1.4)

where dimensional regularization $d=4-2\epsilon$ is used to treat ultra-violet
(UV) divergences, and $\zeta$ 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 ${\vec{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}\to 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
${\vec{k}}_{T}$ integral in $2-2\epsilon$ dimensions.

However, after renormalization, the naive result eq. (1.2) must be broken,

$\displaystyle\int\mathrm{d}^{2}{\vec{k}}_{T}\,f_{i}(x,{\vec{k}}_{T},\mu,\zeta)\neq
f_{i}(x,\mu)\,,$ (1.5)

where $\mu$ is the $\overline{\text{MS}}$ scale. This is immediately clear
from comparing the renormalization group evolution of the TMD and collinear
PDF with respect to $\mu$,

$\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}f_{i}(x,{\vec{k}}_{T},\mu,\zeta)$
$\displaystyle=\gamma_{\mu}^{i}(\mu,\zeta)\,f_{i}(x,{\vec{k}}_{T},\mu,\zeta)$
(1.6a) $\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}f_{i}(x,\mu)$
$\displaystyle=\sum_{j}\int_{x}^{1}\frac{\mathrm{d}z}{z}P_{ij}(z,\mu)f_{j}\Bigl{(}\frac{x}{z},\mu\Bigr{)}\,.$
(1.6b)

Importantly, the $\mu$ 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 $\zeta$ through the Collins-Soper kernel $\gamma_{\zeta}^{i}$,
also known as the rapidity anomalous dimension,

$\displaystyle\zeta\frac{\mathrm{d}}{\mathrm{d}\zeta}f_{i}(x,b_{T},\mu,\zeta)$
$\displaystyle=\frac{1}{2}\gamma_{\zeta}^{i}(b_{T},\mu)\,f_{i}(x,b_{T},\mu,\zeta)$
(1.7)

while the PDF is independent of $\zeta$. 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
$\Gamma_{\mathrm{cusp}}^{i}$,

$\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}\gamma_{\zeta}^{i}(b_{T},\mu)$
$\displaystyle=2\zeta\frac{\mathrm{d}}{\mathrm{d}\zeta}\gamma^{i}_{\mu}(\mu,\zeta)=-2\Gamma_{\mathrm{cusp}}^{i}[\alpha_{s}(\mu)]\,,$
(1.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}\ll 1/\Lambda_{\mathrm{QCD}}$ that are short compared to
the QCD confinement scale $\Lambda_{\mathrm{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],

$\displaystyle
f_{i}(x,b_{T},\mu,\zeta)=\sum_{j}\int_{x}^{1}\frac{\mathrm{d}z}{z}C_{ij}\bigl{(}z,b_{T},\mu,\zeta\bigr{)}\,f_{j}\Bigl{(}\frac{x}{z},\mu\Bigr{)}+\mathcal{O}\bigl{(}\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2}\bigr{)}\,.$
(1.9)

Here the matching coefficients $C_{ij}$ carry the dependence on the short-
distance physics at the scale $\mu\sim 1/b_{T}$. They can be computed
perturbatively in terms of $\alpha_{s}(\mu)\ll 1$ and precisely absorb the
mismatch in renormalization to the working order in $\Lambda_{\mathrm{QCD}}$.

Because the total integral of the TMD PDF over all ${\vec{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}\to 0$ while
holding the renormalization scale $\mu$ fixed. This expresses the TMD PDF in
terms of the PDF at the scale $\mu$ and perturbative corrections starting at
$\mathcal{O}[\alpha_{s}(\mu)]$. However, the perturbative corrections contain
logarithms, $\ln(b_{T}\mu)$, that prohibit setting $b_{T}=0$ to compute the
total integral. Alternatively, one can evaluate the TMD PDF at $\mu_{0}\sim
1/b_{T}$ and use the evolution equations in eqs. (1.6a) and (1.7) to evolve it
back to the overall $\mu$. In this case the perturbative corrections at the
scale $\alpha_{s}(\mu_{0})$ actually vanish for $b_{T}\to 0$ due to asymptotic
freedom. However, this produces a relation with the PDF at asymptotically
large scales $\mu_{0}\to\infty$, 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 ${\vec{k}}_{T}$ with
$\lvert{\vec{k}}_{T}\rvert\leq k_{T}^{\mathrm{cut}}$, where we can think of
$k_{T}^{\mathrm{cut}}$ as an additional UV cutoff for high-energy modes. This
quantity is straightforward to compute in terms of the position-space TMD PDF,

$\displaystyle\int_{\lvert{\vec{k}}_{T}\rvert\leq
k_{T}^{\mathrm{cut}}}\\!\mathrm{d}^{2}{\vec{k}}_{T}\,f_{i}(x,{\vec{k}}_{T},\mu,\zeta)$
$\displaystyle=\int_{\lvert{\vec{k}}_{T}\rvert\leq
k_{T}^{\mathrm{cut}}}\\!\mathrm{d}^{2}{\vec{k}}_{T}\,\int\frac{\mathrm{d}^{2}\,{\vec{b}}_{T}}{(2\pi)^{2}}\,e^{+i{\vec{k}}_{T}\cdot{\vec{b}}_{T}}f_{i}(x,b_{T},\mu,\zeta)$
$\displaystyle=\int_{0}^{k_{T}^{\mathrm{cut}}}\\!\mathrm{d}k_{T}\,k_{T}\int_{0}^{\infty}\\!\mathrm{d}b_{T}\,b_{T}J_{0}(b_{T}k_{T})\,f_{i}(x,b_{T},\mu,\zeta)$
$\displaystyle=k_{T}^{\mathrm{cut}}\int_{0}^{\infty}\\!\mathrm{d}b_{T}J_{1}(b_{T}k_{T}^{\mathrm{cut}})\,f_{i}(x,b_{T},\mu,\zeta)\,.$
(1.10)

Here we inserted the Fourier transform on the first line, performed the
trivial integral over the angle of ${\vec{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.1) is a typical example of a momentum-space quantity that receives
contributions from all $b_{T}$, including a long-distance contribution from
$b_{T}\gtrsim 1/\Lambda_{\mathrm{QCD}}$. Unlike the full integral over all
$k_{T}$, it depends on a physical scale $k_{T}^{\mathrm{cut}}$, which can
guide the choice of the overall scales $\mu,\zeta$ that the TMD PDF is evolved
to, and the behavior of the cumulative distribution for large
$k_{T}^{\mathrm{cut}}$ can be studied explicitly. In ref. [53], cumulative
distributions of the unpolarized, helicity, and transversity TMD PDFs at
$\mu=k_{T}^{\mathrm{cut}}$ 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.1) at perturbative
$k_{T}^{\mathrm{cut}}\gg\Lambda_{\mathrm{QCD}}$ in a model-independent way in
terms of the leading OPE contribution in $b_{T}$ space and quantifiable
corrections.

##  2 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.

###  2.1 Setup and assumptions

We consider an integral functional

$\displaystyle
S[f](k_{T})=k_{T}\int_{0}^{\infty}\\!\mathrm{d}b_{T}\,b_{T}J_{0}(k_{T}b_{T})\,f(b_{T})\,,$
(2.1)

for an arbitrary function $f$ that satisfies the following two assumptions:

$\displaystyle(1)\quad f(b_{T})~{}\text{is continuously differentiable for
all}~{}b_{T}>0\,,$ $\displaystyle(2)\quad\lvert
f(b_{T}\to\infty)\rvert\lesssim b_{T}^{-\rho}~{}\text{and}~{}\lvert
f^{\prime}(b_{T}\to\infty)\rvert\lesssim
b_{T}^{-\rho-1},\quad\rho>\frac{1}{2}\,.$ (2.2)

Here we take the $\lesssim$ sign to mean that there exists a constant $K>0$
and an exponent $\rho>1/2$ such that the bound $\lvert
f(b_{T})\rvert<Kb_{T}^{-\rho}$ is satisfied for all sufficiently large
$b_{T}$, and similarly for the first derivative with
$|f^{\prime}(b_{T})|<Kb_{T}^{-\rho-1}$. In our application, $S[f](k_{T})$
represents the radial momentum spectrum of any position-space function
$f(b_{T})=f(\lvert\vec{b}_{T}\rvert)$ that respects azimuthal symmetry, see
the first two lines of eq. (1.1). 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^{\prime}$ follows from the first.333In addition, we assume that
the original integral in eq. (2.1) _exists_ for any $k_{T}>0$ in the sense of
an improper Riemann integral. Assumption (2) can be made sufficient to
guarantee its existence for any $k_{T}$ in the Lebesgue sense by requiring
$\rho>3/2$, but this is not required for our manipulations below. We also have
not made assumptions about the behavior of $f(b_{T})$ as $b_{T}\to 0$ yet.

We introduce a cutoff, $b_{T}^{\mathrm{cut}}$, which divides the total
functional $S[f]$ into two parts:

$\displaystyle S[f](k_{T})$ $\displaystyle=S_{<}[f](k_{T})+S_{>}[f](k_{T})\,,$
(2.3)

where

$\displaystyle S_{<}[f](k_{T})$
$\displaystyle=k_{T}\int_{0}^{b_{T}^{\mathrm{cut}}}\\!\mathrm{d}b_{T}\,~{}b_{T}J_{0}(k_{T}b_{T})f(b_{T})\,,$
$\displaystyle S_{>}[f](k_{T})$
$\displaystyle=k_{T}\int_{b_{T}^{\mathrm{cut}}}^{\infty}\\!\mathrm{d}b_{T}\,~{}b_{T}J_{0}(k_{T}b_{T})f(b_{T})\,.$
(2.4)

In our applications, $b_{T}^{\mathrm{cut}}<1/\Lambda_{\mathrm{QCD}}$ should be
chosen such that for $S_{<}[f]$ perturbation theory can be used to compute
$f(b_{T}\leq b_{T}^{\mathrm{cut}})$ 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}\gg
1/b_{T}^{\mathrm{cut}}\sim\Lambda_{\mathrm{QCD}}$. Our goal is to approximate
the total functional $S[f](k_{T})$ using only information from the region
$b_{T}\leq b_{T}^{\mathrm{cut}}$. To the zeroth order, we define

$\displaystyle S^{(0)}[f](k_{T},b_{T}^{\mathrm{cut}})\equiv
S_{<}[f](k_{T})\,.$ (2.5)

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.444An
alternative would be to construct approximate Fourier functionals in terms of
the behavior of the function as $b_{T}\to 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 $\mu\to\infty$ as discussed in section
1.1, so we do not consider this option further.  In order to account for the
leading contribution from $S_{>}[f]$, we compute the integral

$\displaystyle S_{>}[f](k_{T})$
$\displaystyle=k_{T}\int_{b_{T}^{\mathrm{cut}}}^{\infty}\\!\mathrm{d}b_{T}\,b_{T}J_{0}(k_{T}b_{T})\,f(b_{T})$
$\displaystyle=-b_{T}^{\rm cut}J_{1}(k_{T}b_{T}^{\rm cut})\,f(b_{T}^{\rm
cut})-\int_{b_{T}^{\rm
cut}}^{\infty}\\!\mathrm{d}b_{T}\,b_{T}J_{1}(k_{T}b_{T})\,f^{\prime}(b_{T})$
$\displaystyle=\sqrt{\frac{2b_{T}^{\rm cut}}{\pi
k_{T}}}\cos\Bigl{(}k_{T}b_{T}^{\rm cut}+\frac{\pi}{4}\Bigr{)}\,f(b_{T}^{\rm
cut})+\frac{1}{k_{T}}\mathcal{O}\bigl{[}(k_{T}b_{T}^{\rm
cut})^{-\frac{1}{2}}\bigr{]}\,.$ (2.6)

Starting with the first exact equality, we used integration by parts and the
fact that the surface term for $b_{T}\to\infty$ vanishes by assumption (2) in
eq. (2.1). To obtain the last line we used the asymptotic form of the Bessel
$J_{1}$ function at large $x=k_{T}b_{T}^{\mathrm{cut}}\gg 1$ to restrict the
surface term to the current working order,

$\displaystyle J_{1}(x\to\infty)=-\sqrt{\frac{2}{\pi
x}}\cos\Big{(}x+\frac{\pi}{4}\Bigr{)}+\mathcal{O}\bigl{(}x^{-\frac{3}{2}}\bigr{)}\,.$
(2.7)

We further used that assumptions (1) and (2) together imply that
$f^{\prime}(b_{T})$ is bounded by $K^{\prime}\,b_{T}^{-\rho-1}$ for _all_
$b_{T}>b_{T}^{\mathrm{cut}}$ for some suitable $K^{\prime}>0$, such that we
can power count away the remaining integral,555The first inequality is in fact
incorrect in general if $f^{\prime}$ changes sign in a correlated fashion with
$J_{1}$. The correct proof proceeds via a Mellin transformation of $f$ and
$f^{\prime}$, as well as a more explicit evaluation of Bessel integrals in the
asymptotic limit, and is given in section 2.2.

$\displaystyle\int_{b_{T}^{\rm
cut}}^{\infty}\\!\mathrm{d}b_{T}\,b_{T}J_{1}(k_{T}b_{T})\,f^{\prime}(b_{T})$
$\displaystyle\lesssim\int_{b_{T}^{\rm
cut}}^{\infty}\mathrm{d}b_{T}~{}b_{T}J_{1}(k_{T}b_{T})\,b_{T}^{-\rho-1}=\frac{(b_{T}^{\mathrm{cut}})^{-\rho}}{k_{T}}\mathcal{O}\bigl{[}(k_{T}b_{T}^{\mathrm{cut}})^{-\frac{1}{2}}\bigr{]}\,.$
(2.8)

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_{T}^{\mathrm{cut}})$.

With this leading correction from $S_{>}[f]$ at hand, we can define

$\displaystyle
S^{(1)}[f](k_{T},b_{T}^{\mathrm{cut}})=S^{(0)}[f](k_{T},b_{T}^{\mathrm{cut}})+\sqrt{\frac{2b_{T}^{\mathrm{cut}}}{\pi
k_{T}}}\cos\Bigl{(}k_{T}b_{T}^{\mathrm{cut}}+\frac{\pi}{4}\Bigr{)}~{}f(b_{T}^{\mathrm{cut}})\,,$
(2.9)

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_{T}^{\mathrm{cut}}$.

###  2.2 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_{T}^{\mathrm{cut}}\to\infty$. In these functionals, we aim to only use
information about the test function from the region $b_{T}\leq
b_{T}^{\mathrm{cut}}$, 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:

$\displaystyle(1^{*})\quad f(b_{T})~{}\text{is $n$ times differentiable for
all}~{}b_{T}>0\,,$
$\displaystyle(2^{*})\quad\frac{\mathrm{d}^{n}f(b_{T})}{\mathrm{d}b_{T}^{n}}\lesssim
b_{T}^{-\rho-n},\quad\rho>1/2\,.$ (2.10)

Assumption $(2^{*})$ again follows from the simpler assumption
$f(b_{T})\lesssim b_{T}^{-\rho}$ if $f(b_{T})$ and its derivatives become
monotonic power laws (or exponentials) for sufficiently large $b_{T}$. The
reverse is always true, i.e.,
$\mathrm{d}f^{k}(b_{T})/\mathrm{d}b_{T}^{k}\lesssim b_{T}^{-\rho-k}$ follows
from assumption $(2^{*})$ for all $0\leq k\leq 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_{T}^{-N}$, with $N>\rho$ a continuous parameter, to make the power
counting of individual terms transparent. This amounts to taking the Mellin
transform,

$\displaystyle\phi(N)=\int_{0}^{\infty}\\!\mathrm{d}b_{T}\,b_{T}^{N-1}f(b_{T})\,,$
(2.11)

so that the function $f(b_{T})$ can be written using the inverse Mellin
transform

$\displaystyle
f(b_{T})=\frac{1}{2\pi\mathrm{i}}\int_{c-\mathrm{i}\infty}^{c+\mathrm{i}\infty}\\!\mathrm{d}N\,b_{T}^{-N}\phi(N)\,.$
(2.12)

Here the real part $c$ of the chosen contour satisfies

$\displaystyle\int_{0}^{\infty}\\!\mathrm{d}b_{T}\,b_{T}^{c-1}\lvert
f(b_{T})\rvert<\infty\,.$ (2.13)

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<\rho$. Below we will need $c>1/2$, and this is why we required a strict
inequality $\rho>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\leq n-1$. The relation is particularly compact for repeated
derivatives with respect to $\ln b_{T}$,

$\displaystyle\Bigl{(}b_{T}\frac{\mathrm{d}}{\mathrm{d}b_{T}}\Bigr{)}^{k}f(b_{T})=\frac{1}{2\pi\mathrm{i}}\int_{c-\mathrm{i}\infty}^{c+\mathrm{i}\infty}\\!\mathrm{d}N\,b_{T}^{-N}\,\phi(N)\,(-N)^{k}\,.$
(2.14)

#### Bessel integral for a single power.

To make contact with the Bessel integral above a cutoff
$x=k_{T}b_{T}^{\mathrm{cut}}\gg 1$, let us consider

$\displaystyle g(\beta,x)$
$\displaystyle\equiv\int_{x}^{\infty}\\!dy\,J_{0}(y)\,y^{1+\beta}$
$\displaystyle=\frac{2^{1+\beta}\,\Gamma\bigl{(}1+\frac{\beta}{2}\bigr{)}}{\Gamma\bigl{(}-\frac{\beta}{2}\bigr{)}}-\frac{x^{2+\beta}}{2+\beta}\,\phantom{}_{1}F_{\,2}\Bigl{(}1+\frac{\beta}{2};1,2+\frac{\beta}{2};-\frac{x^{2}}{4}\Bigr{)}$
$\displaystyle=\sqrt{\frac{2x}{\pi}}x^{\beta}\sum_{k=0}^{n-1}\frac{c_{k}(\beta)}{x^{k}}\cos\Bigl{(}x+\frac{\pi}{4}-\frac{k\pi}{2}\Bigr{)}+\mathcal{O}\bigl{(}x^{\beta-n+\frac{1}{2}}\bigr{)}\,.$
(2.15)

The integral exists for all complex parameters $\beta$ with $\Re(\beta)<-1/2$
and $\phantom{}{}_{1}F_{\,2}$ is the hypergeometric function. To obtain the
last equality we used the asymptotic expansion of
$\phantom{}{}_{1}F_{\,2}\bigl{(}a_{1};b_{1},b_{2};-x^{2}/4\bigr{)}$ for $x\gg
1$ [56], which generalizes the asymptotic expansion of the Bessel function in
eq. (2.7). The coefficients $c_{k}(\beta)$ are polynomials in $\beta$ of order
$k$ and are collected in appendix A.1 up to $k=9$. The first few coefficients
read

$\displaystyle c_{0}=1\,,\qquad c_{1}(\beta)=-\frac{3}{8}-\beta\,,\qquad
c_{2}(\beta)=-\frac{15}{128}-\frac{\beta}{8}+\beta^{2}\,.$ (2.16)

#### 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 $\phi$ of $f$,

$\displaystyle S_{>}[f](k_{T})$
$\displaystyle=k_{T}\int_{b_{T}^{\mathrm{cut}}}^{\infty}\\!\mathrm{d}b_{T}\,b_{T}J_{0}(k_{T}b_{T})f(b_{T})$
$\displaystyle=k_{T}\int_{b_{T}^{\mathrm{cut}}}^{\infty}\\!\mathrm{d}b_{T}\,J_{0}(k_{T}b_{T})\frac{1}{2\pi\mathrm{i}}\int_{c-\mathrm{i}\infty}^{c+\mathrm{i}\infty}\\!\mathrm{d}N\,\phi(N)\,b_{T}^{1-N}$
$\displaystyle=\frac{1}{k_{T}}\frac{1}{2\pi\mathrm{i}}\int_{c-\mathrm{i}\infty}^{c+\mathrm{i}\infty}\\!\mathrm{d}N\,\phi(N)\,k_{T}^{N}\,g(-N,k_{T}b_{T}^{\mathrm{cut}})$
$\displaystyle=\sqrt{\frac{2b_{T}^{\mathrm{cut}}}{\pi
k_{T}}}\sum_{k=0}^{n-1}\frac{\cos\Bigl{(}k_{T}b_{T}^{\mathrm{cut}}+\frac{\pi}{4}-\frac{k\pi}{2}\Bigr{)}}{(k_{T}b_{T}^{\mathrm{cut}})^{k}}\frac{1}{2\pi\mathrm{i}}\int_{c-\mathrm{i}\infty}^{c+\mathrm{i}\infty}\\!\mathrm{d}N\,c_{k}(-N)\,\phi(N)\,(b_{T}^{\mathrm{cut}})^{-N}$
$\displaystyle\qquad+\frac{f(b_{T}^{\mathrm{cut}})}{k_{T}}\mathcal{O}\bigl{[}(k_{T}b_{T}^{\mathrm{cut}})^{-n+\frac{1}{2}}\bigr{]}\,.$
(2.17)

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.2), which is justified because the $b_{T}$ integral exists for each
power $b_{T}^{1-N}$ with $\Re(-N)=-c<1/2$ along the contour individually.
Crucially, on the last line we performed the asymptotic expansion of
$g(-N,k_{T}b_{T}^{\mathrm{cut}})$ 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_{T}^{\mathrm{cut}})^{-n}$ because
the remaining $N$ integral scales as
$(b_{T}^{\mathrm{cut}})^{-c}\leq(b_{T}^{\mathrm{cut}})^{-1/2}$. However, the
coefficient of this additional factor of $(b_{T}^{\mathrm{cut}})^{-1/2}$,
which is encoded in $\phi(N)$, is unrelated to $k_{T}$ and instead comes from
the power law $Kb_{T}^{-\rho}$ 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.2) using eq. (2.14), which
turns powers of $N$ into repeated derivatives on the boundary,

$\displaystyle S_{>}[f](k_{T})$
$\displaystyle=\sqrt{\frac{2b_{T}^{\mathrm{cut}}}{\pi
k_{T}}}\>\sum_{k=0}^{n-1}\>\frac{\cos\Bigl{(}k_{T}b_{T}^{\mathrm{cut}}+\frac{\pi}{4}-\frac{k\pi}{2}\Bigr{)}}{(k_{T}b_{T}^{\mathrm{cut}})^{k}}\,c_{k}\Bigl{(}b_{T}^{\mathrm{cut}}\frac{\mathrm{d}}{\mathrm{d}b_{T}^{\mathrm{cut}}}\Bigr{)}f(b_{T}^{\mathrm{cut}})$
$\displaystyle\quad+\frac{f(b_{T}^{\mathrm{cut}})}{k_{T}}\mathcal{O}\bigl{[}(k_{T}b_{T}^{\mathrm{cut}})^{-n+\frac{1}{2}}\bigr{]}\,.$
(2.18)

This is our key result in this section: At each order $n$ in the power
counting, the long-distance contribution can be expressed entirely in terms of
the test function and its derivatives evaluated on the boundary. Adding the
short-distance contribution $S_{<}[f]=S^{(0)}[f]$, we define

$\displaystyle S^{(n)}[f](k_{T},b_{T}^{\mathrm{cut}})$ $\displaystyle\equiv
k_{T}\int_{0}^{b_{T}^{\mathrm{cut}}}\\!\mathrm{d}b_{T}\,b_{T}J_{0}(k_{T}b_{T})\,f(b_{T})$
$\displaystyle\quad+\sqrt{\frac{2b_{T}^{\mathrm{cut}}}{\pi
k_{T}}}\>\sum_{k=0}^{n-1}\>\frac{\cos\Bigl{(}k_{T}b_{T}^{\mathrm{cut}}+\frac{\pi}{4}-\frac{k\pi}{2}\Bigr{)}}{(k_{T}b_{T}^{\mathrm{cut}})^{k}}\,c_{k}\Bigl{(}b_{T}^{\mathrm{cut}}\frac{\mathrm{d}}{\mathrm{d}b_{T}^{\mathrm{cut}}}\Bigr{)}f(b_{T}^{\mathrm{cut}})\,.$
(2.19)

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 $S^{(n)}[f]$ approximates $S[f]$ up to

$\displaystyle
S[f](k_{T})=S^{(n)}[f](k_{T},b_{T}^{\mathrm{cut}})+\frac{f(b_{T}^{\mathrm{cut}})}{k_{T}}\mathcal{O}\bigl{[}(k_{T}b_{T}^{\mathrm{cut}})^{-n+\frac{1}{2}}\bigr{]}\,.$
(2.20)

Writing out the derivatives with respect to $\ln b_{T}^{\mathrm{cut}}$, we
find for the first few coefficients $0\leq k\leq 2$, as required for $n\leq
3$:

$\displaystyle
c_{0}\Bigl{(}b_{T}^{\mathrm{cut}}\frac{\mathrm{d}}{\mathrm{d}b_{T}^{\mathrm{cut}}}\Bigr{)}f(b_{T}^{\mathrm{cut}})$
$\displaystyle=f(b_{T}^{\mathrm{cut}})\,,$ $\displaystyle
c_{1}\Bigl{(}b_{T}^{\mathrm{cut}}\frac{\mathrm{d}}{\mathrm{d}b_{T}^{\mathrm{cut}}}\Bigr{)}f(b_{T}^{\mathrm{cut}})$
$\displaystyle=-\frac{3}{8}f(b_{T}^{\mathrm{cut}})-b_{T}^{\mathrm{cut}}f^{\prime}(b_{T}^{\mathrm{cut}})\,,$
$\displaystyle
c_{2}\Bigl{(}b_{T}^{\mathrm{cut}}\frac{\mathrm{d}}{\mathrm{d}b_{T}^{\mathrm{cut}}}\Bigr{)}f(b_{T}^{\mathrm{cut}})$
$\displaystyle=-\frac{15}{128}f(b_{T}^{\mathrm{cut}})+\frac{7}{8}b_{T}^{\mathrm{cut}}f^{\prime}(b_{T}^{\mathrm{cut}})+(b_{T}^{\mathrm{cut}})^{2}f^{\prime\prime}(b_{T}^{\mathrm{cut}})\,.$
(2.21)

The first term reproduces the $S^{(1)}[f]$ given in eq. (2.9).

Figure -298: Left: Scaling of error terms of the truncated functionals for the
$k_{T}$ spectrum of the toy test function in eq. (2.2). The error terms
(computed as the difference to the exact result) exactly follow the expected
scaling $\sim k_{T}^{n-1/2}$. Right: Same, but instead applying the truncated
functionals to the _leading “OPE” contribution_ $f^{(0)}_{\rm toy}$ for the
test function. The difference to the full result for $n\geq 2$ now clearly
reveals the presence of quadratic terms $\sim b_{T}^{2}$ in the full spectrum
that were expanded away by making use of the “perturbative” prediction
$f^{(0)}_{\rm toy}$.

In the left panel of figure -298, we test the asymptotic expansion for a
simple, but physically relevant toy test function

$\displaystyle f_{\mathrm{toy}}(b_{T})$
$\displaystyle=\exp\bigl{[}-\lambda_{1}\ln^{2}(b_{T}Q)\big{]}\,\exp[-\lambda_{2}b_{T}^{2}]\,,$
$\displaystyle\lambda_{1}$
$\displaystyle=0.1\,,\qquad\lambda_{2}=1\,\mathrm{GeV}^{2}\,,\qquad
Q=100\,\mathrm{GeV}\,.$ (2.22)

Here the $\lambda_{1}$ term is a proxy for the perturbative Sudakov
exponential appearing in TMDs, and the $\lambda_{2}$ term is a proxy for a
nonperturbative contribution. The full spectrum $S[f_{\mathrm{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_{\mathrm{toy}}]$. We clearly see that compared to the
full spectrum $S[f_{\mathrm{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$.

###  2.3 Marrying the OPE to the truncated functionals

Consider a function $f(b_{T})$ which, apart from satisfying assumptions
$(1^{*})$ and $(2^{*})$ in eq. (2.2), can also be written in terms of an
asymptotic series expansion,

$\displaystyle
f(b_{T})=\sum_{j=0}^{m-1}b_{T}^{j}f^{(j)}(b_{T})+\mathcal{O}\bigl{[}(\Lambda_{\mathrm{QCD}}b_{T})^{m}\bigr{]}\,,$
(2.23)

where $f^{(j)}(b_{T})\sim\Lambda_{\mathrm{QCD}}^{j}$ 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
$\Lambda_{\mathrm{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_{T}^{\mathrm{OPE}}$, such that the limit $m\to\infty$ exists for all
$b_{T}<b_{T}^{\mathrm{OPE}}$, 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 $\sim\exp\bigl{(}-1/(\Lambda_{\mathrm{QCD}}b_{T})\bigr{)}$.
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 $\sim(\Lambda_{\mathrm{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_ $\sum_{j=0}^{m-1}b_{T}^{j}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 $\Lambda_{\mathrm{QCD}}b_{T}$, i.e.,
$f(b_{T})=f^{(0)}(b_{T})+\mathcal{O}[(\Lambda_{\mathrm{QCD}}b_{T})^{2}]$, as
is often the case. After the Fourier transform, we have corrections

$\displaystyle S[f](k_{T})$
$\displaystyle=S^{(n)}[f](k_{T})+\frac{f(b_{T}^{\mathrm{cut}})}{k_{T}}\mathcal{O}\bigl{[}(k_{T}b_{T}^{\mathrm{cut}})^{-n+\frac{1}{2}}\bigr{]}$
(2.24)
$\displaystyle=S^{(n)}[f^{(0)}]+\frac{1}{k_{T}}\mathcal{O}\Bigl{[}(k_{T}b_{T}^{\mathrm{cut}})^{-n+\frac{1}{2}},\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{k_{T}}\Bigr{)}^{2},\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{k_{T}}\Bigr{)}^{2}(k_{T}b_{T}^{\mathrm{cut}})^{-n+\frac{5}{2}}\Bigr{]}\,,$

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_{T}^{\mathrm{cut}})^{2}f^{(2)}(b_{T}^{\mathrm{cut}})$ in the surface term
and the error term of the functional, is suppressed more strongly than the
intrinsic OPE correction $(\Lambda_{\mathrm{QCD}}/k_{T})^{2}$ for functionals
with sufficiently high order $n>\frac{5}{2}$. In particular, there are no
corrections of
$\mathcal{O}[(\Lambda_{\mathrm{QCD}}b_{T}^{\mathrm{cut}})^{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 -298, where we test eq. (2.24) with
the toy function in eq. (2.2). We put only the leading perturbative term
$f_{\mathrm{toy}}^{(0)}$ into the truncated functionals, where

$\displaystyle f_{\mathrm{toy}}^{(0)}$
$\displaystyle=\exp[-\lambda_{1}\ln^{2}(b_{T}Q)]\,,$ $\displaystyle
f_{\mathrm{toy}}$
$\displaystyle=f_{\mathrm{toy}}^{(0)}\bigl{(}1-\lambda_{2}b_{T}^{2}+\mathcal{O}(b_{T}^{4})\bigr{)}\,.$
(2.25)

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
$-\lambda_{2}b_{T}^{2}$. 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.
$f(b_{T})=\sum_{j=0}^{m-1}b_{T}^{j}f^{(j)}(b_{T})+\mathcal{O}[(\Lambda_{\mathrm{QCD}}b_{T})^{m}]$.
We need to compute $S^{(n)}[b_{T}^{m}f^{(m)}]$ to find the leading momentum-
space error term from the truncation of the OPE. Note that for
$b_{T}^{m}f^{(m)}(b_{T})$ itself, the full functional does not exist, since
$b_{T}^{m}f^{(m)}(b_{T})\sim b_{T}^{m}$ with $m>1/2$ at large $b_{T}$.
However, we can construct a function $f^{(m)}_{*}(b_{T})$ that agrees with
$b_{T}^{m}f^{(m)}(b_{T})$ for all $b_{T}\leq b_{T}^{\mathrm{cut}}$ but also
satisfies assumptions $(1^{*})$ and $(2^{*})$ in eq. (2.2), e.g. by applying
an exponential $b^{*}$ prescription (see section 2.5) to
$b_{T}^{m}f^{(m)}(b_{T})$ and multiplying the result with a similarly
constructed exponential model. Then we have

$\displaystyle S^{(n)}[b_{T}^{m}f^{(m)}](k_{T})$
$\displaystyle=S^{(n)}[f^{(m)}_{*}](k_{T})$ (2.26)
$\displaystyle=\underbrace{S[f^{(m)}_{*}](k_{T})}+\underbrace{\frac{f^{(m)}_{*}(b_{T}^{\mathrm{cut}})}{k_{T}}\mathcal{O}\bigl{[}(k_{T}b_{T}^{\mathrm{cut}})^{-n+\frac{1}{2}}\bigr{]}}\,,$
$\displaystyle\mathcal{O}\Bigl{[}\frac{1}{k_{T}}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{k_{T}}\Bigr{)}^{m}\Bigr{]}\qquad\mathcal{O}\Bigl{[}\frac{1}{k_{T}}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{k_{T}}\Bigr{)}^{m}(k_{T}b_{T}^{\mathrm{cut}})^{-n+m+\frac{1}{2}}\Bigr{]}$

where in the last equality we have used that $S[f^{(m)}_{*}](k_{T})$ is
independent of $b_{T}^{\mathrm{cut}}$, and thus by dimensional analysis can
only yield the intrinsic OPE correction
$\mathcal{O}[(\Lambda_{\mathrm{QCD}}/k_{T})^{m}]$. We see that the second
term, which accounts for corrections in
$(\Lambda_{\mathrm{QCD}}b_{T}^{\mathrm{cut}})^{m}$, is suppressed compared to
the first when $n>m+\frac{1}{2}$ without having to assume that
$b_{T}^{\mathrm{cut}}\ll 1/\Lambda_{\mathrm{QCD}}$. Combined with eq. (2.20)
we have

$\displaystyle S[f](k_{T})$
$\displaystyle=\sum_{j=0}^{m-1}S^{(n)}[b_{T}^{j}f^{(j)}]+\frac{1}{k_{T}}\mathcal{O}\Bigl{[}(k_{T}b_{T}^{\mathrm{cut}})^{-n+\frac{1}{2}},\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{k_{T}}\Bigr{)}^{m},\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{k_{T}}\Bigr{)}^{m}(k_{T}b_{T}^{\mathrm{cut}})^{-n+m+\frac{1}{2}}\Bigr{]}\,.$
(2.27)

Note that each term $S^{(n)}[b_{T}^{j}f^{(j)}]$ is individually well defined
despite the growth $\sim b_{T}^{j}$ 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. (2.24), where power corrections of
$\mathcal{O}[(\Lambda_{\mathrm{QCD}}b_{T}^{\mathrm{cut}})^{m}]$ are absent and
the third class of power corrections is suppressed for
$n>m+\frac{1}{2}=\frac{5}{2}$.

Eq. (2.27) is the main result of this section. It allows us to compute a
momentum-space 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.

###  2.4 Extension to cumulative distributions

In physics applications, we are often interested in a slightly different, but
related integral functional,

$\displaystyle K[f](k_{T}^{\mathrm{cut}})$
$\displaystyle\equiv\int_{0}^{k_{T}^{\mathrm{cut}}}\\!\\!\mathrm{d}k_{T}\
S[f](k_{T})$ (2.28)
$\displaystyle=k_{T}^{\mathrm{cut}}\int_{0}^{\infty}\\!\mathrm{d}b_{T}\,J_{1}(b_{T}k_{T}^{\mathrm{cut}})\,f(b_{T})\,,$

which corresponds to the cumulative distribution over all momenta
$\lvert{\vec{k}}_{T}\rvert\leq k_{T}^{\mathrm{cut}}$ in an azimuthally
symmetric position-space function $f(b_{T})$, see eq. (1.1).

We may again introduce a perturbative cutoff $b_{T}^{\mathrm{cut}}$ and divide
the $K[f]$ into $K_{<}[f]$ and $K_{>}[f]$ analogous to eq. (2.1), and apply
the same procedure described in section 2.2 to derive a series of functionals
$K^{(n)}[f]$ that asymptote to $K[f]$ for large $k_{T}^{\mathrm{cut}}$,

$\displaystyle K^{(n)}[f](k_{T}^{\mathrm{cut}},b_{T}^{\mathrm{cut}})$
$\displaystyle\equiv
k_{T}^{\mathrm{cut}}\int_{0}^{b_{T}^{\mathrm{cut}}}\\!\mathrm{d}b_{T}J_{1}(k_{T}^{\mathrm{cut}}b_{T})f(b_{T})$
(2.29) $\displaystyle\quad+\sqrt{\frac{2}{\pi
k_{T}^{\mathrm{cut}}b_{T}^{\mathrm{cut}}}}\sum_{k=0}^{n-1}\frac{\cos\Bigl{(}k_{T}^{\mathrm{cut}}b_{T}^{\mathrm{cut}}-\frac{\pi}{4}-\frac{k\pi}{2}\Bigr{)}}{(k_{T}^{\mathrm{cut}}b_{T}^{\mathrm{cut}})^{k}}\,\tilde{c}_{k}\Bigl{(}b_{T}^{\mathrm{cut}}\frac{\mathrm{d}}{\mathrm{d}b_{T}^{\mathrm{cut}}}\Bigr{)}f(b_{T}^{\mathrm{cut}})\,.$

These truncated cumulative functionals $K^{(n)}[f]$ approximate $K[f]$ up to
an error term given by

$\displaystyle
K[f](k_{T}^{\mathrm{cut}})=K^{(n)}[f](k_{T}^{\mathrm{cut}},b_{T}^{\mathrm{cut}})+\mathcal{O}\bigl{[}(k_{T}^{\mathrm{cut}}b_{T}^{\mathrm{cut}})^{-n-\frac{1}{2}}\bigr{]}\,.$
(2.30)

The expressions above are analogous to eqs. (2.2) and (2.20) for the spectrum.
We give the first few coefficients needed to construct $K^{(n)}[f]$ with eq.
(2.29) for $n\leq 3$:

$\displaystyle\tilde{c}_{0}\Bigl{(}b_{T}^{\mathrm{cut}}\frac{\mathrm{d}}{\mathrm{d}b_{T}^{\mathrm{cut}}}\Bigr{)}f(b_{T}^{\mathrm{cut}})$
$\displaystyle=f(b_{T}^{\mathrm{cut}})\,,$
$\displaystyle\tilde{c}_{1}\Bigl{(}b_{T}^{\mathrm{cut}}\frac{\mathrm{d}}{\mathrm{d}b_{T}^{\mathrm{cut}}}\Bigr{)}f(b_{T}^{\mathrm{cut}})$
$\displaystyle=-\frac{1}{8}f(b_{T}^{\mathrm{cut}})+b_{T}^{\mathrm{cut}}f^{\prime}(b_{T}^{\mathrm{cut}})\,,$
$\displaystyle\tilde{c}_{2}\Bigl{(}b_{T}^{\mathrm{cut}}\frac{\mathrm{d}}{\mathrm{d}b_{T}^{\mathrm{cut}}}\Bigr{)}f(b_{T}^{\mathrm{cut}})$
$\displaystyle=-\frac{9}{128}f(b_{T}^{\mathrm{cut}})+\frac{5}{8}b_{T}^{\mathrm{cut}}f^{\prime}(b_{T}^{\mathrm{cut}})-(b_{T}^{\mathrm{cut}})^{2}f^{\prime\prime}(b_{T}^{\mathrm{cut}})\,.$
(2.31)

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

$\displaystyle K[f](k_{T}^{\mathrm{cut}})$
$\displaystyle=\sum_{j=0}^{m-1}K^{(n)}[b_{T}^{j}f^{(j)}]+\mathcal{O}\Bigl{[}(k_{T}^{\mathrm{cut}}b_{T}^{\mathrm{cut}})^{-n-\frac{1}{2}},\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{k_{T}^{\mathrm{cut}}}\Bigr{)}^{m},\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{k_{T}^{\mathrm{cut}}}\Bigr{)}^{m}(k_{T}^{\mathrm{cut}}b_{T}^{\mathrm{cut}})^{-n+m-\frac{1}{2}}\Bigr{]}\,,$
(2.32)

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.

Figure -297: Same as figure -298, but for the cumulative functionals
$K^{(n)}[f](k_{T}^{\mathrm{cut}})$.

The left panel of figure -297 shows that the cumulative functionals exhibit
the expected power-law scaling of error terms for the toy function in eq.
(2.2). In the right panel of figure -297 we again restrict the input of the
functionals to the leading term at small $b_{T}$,
$K^{(n)}[f_{\mathrm{toy}}^{(0)}]$, and compute the difference to the full
result $K[f_{\mathrm{toy}}]$. Starting from $n=2$, we find that the intrinsic
$\mathcal{O}(b_{T}^{2})$ correction results in a power-law behavior $\sim
k_{T}^{\mathrm{cut}}$ of this difference, as expected. Note that for the
cumulative distributions, $n>m-\frac{1}{2}$ rather than $m+\frac{1}{2}$ is
already sufficient to expose the quadratic behavior of the intrinsic OPE
correction.

###  2.5 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]:666We have taken $b_{\mathrm{min}}\to 0$ in the
original $b^{*}_{\mathrm{Pavia}}$ prescription, which as argued in ref. [26]
only amounts to a modification in the UV (beyond the validity of TMD
factorization). This allows for a more direct comparison to other models which
have $b^{*}(b_{T}\to 0)\to 0$.

$\displaystyle b_{\mathrm{CS}}^{*}(b_{T})$
$\displaystyle=\frac{b_{T}}{\sqrt{1+(b_{T}/b_{\mathrm{max}}^{2}})}=b_{T}\Bigl{(}1-\frac{b_{T}^{2}}{2\,b_{\mathrm{max}}^{2}}\Bigr{)}+\mathcal{O}\bigl{(}b_{T}^{5}\bigr{)}\,,$
$\displaystyle b_{\mathrm{Pavia}}^{*}(b_{T})$
$\displaystyle=b_{\mathrm{max}}\biggl{[}1-\exp\Bigl{(}-\frac{b_{T}^{4}}{b_{\mathrm{max}}^{4}}\Bigr{)}\biggr{]}^{\frac{1}{4}}=b_{T}\Bigl{(}1-\frac{b_{T}^{4}}{8\,b_{\mathrm{max}}^{4}}\Bigr{)}+\mathcal{O}\bigl{(}b_{T}^{9}\bigr{)}\,.$
(2.33)

When evaluating $f$ at $b^{*}(b_{T})$ and expanding for $b_{T}\ll
b_{\mathrm{max}}$, we have

$\displaystyle
f\bigl{(}b_{\mathrm{CS}}^{*}(b_{T})\bigr{)}=f(b_{T})-\frac{b_{T}^{3}}{2\,b_{\mathrm{max}}^{2}}f^{\prime}(b_{T})+\mathcal{O}\bigl{(}b_{T}^{5}\bigr{)}\,,$
$\displaystyle
f\bigl{(}b_{\mathrm{Pavia}}^{*}(b_{T})\bigr{)}=f(b_{T})-\frac{b_{T}^{5}}{8\,b_{\mathrm{max}}^{4}}f^{\prime}(b_{T})+\mathcal{O}\bigl{(}b_{T}^{9}\bigr{)}\,,$
(2.34)

where we note that the first subleading terms in the series differ, with
$b_{T}^{3}$ and $b_{T}^{5}$ scaling respectively.

We can formally assess the _momentum-space_ effect of the residual
$b_{\mathrm{max}}$ dependence by applying our truncated functionals at
sufficiently high orders to $f(b^{*}(b_{T}))$. For example, using $S^{(n)}$
with $n\geq 3$ for CS and $n\geq 5$ for Pavia, we can show using eq. (2.27)
that the leading correction terms in eq. (2.5) carries over to momentum space
with the scaling $\sim 1/(b_{\mathrm{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
-298.

Figure -296: Scaling of the ambiguity from the choice of $b_{\mathrm{max}}$
for different $b^{*}$ prescriptions.

To test the impact of different $b^{*}$ choices we consider $S[\sigma]$ where
$\sigma$ is the leading (double) logarithmic (LL) inclusive $q_{T}$
differential cross section for $pp\to Z/\gamma^{*}$ at the LHC, normalized to
the tree-level cross section. We use the LL solution of the running coupling
starting from $\alpha_{s}(m_{Z})=0.118$, which exhibits a Landau pole at
$\Lambda_{\mathrm{QCD}}^{\mathrm{LL}}\approx 100\,\mathrm{MeV}$. For numerical
stability, we keep the scale of the PDFs in the resummed prediction fixed at
the high scale $\mu_{f}=Q=m_{Z}$, which is an NLL effect. In figure -296, we
show the scaling of the difference between using
$b_{\mathrm{max}}^{A}=1.87\,\mathrm{GeV}^{-1}$ and
$b_{\mathrm{max}}^{B}=1.40\,\mathrm{GeV}^{-1}$ in momentum space for different
$b^{*}$ models. We also include the value of an individual $S[\sigma]$ (dashed
gray curve) using $b^{*}_{\mathrm{Pavia}}$ with
$b_{\mathrm{max}}=1.87\,\mathrm{GeV}^{-1}$ as reference. All curves are
normalized to the Born cross section $\sigma_{\mathrm{LO}}$. We consider the
behavior of $b^{*}_{\mathrm{CS}}$ (dotted black) and $b^{*}_{\mathrm{Pavia}}$
(dashed orange), as well as an exponential $b^{*}$ prescription (solid blue)
defined as

$b^{*}_{\mathrm{exp}}(b_{T})=\begin{cases}b_{T}\Bigl{[}1-\exp\Bigl{(}-\frac{1}{b_{T}/b_{\mathrm{max}}-a}\Bigr{)}\Bigr{]}\,,\quad&b_{T}\geq
a\,b_{\mathrm{max}},\\\ b_{T}\,,\quad&b_{T}<a\,b_{\mathrm{max}}\,.\end{cases}$
(2.35)

The benefit of the $b^{*}_{\mathrm{exp}}$ prescription is that it does not
introduce additional powers beyond those predicted by the OPE. In addition,
for $b_{T}\leq ab_{\mathrm{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^{*}_{\mathrm{exp}}$ is still monotonic. We observe that the scaling of the
$b_{\mathrm{max}}$ ambiguity of the CS and Pavia prescriptions both stabilize
to their corresponding power law ($q_{T}^{-2}$ and $q_{T}^{-4}$,
respectively), as expected from our analytic arguments below eq. (2.5). For
our newly introduced $b^{*}_{\mathrm{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}\to 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 $\mu_{\mathrm{OPE}}=b_{0}/b_{T}+\Delta$ with
$\Delta=2\,\mathrm{GeV}$ at which the OPE of the TMD PDF is performed in ref.
[38], and which we expect to translate to linear power corrections
$\mathcal{O}(\Delta/q_{T})$ in momentum space. In either of these cases, the
total $b_{T}$ derivative in eq. (2.5) is replaced by partial derivatives
$\partial/\partial\mu_{X}$ with respect to the relevant scale(s) $\mu_{X}$,
but the power-law scaling of the artifacts induced by the prescription still
follows from expanding the relevant profile scale function $\mu_{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}\mu_{X}$ may be large, so the $\mu_{X}$ dependence is not
guaranteed to decrease as the perturbative order increases.

###  2.6 Series stability near the Landau pole

Figure -295: Absolute error of the truncated functionals in the presence of a
Landau pole for the cumulative distribution of a single TMD PDF (left) and the
$Z/\gamma^{\ast}$ transverse-momentum spectrum at the LHC (right).

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_{T}^{\mathrm{pole}}$, we
generically expect the derivatives of the test function to grow factorially,
e.g.,

$\displaystyle f(b_{T})=\frac{1}{b_{T}-b_{T}^{\mathrm{pole}}}\,,\qquad
f^{(n)}(b_{T})=\frac{(-1)^{n}n!}{\bigl{(}b_{T}-b_{T}^{\mathrm{pole}}\bigr{)}^{n+1}}\,.$
(2.36)

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
$\Lambda_{\mathrm{QCD}}\approx 100\,\mathrm{MeV}$, but without applying any
$b^{*}$ prescription in this case. We continue to fix the PDF scales at the
high scale $\mu_{f}=\sqrt{\zeta}$ for the TMD PDF and at $\mu_{f}=Q$ for the
cross section, for numerical stability.

In figure -295 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 -298 and -297, 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.

##  3 Calculating the TMD PDFs and their cumulative integral

###  3.1 Evolution and operator product expansion of the TMD PDF

The dependence of the TMD PDF on the renormalization scale $\mu$ and the CS
scale $\zeta$ 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 $(\mu_{0},\zeta_{0})$ to $(\mu,\zeta)$,

$\displaystyle f_{i}(x,b_{T},\mu,\zeta)$
$\displaystyle=f_{i}(x,b_{T},\mu_{0},\zeta_{0})\,\exp\biggl{[}\int_{\mu_{0}}^{\mu}\\!\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\,\gamma^{i}_{\mu}(\mu^{\prime},\zeta_{0})\biggr{]}\,\exp\biggl{[}\frac{1}{2}\gamma^{i}_{\zeta}(b_{T},\mu)\ln\frac{\zeta}{\zeta_{0}}\biggr{]}$
$\displaystyle=f_{i}(x,b_{T},\mu,\zeta_{0})\,\exp\biggl{[}\frac{1}{2}\gamma^{i}_{\zeta}(b_{T},\mu)\ln\frac{\zeta}{\zeta_{0}}\biggr{]}\,.$
(3.1)

Here we have chosen the evolution path as
$(\mu_{0},\zeta_{0})\to(\mu,\zeta_{0})\to(\mu,\zeta)$. We can in particular
use eq. (3.1) to evolve the TMD PDF away from the canonical boundary condition
$\mu_{0}\,,\sqrt{\zeta_{0}}\sim 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 $\Lambda_{\mathrm{QCD}}b_{T}$. In explicit numerics we choose
the canonical boundary scales as

$\displaystyle\mu_{0}=\sqrt{\zeta_{0}}=\frac{b_{0}}{b_{T}}\,,$ (3.2)

where we include a conventional $\mathcal{O}(1)$ numerical coefficient
$b_{0}=2\exp^{-\gamma_{E}}\approx 1.12292$ with $\gamma_{E}$ the Euler-
Mascheroni constant. We will also refer to this boundary condition simply as
the “intrinsic” TMD PDF.777Other 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 $\sqrt{\zeta_{0}}$ is held fixed at a low scale
$Q_{0}=1\,\mathrm{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 $\zeta_{0}$
amount to moving nonperturbative contributions proportional to the rapidity
anomalous dimension in and out of the intrinsic TMD PDF.

We are interested in values of $1/b_{T},\mu_{0},\mu\gg\Lambda_{\mathrm{QCD}}$
where the TMD PDF and the rapidity anomalous dimension $\gamma^{i}_{\zeta}$
obey an operator product expansion [52, 18, 59, 60, 61],

$\displaystyle f_{i}(x,b_{T},\mu,\zeta)$
$\displaystyle=f_{i}^{(0)}(x,b_{T},\mu,\zeta)+\sum_{j=1}^{m-1}b_{T}^{j}f_{i}^{(j)}(x,b_{T},\mu,\zeta)+\mathcal{O}\bigl{[}(\Lambda_{\mathrm{QCD}}b_{T})^{m}\bigr{]}\,,$
$\displaystyle\gamma^{i}_{\zeta}(b_{T},\mu)$
$\displaystyle=\gamma^{(0)}_{\zeta,i}(b_{T},\mu)+\sum_{j=1}^{m-1}b_{T}^{j}\gamma^{(j)}_{\zeta,i}(b_{T})+\mathcal{O}\bigl{[}(\Lambda_{\mathrm{QCD}}b_{T})^{m}\bigr{]}\,,$
(3.3)

with $f_{i}^{(j)},\gamma^{(j)}_{\zeta,i}\sim\Lambda_{\mathrm{QCD}}^{j}$. Most
often the terms with odd $j$ vanish. The $\mu$ 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 $\mu$, which does not need
to be equal to $\mu_{0}$. On the other hand, the $\zeta$ evolution mixes
powers in the $b_{T}$ expansion. Specifically, for the first nonvanishing
orders in a TMD PDF with $f_{i}^{(0)}\neq 0$, $f_{i}^{(1)}=0$,

$\displaystyle\zeta\frac{\mathrm{d}}{\mathrm{d}\zeta}f^{(0)}_{i}(x,b_{T},\mu,\zeta)$
$\displaystyle=\frac{1}{2}\gamma^{(0)}_{\zeta,i}(b_{T},\mu)\,f^{(0)}_{i}(x,b_{T},\mu,\zeta)\,,$
(3.4)
$\displaystyle\zeta\frac{\mathrm{d}}{\mathrm{d}\zeta}f^{(2)}_{i}(x,b_{T},\mu,\zeta)$
$\displaystyle=\frac{1}{2}\Bigl{[}\gamma^{(0)}_{\zeta,i}(b_{T},\mu)\,f^{(2)}_{i}(x,b_{T},\mu,\zeta)+\gamma^{(2)}_{\zeta,i}(b_{T})\,f^{(0)}_{i}(x,b_{T},\mu,\zeta)\Bigr{]}\,.$

The structure of $\gamma_{\zeta}^{i}$ at higher orders is relatively well
understood. The leading perturbative term $\gamma_{\zeta,i}^{(0)}$ is known to
three loops [62, 63, 64], and also contributes the overall additive $\mu$
renormalization of $\gamma_{\zeta}^{i}$ in eq. (1.8). Recently, ref. [61]
demonstrated that the $\mathcal{O}(b_{T}^{2})$ 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],
$\gamma_{\zeta,i}^{(2)}(b_{T})=\gamma_{\zeta,i}^{(2)}\sim\Lambda_{\mathrm{QCD}}^{2}$
is a single fixed number with units of $\,\mathrm{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}\to 0$, and higher-power matrix
elements that provide the power suppression by $\Lambda_{\mathrm{QCD}}^{j}$.
The leading term for the unpolarized quark or gluon TMD PDF takes the form

$\displaystyle f^{(0)}_{i}(x,b_{T},\mu,\zeta)$
$\displaystyle=\sum_{j}\int_{x}^{1}\\!\frac{\mathrm{d}z}{z}\,C_{ij}(z,\,b_{T},\,\mu,\,\zeta)\,f_{j}\Bigl{(}\frac{x}{z},\mu\Bigr{)}$
$\displaystyle=f_{i}(x,\mu)+\mathcal{O}[\alpha_{s}(\mu)]\,,$ (3.5)

where $f_{j}(x/z,\mu)$ is the collinear PDF,
$C_{ij}(z,b_{T},\mu,\zeta)=\delta_{ij}\delta(1-z)+\mathcal{O}[\alpha_{s}(\mu)]$
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,\mu)$ was determined to all powers in $b_{T}$ and at tree level in
$\alpha_{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 $\mathcal{O}(b_{T}^{j})$ in this case; nevertheless, they
only constitute a subset of the full OPE at any given order in
$b_{T}$.888Interestingly, 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
$\displaystyle K\bigl{[}f_{i\,\text{tree-
level}}^{\text{ref.~{}\cite[cite]{[\@@bibref{Number}{Moos:2020wvd}{}{}]}}}(x,b_{T})\bigr{]}(k_{T}^{\mathrm{cut}})=f_{i}(x,\mu)-\frac{1}{[1+(k_{T}^{\mathrm{cut}}/xM)^{2}]}f_{i}\Bigl{(}x+\frac{(k_{T}^{\mathrm{cut}})^{2}}{xM^{2}},\mu\Bigr{)}\,.$
(3.6) In particular, the $k_{T}$ dependence of the TMD PDF in this
approximation is cut off at $k_{T}^{2}\leq 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}\geq 0$ with $P^{\mu}$ the hadron and $k^{\mu}$
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. For the purposes of this paper, we will consider the
most general form that the quadratic correction can have based on the RG
properties,

$\displaystyle
b_{T}^{2}f^{(2)}_{i}(x,b_{T},\mu,\zeta)=b_{T}^{2}\,\Lambda_{i}^{(2)}(x,b_{T},\zeta)\,f^{(0)}_{i}(x,b_{T},\mu,\zeta)\Bigl{[}1+\mathcal{O}\bigl{[}\alpha_{s}(\mu)\bigr{]}\Bigr{]}\,,$
(3.7)

where $\Lambda_{i}^{(2)}\sim\Lambda_{\mathrm{QCD}}^{2}$ has units of
$\,\mathrm{GeV}^{2}$. To see this, recall that the $\mu$ 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 $\Lambda_{i}^{(2)}$ is renormalized as

$\displaystyle\zeta\frac{\mathrm{d}}{\mathrm{d}\zeta}\Lambda_{i}^{(2)}(x,b_{T},\zeta)=\frac{1}{2}\gamma_{\zeta,i}^{(2)}(b_{T})\,.$
(3.8)

Evaluating $\Lambda_{i}^{(2)}$ at our reference boundary scale
$\sqrt{\zeta_{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 leading-
order Wilson coefficients and ignore the running of the higher-twist matrix
elements and soft condensates from $\Lambda_{\mathrm{QCD}}$ to $1/b_{T}$. (The
running in general would dress the leading scaling $\sim b_{T}^{2}$ with an
anomalous exponent, and similarly for Wilson coefficients that start at
$\mathcal{O}[\alpha_{s}(1/b_{T})]$.) This leaves us with a function
$\Lambda_{i}^{(2)}(x,b_{T},\zeta_{0})=\Lambda_{i}^{(2)}(x)$ of a single
variable for each flavor.

Combining both sources of $\mathcal{O}(b_{T}^{2})$ nonperturbative corrections
and fixing $\zeta_{0}$, we have

$\displaystyle
f_{i}(x,b_{T},\mu,\zeta)=f_{i}^{(0)}(x,b_{T},\mu,\zeta)\Bigl{\\{}1+b_{T}^{2}\Bigl{(}\Lambda_{i}^{(2)}(x)+\frac{1}{2}\gamma^{(2)}_{\zeta,i}L_{\zeta}\Bigr{)}\Bigr{\\}}+\mathcal{O}\bigl{[}(\Lambda_{\mathrm{QCD}}b_{T})^{4}\bigr{]}\,,$
(3.9)

where we defined $L_{\zeta}\equiv\ln(b_{T}^{2}\zeta/b_{0}^{2})$ as a
shorthand. We also absorbed the leading-power rapidity evolution factor
$\exp\bigl{[}\frac{1}{2}\gamma^{(0)}_{\zeta,i}L_{\zeta}\bigr{]}$ into
$f_{i}^{(0)}$, which evolves it back to the overall $\zeta$. Importantly, the
impact of $\gamma^{(2)}_{\zeta,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
$\gamma_{\zeta,i}^{(2)}$ out of the Fourier transform. This setup will allow
us to study the linear impact of the OPE coefficients $\Lambda_{i}^{(2)}$ and
$\gamma^{(2)}_{\zeta,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.

###  3.2 Approximating the cumulative TMD PDF with truncated functionals

The cumulative integral of an azimuthally symmetric TMD PDF over
$\lvert{\vec{k}}_{T}\rvert\leq k_{T}^{\mathrm{cut}}$,

$\displaystyle\int_{\lvert{\vec{k}}_{T}\rvert\leq
k_{T}^{\mathrm{cut}}}\\!\mathrm{d}^{2}{\vec{k}}_{T}\,f_{i}(x,{\vec{k}}_{T},\mu,\zeta)$
$\displaystyle=K[f_{i}(x,b_{T},\mu,\zeta)](k_{T}^{\mathrm{cut}})\,,$ (3.10)

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_{T}^{\mathrm{cut}}\gg\Lambda_{\mathrm{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.1), we can calculate the resummed leading OPE term
$f_{i}^{(0)}$ 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_{i}^{(0)}$ are given in section 3.3. We evaluate
$f_{i}^{(0)}$ at N3LL order, which requires the two-loop matching coefficients
for the TMD PDF [68, 69, 62], 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
$\alpha_{s}(m_{Z})=0.118$.

Applying our truncated functionals to the leading OPE term,
$K^{(n)}[f_{i}^{(0)}]$, 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_{T}^{\mathrm{cut}}=b_{0}/(1.2\,\mathrm{GeV})$ such that all scales in the
resummed TMD PDF stay above $1.2\,\mathrm{GeV}$, i.e., well away from the
Landau pole. We use standard adaptive integration 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_{T}^{\mathrm{cut}}$,
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}<\infty$ using a double-exponential
method for oscillatory integrals [79, 80, 81] as implemented in SCETlib [66].

Figure -294: 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_{T}^{\mathrm{cut}}$ (left) and of the
momentum fraction $x$ (right). The blue band shows the uncertainty from
varying the position-space integral cutoff $b_{T}^{\mathrm{cut}}$, and the
yellow band shows the uncertainty from varying the quadratic OPE coefficients
$\gamma_{\zeta,q}^{(2)}$ and $\Lambda_{q}^{(2)}$. 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 $K^{(3)}[f_{i}^{(0)}]$, which is purely
perturbative. In both figures, the evolution scales are set to be
$\mu=\sqrt{\zeta}=k_{T}^{\mathrm{cut}}$, and the TMD PDF is for a down quark.
See figure -288 in appendix B for the results for other light quark flavors
including $u,\bar{d},\bar{u}$, and $s$.

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 -294, we show the normalized deviation of $K^{(3)}[f_{i}^{(0)}]$
from the collinear PDF as a function of the transverse-momentum cutoff
$k_{T}^{\mathrm{cut}}$ while keeping $x=0.01$ fixed. In the right panel of
figure -294, we show the same deviation from the collinear PDF as a function
of the momentum fraction $x$, while keeping
$k_{T}^{\mathrm{cut}}=10\,\mathrm{GeV}$ fixed. In both cases we set
$\mu=\sqrt{\zeta}=k_{T}^{\mathrm{cut}}$. 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_{T}^{\mathrm{cut}}\geq 5\,\mathrm{GeV}$, and for $x=0.01$
is essentially equal to the PDF for $k_{T}^{\mathrm{cut}}\geq
10\,\mathrm{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 $\Delta_{\mathrm{cut}}$, shown in blue in
figure -294, is estimated by varying $b_{T}^{\mathrm{cut}}$ from
$b_{0}/(1.5\,\mathrm{GeV})$ to $b_{0}/(1\,\mathrm{GeV})$. Since the residual
$b_{T}^{\mathrm{cut}}$ 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 $\lesssim b_{T}^{\mathrm{cut}}$, 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 $\gamma_{\zeta,i}^{(2)}$ and $\Lambda_{i}^{(2)}$ 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,

$\displaystyle
K^{(3)}[f_{i}](k_{T}^{\mathrm{cut}})=K^{(3)}[f_{i}^{(0)}](k_{T}^{\mathrm{cut}})+\Lambda_{i}^{(2)}K^{(3)}[b_{T}^{2}f_{i}^{(0)}]+\frac{1}{2}\gamma^{(2)}_{\zeta,i}K^{(3)}[b_{T}^{2}L_{\zeta}f_{i}^{(0)}]+\mathcal{O}\Bigl{[}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{k_{T}^{\mathrm{cut}}}\Bigr{)}^{4}\Bigr{]}.$
(3.11)

The impact of $\gamma_{\zeta,i}^{(2)}$ and $\Lambda_{i}^{(2)}$ 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
$\Lambda_{q}^{(2)}(x)\in[-1\,\mathrm{GeV}^{2},1\,\mathrm{GeV}^{2}]$
independent of $x$, which yields an uncertainty $\Delta_{\mathrm{NP}}^{f}$. We
stress again that this assumption is straightforward to relax point by point
in the right panel of figure -294 by a trivial rescaling. To estimate the
uncertainty due to higher OPE corrections to the CS kernel, we adopt a very
conservative range of
$\gamma_{\zeta,q}^{(2)}\in[-(0.5\,\mathrm{GeV})^{2},(0.5\,\mathrm{GeV})^{2}]$,
yielding the second OPE-related uncertainty component
$\Delta_{\mathrm{NP}}^{\zeta}$; recent model-based global fits [38, 39] and
the condensate-based estimate from ref. [61] all point to substantially lower
values of $\lvert\gamma_{\zeta,q}^{(2)}\rvert\lesssim
0.05\,\mathrm{GeV}^{2}\approx(0.22\,\mathrm{GeV})^{2}$. For readability, both
nonperturbative uncertainty components are added in quadrature in figure -294,
$\Delta_{\mathrm{NP}}^{\zeta}\oplus\Delta_{\mathrm{NP}}^{f}\equiv\bigl{[}(\Delta_{\mathrm{NP}}^{\zeta})^{2}+(\Delta_{\mathrm{NP}}^{f})^{2}\bigr{]}^{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_{T}^{\mathrm{cut}}$ 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_{T}^{\mathrm{cut}}$ increases, as expected
from their respective power counting.

Finally we also include in figure -294 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 N3LL
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}\approx 2\,\mathrm{GeV}$, and up to $10\%$ at $b_{0}/b_{T}\approx
1\,\mathrm{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 N3LL
Pavia fit as provided in ref. [84], and combined it with our implementation of
$f_{i}^{(0)}$ in SCETlib evaluated with the $b^{*}_{\mathrm{Pavia}}$
prescription of ref. [26], with $b_{\mathrm{min}}\to 0$ as in eq. (2.5). 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_{\mathrm{min}}=0$.999The 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 $\zeta$ prescription used in arTeMiDe.
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 -288 in appendix B, we show the same plots as in figure -294 for
other light quark flavors ($u,\bar{d},\bar{u}$, and $s$). We find that the
results are very similar to the discussion for the down quark in this section.

###  3.3 Perturbative uncertainty

Since we rely on the accuracy of the perturbative input $f_{i}^{(0)}$, 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 $\zeta$, which may be parametrically separate from
$\zeta_{0}\sim 1/b_{T}^{2}$, is a multi-scale problem. One way of treating it
in resummed perturbation theory is to evaluate it at
$\mu_{0}\sim\sqrt{\zeta_{0}}\sim 1/b_{T}$, where it is free of large
logarithms, and use eq. (3.1) to evolve it back to the overall $(\mu,\zeta)$.
The uncertainty estimate in that case could be estimated by varying the
boundary scales $(\mu_{0},\zeta_{0})$.

Alternatively, one can consider a factorization of the TMD PDF at generic
$\zeta$ into purely collinear and purely soft matrix elements,

$f_{i}(x,b_{T},\mu,\zeta)=B_{i}(x,b_{T},\mu,\nu/\sqrt{\zeta})\sqrt{S^{i}(b_{T},\mu,\nu)}\,,$
(3.12)

which individually are single-scale problems and feature an additional
rapidity renormalization scale $\nu$. In the SCET literature, these matrix
elements are known as beam and soft functions.101010We 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].
In this setup, the running in $\nu$ from the canonical beam scale
$\nu_{B}\sim\sqrt{\zeta}$ to the soft scale $\nu_{S}\sim 1/b_{T}$ resums the
large logarithms of $\zeta b_{T}^{2}$. 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, once the possible variations of the boundary
scales are considered. In particular, the factorized form in eq. (3.12) allows
for varying all of $\mu_{B,S}$ and $\nu_{B,S}$ individually, providing a more
robust estimate of the perturbative uncertainty that makes use of the full
available information about the factorization structure. (In practice, to
estimate the perturbative uncertainty on the leading OPE term $f_{i}^{(0)}$ we
of course only need to evaluate eq. (3.12) at leading power in
$\Lambda_{\mathrm{QCD}}b_{T}$.)

The procedure of ref. [67] is to consider the variations111111Note that in
ref. [67], hybrid profile scales were used to turn off the resummation as
$q_{T}\to Q$ to implement the matching to power corrections of
$\mathcal{O}(q_{T}^{2}/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 $x_{1,2,3}\to\infty$ in the profile scales
in ref. [67].

$\displaystyle\mu_{B}=2^{v_{\nu_{S}}}\,\frac{b_{0}}{b_{T}}\,,\qquad\nu_{B}=2^{v_{\nu_{S}}}\,\sqrt{\zeta}\,,\qquad\mu_{S}=2^{v_{\mu_{S}}}\,\frac{b_{0}}{b_{T}}\,,\qquad\nu_{S}=2^{v_{\nu_{S}}}\,\frac{b_{0}}{b_{T}}\,,$
(3.13)

where the four variation exponents can be $v_{i}=\\{-1,0,1\\}$. The central
scale choices are given by
$(v_{\mu_{S}},v_{\nu_{S}},v_{\mu_{B}},v_{\nu_{B}})=(0,0,0,0)$. A priori there
are 80 possible different combined variations of the $v_{i}$. Since the
arguments of the resummed logarithms are ratios of scales, e.g.
$\nu_{B}/\nu_{S}$, some combinations of the $v_{i}$’s lead to variations of
these arguments by more than the conventional factor $2$ and are therefore
excluded, leaving us with 36 different combined scale variations. We calculate
$K^{(3)}[f_{i}^{(0)}]|_{v_{i}}$ for each of them and then take the maximum
envelope as our estimate of the perturbative resummation uncertainty
$\Delta_{\mathrm{res}}$.

Figure -293: Similar to figure -294, but for the perturbative uncertainties
from resummation orders, which are estimated from varying the scales in the
boundary condition of the TMD PDF evolution. When increasing the resummation
order, the uncertainty bands converge. The central values at all orders are
within a few percent from the collinear PDF. See figure -287 in appendix B for
the results for other light quark flavors including $u,\bar{d},\bar{u}$, and
$s$.

In figure -293 we show $\Delta_{\mathrm{res}}$ for the cumulative TMD PDF
calculated at NLL, NNLL, and N3LL as a function of $k_{T}^{\mathrm{cut}}$ 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 $\Delta_{\mathrm{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_{T}^{\mathrm{cut}}\gtrsim 10\,\mathrm{GeV}$, and the cumulative TMD PDF is
compatible with the collinear PDF well within the perturbative uncertainty for
all $k_{T}^{\mathrm{cut}}$ and $x$. Notice that the perturbative uncertainty
in figure -293 is much larger than the nonperturbative uncertainties shown in
figure -294.

Figure -287 in appendix B shows the same plots as in figure -293 for other
light quark flavors. The results are once again very similar to the discussion
for the down quark in this section.

###  3.4 Impact of evolution effects

In the discussions above, we have chosen
$\mu=\sqrt{\zeta}=k_{T}^{\mathrm{cut}}$, which is the natural scale choice
such that the left-hand side of eq. (3.10) only depends on one scale
$k_{T}^{\mathrm{cut}}$. 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 -292, we show the normalized deviation of the cumulative TMD PDF at
$k_{T}^{\mathrm{cut}}=10\,\mathrm{GeV}$ from the collinear PDF as a function
of both the renormalization scale $\mu$ and the CS scale $\zeta$. We find
that, in the vicinity of the natural scale choice
$\mu=\sqrt{\zeta}=k_{T}^{\mathrm{cut}}$, our conclusion is robust against
evolution effects and the deviation remains within percent level. In
particular, along the line of fixed $\mu=k_{T}^{\mathrm{cut}}$, we find that
the $\zeta$ evolution is negligible. This has the intriguing physical
interpretation that the $\zeta$ evolution of the momentum-space 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 $\mu=k_{T}^{\mathrm{cut}}$, 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
$(\mu_{a},\zeta_{a})=(k_{T}^{\mathrm{cut}},k_{T}^{\mathrm{cut}})\to(\mu_{a},\zeta_{b})\to(\mu_{b},\zeta_{b})$.
Since the $\zeta$ evolution at fixed $\mu_{a}=k_{T}^{\mathrm{cut}}$ 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
$\mu$ evolution factor resums large double logarithms of the form

$\displaystyle\exp\biggl{[}\int_{\mu_{a}}^{\mu_{b}}\\!\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\,\gamma^{i}_{\mu}(\mu^{\prime},\zeta_{b})\biggr{]}=1+\frac{\alpha_{s}(\mu_{a})}{2\pi}\Bigl{[}\frac{\Gamma_{0}^{i}}{4}\ln\frac{\mu_{b}^{2}}{\mu_{a}^{2}}\ln\frac{\mu_{a}^{2}}{\zeta_{b}}+\cdots\Bigr{]}+\mathcal{O}(\alpha_{s}^{2})\,,$
(3.14)

where the ellipses indicate single logarithms and
$\Gamma_{\mathrm{cusp}}[\alpha_{s}]=\frac{\alpha_{s}}{4\pi}\Gamma_{0}^{i}+\mathcal{O}(\alpha_{s}^{2})$
is the one-loop coefficient of the cusp anomalous dimension. Thus the effect
of the $\mu$ evolution is stronger for a larger ratio of
$\sqrt{\zeta_{b}}/\mu_{a}$ at the starting point, which is responsible for the
saddle point structure observed in figure -292. 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 $\mu,\sqrt{\zeta}\sim k_{T}^{\mathrm{cut}}$ is not
in contradiction with the renormalization, and is in fact observed in figure
-292 to also extend to parametrically separate values of $\zeta$ at fixed
$\mu=k_{T}^{\mathrm{cut}}$.

We conclude that, even when we take the evolution effects in the vicinity of
$\mu=k_{T}^{\mathrm{cut}}$ into account, our intuitive result that the
transverse-momentum integral of TMD PDF is within few percent of the
corresponding collinear PDF remains robust.

Figure -292: The normalized deviation of the transverse-momentum integral of
the TMD PDF from the collinear PDF as a function of the evolution scales $\mu$
and $\zeta$. The green star marks $\mu=\sqrt{\zeta}=k_{T}^{\mathrm{cut}}$. The
countour lines indicate increments of five percent points. There are only
percent-level deviations in the vicinity of
$\mu=\sqrt{\zeta}=k_{T}^{\mathrm{cut}}$. For $\mu=k_{T}^{\mathrm{cut}}$, the
evolution of $\zeta$ is negligible and only a shape effect. There are sizable
corrections from evolution away from the central regions. We approximate the
transverse-momentum integral of the TMD PDF with $K^{(3)}[f_{i}^{(0)}]$, which
is purely perturbative. In this figure, we use
$k_{T}^{\mathrm{cut}}=10\,\mathrm{GeV}$ and $x=0.01$, and the TMD PDF is for a
down quark. See figure -286 in appendix B for the results for other light
quark flavors including $u,\bar{d},\bar{u}$, and $s$.

Figure -291: The normalized deviation of the transverse-momentum integral of
the TMD PDF from the collinear PDF as a function of the evolution scale
$\sqrt{\zeta}$. We also include the SV and Pavia global fits in the left
panel. We approximate the transverse-momentum integral of the TMD PDF with
$K^{(3)}[f_{i}^{(0)}]$, which is purely perturbative. In both figures, we use
fixed $\mu=k_{T}^{\mathrm{cut}}=10\,\mathrm{GeV}$ and $x=0.01$, and the TMD
PDF is for a down quark. See figure -285 in appendix B for the results for
other light quark flavors including $u,\bar{d},\bar{u}$, and $s$.

It is also interesting to look at the $\zeta$ evolution alone. In figure -291,
we show the deviation of the transverse-momentum integral of a TMD PDF from
the collinear PDF as a function of $\sqrt{\zeta}$, while keeping $x=0.01$ and
$\mu=k_{T}^{\mathrm{cut}}=10\,\mathrm{GeV}$ fixed. In the left (right) panel
of figure -291, we consider nonperturbative (perturbative) uncertainties
$\Delta_{\mathrm{cut}}$, $\Delta_{\mathrm{NP}}^{\zeta}$, and
$\Delta_{\mathrm{NP}}^{\zeta}\oplus\Delta_{\mathrm{NP}}^{f}$
($\Delta_{\mathrm{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 $\sqrt{\zeta}\sim 1\,\mathrm{TeV}$. The
$\Delta_{\mathrm{cut}}$ uncertainty (light blue in the left panel) likewise
reduces at large $\sqrt{\zeta}$, which is expected as the Sudakov suppression
of the long-distance region increases, while the
$\Delta_{\mathrm{NP}}^{\zeta}$ (light green) component comes to dominate the
total OPE uncertainty (yellow) both at very small, $\zeta\to 1\,\mathrm{GeV}$,
and very large $\zeta\to 1\,\mathrm{TeV}$, i.e., far away from the natural
choice $\zeta\sim k_{T}^{\mathrm{cut}}$.

Figure -286 and figure -285 in appendix B show the same plots as in figure
-292 and figure -291 for other light quark flavors. The results are very
similar to the discussion for the down 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.

##  4 Model-independent constraints on nonperturbative physics from data

###  4.1 Outline of the method

Figure -290: Impact plots of quadratic “OPE” parameters in the toy function in
eq. (2.3) interpreted as nonperturbative uncertainties (left) or fits to data
(right). We consider the impact of the “nonperturbative” parameter
$\lambda_{2}$ on the toy model in eq. (2.2) with the expansion in eq. (2.3).
To produce “data” (solid lines on the right) for coefficients $\lambda_{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 $\Delta_{\mathrm{cut}}$, which
for clarity is only shown for $\lambda_{2}=0$ in the right panel.

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 $\sigma=f_{\mathrm{toy}}$ in eq.
(2.3). In the left panel of figure -290 we show the transverse-momentum
spectrum $S[f_{\mathrm{toy}}]$, where we vary the quadratic coefficient
$\lambda_{2}$ within two different ranges of values, as indicated by the
yellow bands. We normalize the impact of $\lambda_{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 $\lambda_{2}$ is slightly larger than expected from power counting, e.g.,
it amounts to $\sim 2\%$ at $k_{T}=10\,\mathrm{GeV}$ for
$\lambda_{2}=1\,\mathrm{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.3),

$\displaystyle
S^{(3)}[f_{\mathrm{toy}}^{(0)}+f_{\mathrm{toy}}^{(2)}]=S^{(3)}[f_{\mathrm{toy}}^{(0)}]-\lambda_{2}S^{(3)}[b_{T}^{2}f_{\mathrm{toy}}^{(0)}]\,,$
(4.1)

where $\lambda_{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 -290, where the dashed lines are the combined fit
template for signal and background for different values of $\lambda_{2}$. To
generate illustrative “data” (solid lines) for $\lambda_{2}<0$, where the
original toy function in eq. (2.2) has a divergent Fourier transform, we
modify it as

$\displaystyle
f_{\mathrm{toy}}(b_{T})=f^{(0)}_{\mathrm{toy}}(b_{T})\,\begin{cases}\exp(-\lambda_{2}b_{T}^{2})\,,\quad&\lambda_{2}>0\,,\\\
(1+2\lvert\lambda_{2}\rvert b_{T}^{2})\exp(-\lvert\lambda_{2}\rvert
b_{T}^{2})\,,\quad&\lambda_{2}<0\,,\end{cases}$ (4.2)

where $f_{\mathrm{toy}}^{(0)}=\exp[-\lambda_{1}\ln^{2}(b_{T}Q)]$ is unchanged.
This ensures that the full integral $S[f_{\mathrm{toy}}]$ exists for all
$\lambda_{2}$ while maintaining $f^{(2)}(b_{T})=-\lambda_{2}b_{T}^{2}$ for
either sign. We see that the fit template predicted by the truncated
functionals closely follows the mock data in figure -290 for all values of
$k_{T}$ and $\lambda_{2}$, with the exception of
$\lambda_{2}=-(1\,\mathrm{GeV})^{2}$ where our modification in eq. (4.2)
results in a strong modification of $f_{\mathrm{toy}}$ at large distances,
leading to a larger oscillatory error term.

###  4.2 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\to
Z/\gamma^{*}\to\ell^{+}\ell^{-}$, which are some of the most precise data
taken at the LHC to date [42, 43]. At leading power in $q_{T}\ll Q$, but to
all orders in $\Lambda_{\mathrm{QCD}}/q_{T}$, the fiducial cross section
differential in the dilepton transverse momentum reads

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}q_{T}}=\int\\!\mathrm{d}\Phi_{B}\,A_{B}(\Phi_{B})\,S\bigl{[}\sigma(Q,Y,b_{T})\bigr{]}(q_{T})\,\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{q_{T}}{Q}\Bigr{)}\Bigr{]}\,,$
(4.3)

where $\Phi_{B}=\\{Q,Y,\cos\theta\\}$ 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 $\theta$ 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}(\Phi_{B})$. For a
typical fiducial volume at the LHC,

$\displaystyle\mathrm{d}\Phi_{B}$
$\displaystyle=\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}\cos\theta\,\frac{3}{8}\bigl{(}1+\cos^{2}\theta\bigr{)}\,,$
$\displaystyle A_{B}(\Phi_{B})$
$\displaystyle=\Theta\bigl{(}Q_{\mathrm{min}}\leq Q\leq
Q_{\mathrm{max}}\bigr{)}\,\Theta\bigl{(}p_{T}^{1,2}\geq
p_{T}^{\mathrm{cut}}\bigr{)}\,\Theta\bigl{(}\lvert\eta_{1,2}\rvert\leq\eta_{\mathrm{cut}}\bigr{)}\,,$
(4.4)

where $Q_{\mathrm{min},\mathrm{max}}$ define an invariant-mass bin (e.g.
around the $Z$ resonance) and the lepton transverse momenta
$p_{T}^{1}=p_{T}^{2}=\sin\theta\,Q/2$ and pseudorapidities
$\eta_{1,2}=Y\mp\ln\bigl{(}\tan\theta/2\bigr{)}$ 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[\dots](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],

$\displaystyle\sigma(Q,Y,b_{T})=H(Q^{2},\mu)\sum_{i,j}\sigma^{B}_{ij}(Q)\,f_{i}(x_{a},b_{T},\mu,\zeta_{a})\,f_{j}(x_{b},b_{T},\mu,\zeta_{b})$
(4.5)

where $\sigma^{B}_{ij}(Q)$ is the Born cross section for the hard production
process $ij\to Z/\gamma^{*}$ that includes the $Q$-dependent $Z$ lineshape,
$H(Q^{2},\mu)=1+\mathcal{O}(\alpha_{s})$ is the so-called hard function that
encodes virtual corrections to the hard process,121212Here we have neglected
singlet contributions to the quark axial and vector form factors starting at
$\mathcal{O}(\alpha_{s}^{2})$ and $\mathcal{O}(\alpha_{s}^{3})$, 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 $\overline{\Lambda}^{(2)}$ coefficient that we outline below. 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

$\displaystyle x_{a,b}=\frac{Q}{E_{\mathrm{cm}}}\,e^{\pm
Y}\,,\qquad\zeta_{a,b}=Q^{2}$ (4.6)

with $E_{\mathrm{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
$\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2})$,

$\displaystyle\sigma(Q,Y,b_{T})$
$\displaystyle=H(Q^{2},\mu)\sum_{i,j}\sigma^{B}_{ij}(Q)\,f_{i}^{(0)}(x_{a},b_{T},\mu,\zeta_{a})\,f_{j}^{(0)}(x_{b},b_{T},\mu,\zeta_{b})$
$\displaystyle\quad\times\Bigl{\\{}1+b_{T}^{2}\Bigl{(}\Lambda_{i}^{(2)}(x_{a})+\Lambda_{j}^{(2)}(x_{b})+\gamma^{(2)}_{\zeta,q}L_{Q^{2}}\Bigr{)}\Bigr{\\}}+\mathcal{O}\bigl{[}(\Lambda_{\mathrm{QCD}}b_{T})^{4}\bigr{]}\,,$
(4.7)

where $L_{Q^{2}}=\ln(b_{T}^{2}Q^{2}/b_{0}^{2})$.

If we try to use eq. (4.2) 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 $\Lambda_{i,j}^{(2)}(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 $\overline{\Lambda}^{(2)}$ per fiducial volume that can
be determined by a two-parameter fit with $\gamma_{\zeta,i}^{(2)}$ to the data
in that single fiducial volume. We also demonstrate that the measured central
value, or limits, for $\overline{\Lambda}^{(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,

$\displaystyle\sigma(b_{T})\equiv\int\\!\mathrm{d}\Phi_{B}\,A(\Phi_{B})\,\sigma(Q,Y,b_{T})\,,$
(4.8)

which is possible because we consider the $q_{T}$-independent Born acceptance
$A(\Phi_{B})$. We then simply have
$\mathrm{d}\sigma/\mathrm{d}q_{T}=S[\sigma(b_{T})](q_{T})$. Our goal in the
following is to break down the quadratic contribution $\sigma^{(2)}(b_{T})$ to
$\sigma(b_{T})$ into two terms
$\sigma^{(2,\gamma_{\zeta})}+\sigma^{(2,\Lambda)}$ with a single, overall
coefficient each encoding the nonperturbative physics. This is straightforward
for the contribution from the rapidity anomalous dimension,

$\displaystyle\sigma^{(2,\gamma_{\zeta})}(b_{T})=\gamma_{\zeta,q}^{(2)}b_{T}^{2}\int\\!\mathrm{d}\Phi_{B}\,L_{Q^{2}}\,A(\Phi_{B})\,\sigma^{(0)}(Q,Y,b_{T})\,,$
(4.9)

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

$\displaystyle\sigma^{(2,\Lambda)}(b_{T})=2\overline{\Lambda}^{(2)}b_{T}^{2}\,\sigma^{(0)}(b_{T})\,,$
(4.10)

where to match eq. (4.2) the average coefficient must be defined by

$\displaystyle\overline{\Lambda}^{(2)}\\!=\\!\frac{\int\\!\\!\mathrm{d}\Phi_{\\!B}A(\Phi_{\\!B})H(Q^{2},\mu)\sum_{i,j}\limits\sigma^{B}_{ij}(Q)f_{i}^{(0)}\\!(x_{a},b_{T},\mu,\zeta_{a})f_{j}^{(0)}\\!(x_{b},b_{T},\mu,\zeta_{b})\bigl{[}\Lambda_{i}^{(2)}(x_{a})\\!+\\!\Lambda_{j}^{(2)}(x_{b})\bigr{]}}{\int\\!\\!\mathrm{d}\Phi_{\\!B}A(\Phi_{\\!B})H(Q^{2},\mu)\sum_{i,j}\limits\sigma^{B}_{ij}(Q)f_{i}^{(0)}\\!(x_{a},b_{T},\mu,\zeta_{a})f_{j}^{(0)}\\!(x_{b},b_{T},\mu,\zeta_{b})}.$
(4.11)

The fact that it is simply a ratio of (weighted) leading-power integrals is
due to the linear impact of $\Lambda_{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
$\overline{\Lambda}^{(2)}$.

However, given the expected sensitivity to $\overline{\Lambda}^{(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)
$\overline{\Lambda}^{(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 $\mu_{0}\sim 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

$\displaystyle\overline{\Lambda}^{(2)}\stackrel{{\scriptstyle\text{NLL}}}{{=}}\frac{\int\\!\mathrm{d}\Phi_{B}\,A(\Phi_{B})\,\sum_{i,j}\sigma^{B}_{ij}(Q)\,f_{i}^{(0)}(x_{a},\mu_{0})\,f_{j}^{(0)}(x_{b},\mu_{0})\bigl{[}\Lambda_{i}^{(2)}(x_{a})+\Lambda_{j}^{(2)}(x_{b})\bigr{]}}{2\int\\!\mathrm{d}\Phi_{B}\,A(\Phi_{B})\,\sum_{i,j}\sigma^{B}_{ij}(Q)\,f_{i}^{(0)}(x_{a},\mu_{0})\,f_{j}^{(0)}(x_{b},\mu_{0})}\,.$
(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 $\Lambda_{i,j}^{(2)}(x_{a,b})$ in the numerator. The
residual dependence on $\mu_{0}$ can be used as a handle on the theory
uncertainty during this recasting step, where $\mu_{0}\sim b_{0}/b_{T}$ should
be picked near the $b_{T}$-space Sudakov peak, $b_{T}\sim
0.2\,\mathrm{GeV}^{-1}$.

Figure -289: 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 N3LL cross section. Our averaging procedure for the intrinsic TMD PDF
coefficient $\overline{\Lambda}^{(2)}$ is defined in eq. (4.11) and discussed
there. We generate mock data (solid lines in bottom row) as described around
eq. (4.15). For clarity the uncertainty band for $\Delta_{\mathrm{cut}}$
(light blue band in bottom and middle row, gray in bottom row) is only shown
for the central fit template in the bottom row. For nonzero values of the
nonperturbative parameters of interest, the band looks essentially identical.
All OPE results use a third-order functional $S^{(3)}$.

In total, our $b_{T}$-space input for a given fiducial volume integrated over
hard phase space reads

$\displaystyle\sigma(b_{T})=\sigma^{(0)}(b_{T})\Bigl{\\{}1+b_{T}^{2}\Bigl{(}2\overline{\Lambda}^{(2)}+\gamma^{(2)}_{\zeta,q}L_{Q^{2}}\Bigr{)}\Bigr{\\}}+\mathcal{O}\bigl{[}(\Lambda_{\mathrm{QCD}}b_{T})^{4}\bigr{]}\,.$
(4.13)

Applying our truncated functionals at order $n=3$ and exploiting their
linearity, we have,

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}q_{T}}$
$\displaystyle=S^{(3)}\big{[}\sigma^{(0)}(b_{T})\big{]}(q_{T})+2\overline{\Lambda}^{(2)}\,S^{(3)}\big{[}b_{T}^{2}\sigma^{(0)}(b_{T})\big{]}(q_{T})+\gamma_{\zeta,q}^{(2)}\>S^{(3)}\big{[}b_{T}^{2}L_{Q^{2}}\sigma^{(0)}(b_{T})\big{]}(q_{T})$
$\displaystyle\quad+\frac{1}{q_{T}}\mathcal{O}\Bigl{[}(q_{T}b_{T}^{\mathrm{cut}})^{-\frac{5}{2}},\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{q_{T}}\Bigr{)}^{4}\,\Bigr{]}\,.$
(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

$\displaystyle\mathrm{d}\Phi_{B}=\mathrm{d}Q\,\mathrm{d}Y\mathrm{d}\cos\theta\,\frac{3}{8}\bigl{(}1+\cos^{2}\theta\bigr{)}\,\delta(Q-m_{Z})\,\delta(Y)\,,\qquad
A_{B}=1\,.$ (4.15)

The recasting relation in eq. (4.12) in this scenario reduces to a weighted
flavor average at $x_{a,b}=m_{Z}/E_{\mathrm{cm}}$. We plot the impact of
$\gamma_{\zeta}^{(2)}$ ($\overline{\Lambda}^{(2)}$) in the left (right) column
of figure -289, where in the first two rows we again interpret the impact of
the parameters as a model-independent nonperturbative uncertainty estimate and
compare it to two recent model-based 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 $\Delta_{\mathrm{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].

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 -289 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
$\gamma_{\zeta,q}^{(2)}$ on the bottom-left and $\overline{\Lambda}^{(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.2), shown by solid lines in the bottom row of figure -289. To do so
we use a $b^{*}_{\mathrm{Pavia}}$ prescription on $\sigma^{(0)}$, which as
demonstrated in section 2.5 modifies the cross section beyond our working
order in eq. (4.2). 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,

$\displaystyle f_{\mathrm{NP}}(Q,b_{T})=\begin{cases}\exp\big{(}-\lambda
L_{Q^{2}}^{k}b_{T}^{2}\big{)}\,,\quad&\lambda>0\,,\\\
\big{(}1+2\lvert\lambda\rvert
L_{Q^{2}}^{k}b_{T}^{2}\big{)}\exp\big{(}-\lvert\lambda\rvert
L_{Q^{2}}^{k}b_{T}^{2}\big{)}\,,\quad&\lambda<0\,,\end{cases}$ (4.16)

where $\lambda=-\overline{\Lambda}^{(2)}$ and $k=0$ for the intrinsic TMD PDF,
and $\lambda=-\gamma^{(2)}_{\zeta,q}$ and $k=1$ for the rapidity anomalous
dimension.

Comparing the dashed and solid lines in figure -289, we find that the proposed
fit template approximates the model results well down to $q_{T}\sim
8\,\mathrm{GeV}$ before nonlinearities and oscillatory error terms kick in,
which suggests $q_{T}\geq 8\,\mathrm{GeV}$ as a possible lower bound for the
fit window. This is similar in spirit to searches for nonzero $d\geq 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 $\gamma_{\zeta,q}^{(2)}$ enters linearly in
both.

###  4.3 Overview of theory systematics

From figure -289, we can read off that to be sensitive to
$\lvert\gamma_{\zeta,q}^{(2)}\rvert\lesssim(0.35\,\mathrm{GeV})^{2}$ and
$\lvert\overline{\Lambda}^{(2)}\rvert\lesssim(0.7\,\mathrm{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 $\sigma^{(0)}(Q,Y,b_{T})$.

Based on the perturbative uncertainty estimates in ref. [67], the N3LL
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 N3LL′ accuracy, which we expect to lower the perturbative QCD
uncertainty by a factor of $\approx 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 $\mathcal{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.2), where the effect of
$\gamma_{\zeta,q}^{(2)}$ and $\Lambda_{i,j}(x_{a},x_{b})$ can be dropped as
higher order cross terms. In addition, there are genuine QCD corrections of
$\mathcal{O}(q_{T}^{2}/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\lesssim 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 $\mathcal{O}(\alpha_{\mathrm{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
$\sim\alpha_{em}\Gamma_{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 $\mathcal{O}(b_{T}^{2})$.) 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 $\mathcal{O}(\alpha_{s}^{2}m_{Q}b_{T}^{2})$, 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_{T}^{2}$ 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}\lesssim 1/b_{T}\lesssim 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
$\mathcal{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.

##  5 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 short-distance 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_{T}^{\mathrm{cut}}<1/\Lambda_{\mathrm{QCD}}$ in the Fourier transform, such
that the momentum-space TMDs at perturbative momenta
$k_{T}\gg\Lambda_{\mathrm{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_{T}^{\mathrm{cut}}$. 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_{T}^{\mathrm{cut}}$. 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_{T}^{\mathrm{cut}}$. 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
$\mu=\sqrt{\zeta}=k_{T}^{\mathrm{cut}}$, performing the renormalization group
evolution of the TMD PDF to these scales in $b_{T}$ space at N3LL 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_{T}^{\mathrm{cut}}\gtrsim
10\,\mathrm{GeV}$. In particular, the deviation as a function of the momentum
fraction $x$ is within $\sim 2\%$ for $10^{-3}\leq x\leq 10^{-2}$. These N3LL
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 uncertainty 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 $\mu,\zeta$, we find that the
deviation remains within a few percent in the vicinity of
$\mu=\sqrt{\zeta}=k_{T}^{\mathrm{cut}}$. Furthermore, the evolution of $\zeta$
at fixed $\mu=k_{T}^{\mathrm{cut}}$ 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 $\mathcal{O}(b_{T}^{2})$ nonperturbative contributions to the
Collins-Soper kernel, $\gamma_{\zeta,q}^{(2)}$, and the intrinsic TMD PDF,
${\overline{\Lambda}_{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
$\mathcal{O}[(1\,\mathrm{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(\lvert\vec{r}\rvert)$
(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<\lvert\vec{r}\,\rvert<r_{\mathrm{cut}}$, and provide a systematic expansion
of the Fourier transform at large $\vec{k}$ in inverse powers of
$r_{\mathrm{cut}}\lvert\vec{k}\rvert\gg 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_{T}^{\mathrm{cut}}$ 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_{T}^{\mathrm{cut}}$. 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\bigl{[}-(b_{T}/b_{\mathrm{max}}-a)^{-1}\bigr{]}$. 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\leq b_{T}<\infty$. 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 $\mu=k_{T}^{\mathrm{cut}}$ and large
enough $k_{T}^{\mathrm{cut}}$. 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