Towards Quasi-Transverse Momentum Dependent PDFs Computable on the Lattice

Transverse momentum dependent parton distributions (TMDPDFs) which appear in factorized cross sections involve infinite Wilson lines with edges on or close to the light-cone. Since these TMDPDFs are not directly calculable with a Euclidean path integral in lattice QCD, we study the construction of quasi-TMDPDFs with finite-length spacelike Wilson lines that are amenable to such calculations. We define an infrared consistency test to determine which quasi-TMDPDF definitions are related to the TMDPDF, by carrying out a one-loop study of infrared logarithms of transverse position $b_T\sim \Lambda_{\rm QCD}^{-1}$, which must agree between them. This agreement is a necessary condition for the two quantities to be related by perturbative matching. TMDPDFs necessarily involve combining a hadron matrix element, which nominally depends on a single light-cone direction, with soft matrix elements that necessarily depend on two light-cone directions. We show at one loop that the simplest definitions of the quasi hadron matrix element, the quasi soft matrix element, and the resulting quasi-TMDPDF all fail the infrared consistency test. Ratios of impact parameter quasi-TMDPDFs still provide nontrivial information about the TMDPDFs, and are more robust since the soft matrix elements cancel. We show at one loop that such quasi ratios can be matched to ratios of the corresponding TMDPDFs. We also introduce a modified"bent"quasi soft matrix element which yields a quasi-TMDPDF that passes the consistency test with the TMDPDF at one loop, and discuss potential issues at higher orders.

A characteristic feature of TMDs is that they depend on both the longitudinal momentum fraction x and transverse momentum q T carried by the struck parton. Currently, TMDs are only directly calculable for the perturbative regime q T Λ QCD , where they can be obtained in terms of collinear parton distribution functions (PDFs). In contrast, nonperturbative TMDs with small q T ∼ Λ QCD have only been extracted from measurement by performing global fits to a variety of experimental data sets, see e.g. Refs. [29][30][31][32][33][34]. Imprecise knowledge of these TMDs limits the predictive power at small transverse momenta and constitutes a significant theoretical uncertainty. At small q T the TMDs also provide a crucial window into the structure of the proton. It is therefore desirable to find a method to calculate them from first principles.
For nonperturbative objects, lattice QCD currently provides the only practical means for first-principle calculations, and studies have been performed in Refs. [35][36][37][38][39], where ratios of x-moments of quark TMDPDFs were determined. This gets around the fact that the x-dependence of the TMDPDFs is not directly accessible from the Euclidean lattice, since it involves light-cone correlations which depend on the Minkowski time, in direct parallel to the same issue for calculating the x-dependence of the longitudinal PDFs. These analyses use Lorentz invariance to connect the TMD nucleon matrix element of interest to spatial matrix elements that are accessible from the lattice.
In Refs. [40,41] the large momentum effective theory (LaMET) was proposed as a method to overcome the hurdle of calculating light-cone quantities by relating them to time-independent quasi observables in a large-momentum nucleon state. The latter can be directly calculated on the Euclidean lattice and matched to the corresponding light-cone quantity through a factorization theorem that is based on a systematic expansion in the nucleon momentum. For example, the unpolarized collinear PDF of a proton moving along the n-direction is usually defined in the MS scheme as where i is the flavor index and the square bracket with subscript µ denotes that the operator is renormalized at scale µ. For later use we introduce the lightlike reference vectors n µ = (1, 0, 0, 1) andn µ = (1, 0, 0, −1). The light-cone coordinates are defined as b ± = b 0 ∓ b 3 , and p(P ) denotes a proton state with momentum P µ = (E, 0, 0, P z ). The path-ordered light-cone collinear Wilson line is To calculate the PDF f i (x, µ) using LaMET, one starts from a quasi-PDF, which is defined from equal-time correlation functions, where the operator is renormalized with a lattice friendly renormalization scheme andμ denotes a corresponding scale. The spacelike Wilson line is and Γ = γ 0 or Γ = γ 3 as they both belong to the same universality class of operators that can be related to the PDF through an infinite Lorentz boost along the z direction [42]. Unlike the PDF f i (x, µ) that is boost invariant, the quasi-PDFf i (x, P z ,μ) depends nontrivially on the nucleon momentum P z . When P z is large compared to Λ QCD as well as the nucleon mass M , the quasi-PDF satisfies the following factorization theorem [40,41,[43][44][45], (1. 5) where C ij are perturbative matching coefficients, and O(Λ 2 QCD /P 2 z , M 2 /P 2 z ) terms are power corrections. Here f j (y, µ) for −1 < y < 0 corresponds to the anti-quark PDF.
Due to the interest in TMDPDFs it is natural to consider the extension of the LaMET approach to transverse momentum observables. Due to the required focus on spatial matrix elements for TMDPDFs, studies based on LaMET are actually related to the lattice methods developed in Refs. [37][38][39]. While applying LaMET to TMDPDFs might seem straightforward, the richer structure of TMD factorization, which we review in Sec. 2, actually makes this quite non-trivial. In contrast to the case for collinear factorization, TMD physics is plagued by so-called rapidity divergences and the need for combining collinear 1 proton matrix elements with soft vacuum matrix elements. Such soft matrix elements retain a minimal amount of information about both incoming protons (their direction and the color charge of the probing parton). The importance of the soft matrix elements to cancel the analog of rapidity divergences in the spatial matrix elements for TMDPDFs has been discussed in Ref. [92], which was aimed at constructing a new TMD factorization theorem for the Drell-Yan process in terms of the so-called quasi-TMDPDFs in LaMET. More recently, the matching relationship between the quasi-TMDPDF on a finite-volume lattice and standard TMDPDF was studied at one-loop order in Ref. [93], where it was shown that the finite lattice size regulates these divergences, thus eliminating the need to introduce a dedicated regulator in the lattice calculation. However, we argue here that the final result of Ref. [93], which agrees at one loop with one of the results studied here, cannot be used for nonperturbative b T or q T , for reasons that will be discussed in detail.
In this work we consider the problem of constructing a quasi-TMDPDF that can be used to study the TMDPDF for nonperturbative q T ∼ Λ QCD . Here the quasi-TMDPDF must be chosen such that it can be calculated with lattice QCD, agrees with the physical TMDPDF in the infrared, and differs only by short distance contributions from the ultraviolet. This is required in order for a matching equation connecting the quasi-TMDPDF and TMDPDF to exist, analogous to Eq. (1.5), with a coefficient C that is not sensitive to infrared physics. These requirements still leave some freedom in the construction of a suitable quasi-TMDPDF. We therefore propose to test which quasi-TMDPDF definitions are feasible by carrying out a perturbative study of the infrared logarithms of q T , or equivalently the transverse position b T , which must agree order by order in α s with the same logarithms in the TMDPDF. Collinear infrared divergences related to the collinear momentum fraction x must of course also agree between TMDPDF and quasi-TMDPDF. We construct our quasi-TMDPDFs from distinct collinear proton and soft vacuum matrix elements, in direct analogy with the TMDPDF, where the quasi adjustment is made through the form of the operators appearing in these matrix elements. We carry out our study of infrared logarithms separately for the collinear and soft matrix elements, thus enabling us to separately probe the form of the corresponding collinear and soft operators. We also discuss in detail the role of rapidity divergences, which play an important role in the construction of TMDPDFs on the light-cone, which for the quasi-TMDPDF are regulated by having Wilson lines of finite length L, as required on lattice. The cancellation of the leftover L-dependence constrains how one must combine quasi-soft and quasi-collinear matrix elements. This paper is structured as follows. In Sec. 2, we review the TMD factorization theorem, discussing in particular the role of rapidity divergences and regulators, and how the TMDPDF is constructed from combining a proton matrix element and a vacuum soft matrix element. We then discuss in Sec. 3 the expected form of a matching relation between TMDPDF and quasi-TMDPDF, how Wilson lines of finite length naturally regulate the analog of rapidity divergences in lattice calculations, and discuss the quasi-construction of the matrix elements required to define the quasi-TMDPDF. In Sec. 4, we explicitly test the suggested quasi-TMDPDF by comparing to the TMDPDF at one loop, showing that the simplest quasi-collinear proton matrix element fails the infrared test. The simplest attempt of constructing a quasi-soft matrix element also fails this test. Finally, combining these into the simplest quasi-TMDPDF also gives a result that fails this test. To resolve this issue we introduce a bent quasi soft function, which leads to a quasi-TMDPDF whose infrared divergences properly match those of the TMDPDF at one-loop. Discussion of the implications of these results are given in Sec. 5, including issues that may still spoil the bent quasi soft function construction at higher orders, and how the main issues can be avoided by only studying ratios of TMDPDFs in impact parameter space. We conclude in Sec. 6. Further details that are important for our analysis are provided in appendices. In appendix A we summarize our notations and conventions for light-cone coordinates and MS, and contrast them with another popular convention used in the literature. In appendix B we discuss different schemes for TMDPDFs that are used in the literature, demonstrating that they all satisfy the general functional form for the TMDPDF that we use for our analysis. In appendix C we provide details on the perturbative calculation of the quasi proton matrix element.

Review of TMD Factorization
A precise understanding of the TMD factorization theorem and its different formulations in the literature is important to properly connect a lattice determination of the TMD-PDF to the phenomenological TMDPDF. In this section, we provide a detailed review of TMD factorization and set up a general notation for the definition of the TMDPDF that encompasses most of the available definitions in the literature.
For simplicity we consider TMD factorization in the context of the production of a color-singlet final state F with invariant mass Q, rapidity Y and transverse momentum q T = | q T | in the scattering of two energetic protons moving along the n µ = (1, 0, 0, 1) and n µ = (1, 0, 0, −1) directions. Examples of such cross sections include Drell-Yan, W , and Higgs production. In the limit q T Q the cross section can be factorized as The factorized cross section receives power corrections in q 2 T /Q 2 and Λ 2 QCD /Q 2 , and first studies of their perturbative and nonperturbative structure have been performed in Refs. [117][118][119]. Importantly, Eq. (2.1) remains valid for nonperturbative q T ∼ Λ QCD . The factorization is most elegantly written in impact parameter space, with b T being Fourier conjugate to the measured transverse momentum q T . In Eq. (2.1) i, j are parton indices for quark flavors or gluons, 2 and H ij is the hard function containing virtual corrections to the underlying hard process ab → F . The x a,b are the fractions of the proton momenta carried by the partons i and j participating in the hard collision, so x a = Qe Y /E cm and x b = Qe −Y /E cm where E cm is the center of mass energy of the pp collision. Finally, ζ a and ζ b are related to the momenta of the struck partons and always satisfy ζ a ζ b = Q 4 , and hence can be written as where P µ a = P − a n µ 2 are P µ b = P + bn µ 2 are the momenta of the protons, such that P − a P + b = E 2 cm , and y n is a parameter that encodes scheme dependence. The small transverse momentum q T Q of the final state is generated from soft and collinear radiation off the incoming protons. In Eq. (2.1a), the collinear and soft radiation are described separately. B i and B j are collinear proton matrix elements that measure the transverse momentum originating from energetic radiation close to the nandn-collinear protons, and to distinguish them from the final TMDPDF we will follow the language of Ref. [120] and refer to B i and B j as beam functions. The soft function S i ≡ S ij is a vacuum matrix element that encodes soft exchange between the incoming partons. It only differs for gluons, S g ≡ S gg , and quarks, S q ≡ S qq , but is independent of the light quark flavor. S i tracks the direction of both incoming partons, and hence depends on both light-cone directions n andn.
Note that H ij only depends on the renormalization scale µ, while B i , B j and S depend on µ and an additional rapidity renormalization scale ν. 3 This arises because the matrix elements not only suffer from UV divergences, but also from so-called rapidity divergences that require a dedicated regulator [94,109,112,115,[121][122][123]. We will discuss their physical origin in more detail in Sec. 2.2. Denoting the rapidity regulator generically as τ and the associated scale by ν, one can define UV and rapidity-renormalized beam and soft functions as For the soft function, this is straightforward. For the beam function, we note that the renormalized beam function depends on ζ ∝ (xP − ) 2 , where the proportionality is scheme-dependent, i.e. depends on the precise definition of τ . Secondly, we define B i as a purely collinear matrix element, while in practical calculations it often contains a soft contribution.
To avoid double counting with the soft function, one has to subtract out this overlap from the unsubtracted beam function B (unsub) i , and the subtraction factor is denoted as S 0 i in Eq. (2.3). (In SCET, this is referred to as soft zero-bin subtraction [124].) Since S 0 i describes soft physics, it cannot depend on the large momentum xP − . As indicated by the notation for S 0 i (b T , , τ ), this factor is often equivalent or closely related to the soft function S i itself, but this is not always the case. In particular, one can find rapidity renormalization schemes which (a) have no zero-bin subtraction, S 0 i = 1 [115], (b) where the zero-bin is equal to the soft function, S 0 i = S i [112,125], and (c) where the subtraction is given by a combination of soft functions in different non-lightlike directions [100].
In practice, one often combines the beam and soft functions to yield two separate TMDPDFs f TMD i and f TMD j , as used in the factorized cross section in Eq. (2.1b). This can be achieved either by combining the renormalized beam and soft functions, 5) or by combining the bare functions and performing the UV renormalization afterwards, In Eq. (2.5), the rapidity renormalization scale ν cancels between B i and S i , leaving only a dependence on ζ in the TMDPDF. Likewise, in Eq. (2.6) the dependence on the regulator τ cancels between B i and S i . Note that Eq. (2.6) has been written in terms of the unsubtracted beam function B (unsub) i , as is common in the literature, and for ease of notation we have combined the zero-bin subtraction S 0 i and the soft function S i to define the soft factor (2.7) As before, the ζ scale arises as a rapidity-regulator dependent function with the general form ζ ∝ (xP − ) 2 . The need to regulate (and renormalize) the rapidity divergences in the beam and soft functions has led to several different definitions in the literature. The original derivation by Collins and Soper used a non-lightlike axial gauge [94], whereas in more recent work Wilson lines off the light cone are employed [100,126]. The regulators used in the SCETbased approaches are the analytic regulator acting on eikonal propagators [127][128][129], the η-regulator inserted into Wilson lines on the light-cone [115,123], the δ-regulator which adds mass-like terms to eikonal propagators [112,130], and the exponential regulator inserted into the phase space [125]. The definition of the TMDPDF in terms of bare beam and soft functions, Eq. (2.6), is used in Refs. [100,112,113,126], where renormalized beam and soft functions were not defined. Both the bare and renormalized forms can be used in Refs. [115,125]. In the analytic regulator approach of Refs. [109,111], the soft function vanishes and hence one can only define a product of rapidity-finite TMDPDFs, not individual TMDPDFs. We provide a more detailed discussion of these regulators in appendix B, including explicit one-loop results for the quark beam function illustrating the combination of beam and soft functions into the TMDPDF.
Note that in the formulation of Refs. [94,126], the TMDPDF depends on an additional parameter ρ that cancels against the hard function. This definition is not included in our general setup here, since we have fixed the hard function to be in the MS scheme. This choice then fixes the scheme for the product of TMDPDFs in order to obtain a schemeindependent cross section, and hence does not allow room for an extra parameter like ρ. However, these definitions can be perturbatively related to those covered by our discussion, see Refs. [116,131].
In this work, we focus for simplicity only on the quark TMDPDF. For a hadron h moving along the n direction with momentum P , the bare beam and soft function are defined as Here the bracket [· · · ] τ denotes that the operator inside is considered by implementing the rapidity regulator τ . For clarity, we denote the Wilson lines in the n-collinear and soft matrix elements by W and S X , respectively, and both are defined by path-ordered exponentials. One needs both lines of infinite length along the light-cone, as well as finite-length gauge links with transverse paths, The Wilson paths for matrix elements in Eqs. (2.8) and (2.9) are shown in the (b + , b − , b T ) plane in Fig. 1. Note that the transverse gauge links at light-cone infinity create a closed path for the soft function, and a connected path between quark fields for the beam function, thus yielding matrix elements that are gauge invariant for both B q and S q . In nonsingular Figure 1: Graphs of the Wilson line structure of the n-collinear beam function B q (a) and the soft function S q (b), defined in Eqs. (2.8) and (2.9). The Wilson lines (solid) extend to infinity in the directions indicated. Adapted from Ref. [125].
gauges such as Feynman gauge where the gluon field strength vanishes at infinity, the transverse gauge links can often be neglected, but are known to be important in certain singular gauges, see e.g. Refs. [132][133][134][135]. These gauge links are also important for constructing the analogues of Eqs. (2.8) and (2.9) on a finite-size lattice.
Lastly, note that extracting TMDPDFs from lattice QCD (or experiment) is only necessary for nonperturbative For perturbative values, one can instead perform an OPE to match the TMDPDF, or equivalently the beam function, onto the collinear PDF, Here, C ij are perturbative matching kernels that are known to NNLO [136][137][138][139][140][141], and even to N 3 LO for the soft contribution [142]. The nonperturbative input is given solely by the standard longitudinal PDFs. Throughout this work, we will limit our discussion to the case b T ∼ Λ −1 QCD , where Eq. (2.12) cannot be applied.

Rapidity Divergences in TMDs
Quantum corrections to the beam and soft function defined in Eqs. (2.8) and (2.9) lead to two types of divergences: ordinary UV and IR divergences that can be regulated using dimensional regularization (and if desired a different IR regulator), and so-called rapidity divergences requiring a dedicated regulator [94,109,112,115,[121][122][123][124]. Rapidity divergences arise because nandn-collinear beam functions and the soft function, or equivalently the collinear and soft Wilson lines, are defined to describe modes with momenta scaling as where we use the light-cone notation p = (p + , p − , p T ), see appendix A, and λ ∼ q T /Q 1. This is illustrated in Fig. 2, where the orange dots denote the dominant region for the nandn-collinear modes, and the green dots denotes that for the soft modes. Hard modes mediating the underlying hard process are shown in blue. Since collinear and soft momenta have the same virtuality which in Fig. 2 is shown by lying on the same hyperbola, there is an "overlap" between soft and collinear momentum regions. This gives rise to divergences which cannot be regulated purely in dimensional regularization, and hence requires a dedicated regulator. This is indicated by the dashed lines in Fig. 2 that split the hyperbola. These divergences can also be understood as a conformal mapping from UV divergences in soft multi-parton scatterings [143,144], To give a concrete example of how rapidity divergences appear in perturbative calculations, consider the following integral that appears in calculations of the soft factor, Here the integrand only depends on the product k + k − . Singularities as k + k − → 0 or k + k − → ∞ are clearly regulated by dimensional regularization, while the integration over k − /k + is unconstrained, leading to a divergence that requires a separate regulator. Since k − /k + = e 2y is directly related to the rapidity y of k, these are often referred to as rapidity divergences that require an additional rapidity regulator. Alternatively, these are sometimes referred to as light-cone divergences, as they can also be avoided by displacing the light-cone propagators 1/k + and 1/k − away from the light-cone. To render the beam and soft functions well defined, integrals such as Eq. (2.15) require an additional regulator. A large variety of regulators has been suggested in the literature, see Refs. [94,100,112,115,119,123,[125][126][127][128][129][130]. The key idea in all regulators is to break the symmetry in k + ↔ k − . For example, in the η regulator of Refs. [115,123] one regulates the integral Eq. (2.15) by inserting a factor |2k z /ν| −η , such that (2.16) Divergences as k − /k + → 0 or k − /k + → ∞ are thus made manifest as poles in η. Similar to dimensional regularization, one introduces a new dimensionful scale ν to keep the regulator scaleless, and a parameter η to be taken to zero at the end of the calculation, analogous to the limit → 0. In this way, renormalized beam and soft functions acquire an additional Figure 2: Illustration of the collinear (orange) and soft (green) modes contribution to the q T measurement. The hard modes (blue) describe offshell modes producing the energetic final state. The dashed lines indicate that the degeneracy in p + p − has to be resolved to properly define separate collinear and soft functions.
scale dependence, namely ν, or equivalently the bare functions entering Eq. (2.6) depend on η. The chosen rapidity regulator determines the precise form of this dependence on ν and η, respectively. As discussed, this dependence cancels in the combination f TMD i = B i √ S i , but the cancellation leaves a residual dependence on the scale ζ. Intuitively, the definitions of ζ a and ζ b in Eq. (2.2) reflects how the hyperbola in Fig. 2 is split between soft and collinear modes, see e.g. Ref. [113] for more details.

Evolution of TMDPDFs
The cross section in Eq. (2.1) must be independent of the unphysical scale µ and of the precise choice of ζ a and ζ b as long as ζ a ζ b = Q 4 . This induces renormalization group equations (RGEs) that encode the dependence of hard, beam and soft functions, or equivalently hard functions and TMDPDF, on these scales. Here, we only discuss the RGEs for the TMDPDF, which is the object of interest in this work. Details on the separate evolution of beam and soft functions can be found e.g. in Refs. [115,145] and in appendix B.3. The RGEs for f TMD where Γ i cusp [α s ] is the cusp anomalous dimension. The second equation in Eq. (2.17) is known as the Collins-Soper equation [94,95]. The subscripts µ and ζ on the anomalous dimensions γ i µ and γ i ζ denote the scale evolution they govern, and the superscript i = q, g distinguishes quarks from gluons. These anomalous dimensions are defined as Here, we used that the ν-dependence in f TMD = B √ S cancels and that B only depends on the combination ν 2 /ζ, and thus γ i ζ can be equivalently obtained from either f TMD or B or S. Since the definition of S, Eq. (2.9), is independent of the hadron state entering the TMDPDF and the light quark flavor, this immediately implies that γ i ζ is independent of the choice of hadron state and light quark flavor as well.
The all-order forms of the anomalous dimensions are given by where γ i [α s ] denotes the noncusp piece. Note that γ i ζ (µ, b T ) has an intrinsically nonperturbative component for b T ∼ Λ −1 QCD , independently of the scale µ, as is clear from Eq. (2.19). However, once it has been nonperturbatively determined at a scale µ 0 Λ QCD , it can be perturbatively evolved to any other scale µ Λ QCD using (2.20) The boundary term γ i ζ (µ 0 , b T ) is known to three loops for perturbative b T and µ 0 [125,142,143].
Combining the solutions to Eq. (2.17), the TMDPDF can be evolved from arbitrary initial scales (µ 0 , ζ 0 ) to the desired final scales (µ, ζ) through (2.21) Here, we have chosen to first evolve ζ 0 → ζ at fixed µ 0 , and then µ 0 → µ. Since the evolution must be path independent, other choices are also possible, see also Refs. [115,146] for a discussion of this path independence. Also note that due to the b T dependence of γ ζ , the TMD evolution is severely more complicated in momentum ( q T ) space [145], which is part of the reason why the factorization is commonly written in b T space. Eq. (2.21) is crucial to relate the TMDPDF at reference scales (µ 0 , ζ 0 ), where they are either measured or determined from lattice QCD, to the phenomenological scales (µ, ζ).
The boundary term f TMD thereby reducing all the nonperturbative input to these more well-studied PDFs. In this case, the RG-evolved TMDPDF Eq. (2.21) serves to resum large logarithms ln(µb T ) and ln(ζb 2 T ), which otherwise spoil the perturbative convergence of the matching kernel C ij in Eq. (2.12).
In this paper, we are instead interested in obtaining nonperturbative TMDPDFs valid also for the case b T ∼ Λ −1 QCD , where the boundary term f TMD i (x, b T , µ 0 , ζ 0 ) in Eq. (2.21) is obtained from lattice QCD (or equivalently is measured from experiment). In this case, the values µ 0 and ζ 0 are fixed by the lattice calculation (or measurement), and the role of Eq. (2.21) is to evolve f TMD i (x, b T , µ 0 , ζ 0 ) to the scales (µ, ζ) required in the phenomenological application. This is completely analogous to the determination of collinear PDFs, which are extracted at a reference scale µ 0 and then evolved via the DGLAP evolution. Similar to the DGLAP evolution, the µ evolution of the TMDPDF encoded in the first exponential in Eq. (2.21) is perturbative as long as µ 0 and µ are perturbative. In contrast, the ζ evolution intrinsically contains a nonperturbative component γ i ζ (µ, b T ) for b T ∼ Λ −1 QCD , even if the TMDPDF is extracted at perturbative (µ 0 , ζ 0 ). Thus one needs to determine both f TMD i (x, b T , µ 0 , ζ 0 ) and γ ζ (µ, b T ) nonperturbatively. In particular, the nonperturbative determination of γ ζ (µ, b T ) is a phenomenologically relevant task on its own, as it provides valuable information on its all-order structure. For example, it is known to suffer from renormalons at large b T [147]. A dedicated discussion of the extraction of γ ζ (µ, b T ) from lattice QCD, exploiting some of the results discussed in this paper, has been given in Ref. [148].

Towards Constructing quasi-TMDPDFs
The goal of this work is to define a quasi-TMDPDFf TMD which involves matrix elements that are calculable with lattice QCD, and which can be matched onto the TMDPDF f TMD that is relevant for collider phenomenology. As reviewed in Sec. 2, TMDPDFs are constructed from hadronic and vacuum matrix elements, namely the beam function B i and the soft function S i . Both of these functions are sensitive to infrared physics for b T ∼ Λ −1 QCD . According to the LaMET approach, to calculate the TMDPDF in lattice QCD we need to construct quasi observables with the same infrared physics, which we refer to as a "quasi beam function"B i and "quasi soft function"S i .
In general both beam and soft functions can be rapidity divergent. On the lattice, we will see that the analog of rapidity divergences are regulated by a finite length L of the Wilson lines, while in the lightlike case many regulators have been suggested (see Sec. 2 and appendix B for more details).
In principle, one could envision separately matching the quasi beam and quasi soft functions onto the lightlike beam and lightlike soft functions, and then combining the matched results into the physical TMDPDF. However, such individual matching results necessarily depend on the choice of rapidity regulators, and furthermore since these regulators break boost invariance, some of them will likely spoil the boost relation between quasi and lightlike functions. For the physical TMDPDF the choice of the rapidity regulator has a significant impact on the form of the beam and soft functions, including their infrared logarithms, so at minimum the individual matching would have to be worked out separately for different schemes. Finally, due to the fact that L plays the role of a rapidity regulator for quasi distributions on the lattice, there could be non-trivial L dependence in intermediate stages of the matching result. For these reasons taking an approach of individually matching beam and soft functions is not preferred.
As discussed in Sec. 2.1, it is equivalent to instead combine the unrenormalized beam and soft functions, in which case the rapidity divergences cancel, and then perform the UV renormalization. This should hold for the quasi functions as well. Hence the more straightforward approach, adopted here, is to combine the unrenormalized quasi beam and quasi soft functions into a quasi-TMDPDF, in analogy to Eq. (2.6). This approach also has the advantage of canceling out all L dependence in the quasi-TMDPDF calculation, up to power suppressed terms. Of course, if the resulting quasi-TMDPDF fails the infrared consistency test, by not yielding infrared logarithms that are consistent with those in the TMDPDF, then, as we will see, it will still be advantageous to examine the contributing quasi-beam and quasi-soft functions to determine where the issues lie.
An additional complication in the lattice computation is the appearance of Wilson line self energies. For a Wilson line along the direction v, the self energy is proportional to v 2 and thus vanishes in the light-cone case where n 2 =n 2 = 0. On the lattice, one works in Euclidean space where v 2 = 0, so one obtains nonvanishing Wilson line self energies that have to subtracted. Since the self energies are also proportional to the length of the Wilson lines, and L dependence cancels when combiningB i andS i , this gives a b zdependent contribution to be removed by the UV renormalization. Here b z is the Fourier transform variable to xP z . Since this contribution is multiplicative in b z position space, it is most natural to perform the UV renormalization in position space. We thus define the quasi-TMDPDF in the MS scheme analogous to Eq. (2.6) as Here,B i ≡B (unsub) i and∆ i S are the quasi beam and quasi soft function, which remain to be constructed. They are the analogs of the unsubtracted beam function and the soft factor in Eq. (2.6). We will always consider the unsubtracted beam function and absorb the zero-bin subtraction factor into∆ i S , so for simplicity we drop the superscript "(unsub)".Z uv is the lattice renormalization constant, andZ converts from the lattice renormalization scheme to the MS scheme. These schemes are typically distinct, and hence we distinguish the lattice renormalization scaleμ from the MS scale µ. The finite lattice spacing a takes the role of as an UV regulator. On the lattice, one also has to truncate the Wilson lines that enterB i and∆ i S at a finite length L. As we will discuss in Sec. 3.2, having a finite L regulates the analog of rapidity divergences on lattice and thus the τ dependence is replaced by additional dependence on L. The L dependence associated with both the rapidity regularization and finite length of Wilson lines cancels betweenB i and∆ i S , and hencef TMD i does not depend on L. (In practice, there will be a finite L dependence from power corrections that vanishes in the limit L → ∞, and we suppress these in our notation.) Finally,f TMD i also depends on the proton momentum P z , which encodes short distance ultraviolet effects like in the quasi-PDF, and in addition acts as the analog of its ζ scale.

Derivation of a schematic form of the matching relation
The main focus of this paper is to give constructions ofB i and∆ i S such thatf TMD i can be perturbatively matched onto the TMDPDF f TMD i . In this section we will assume such a perturbative matching exists in order to constrain the general form it has to take.
Given the physical scales present in the calculation, the matching relation is expected to take the schematic form Given the physical picture underlying LaMET, we can also write down the expected power corrections to the matching Eq. (3.2), obtained from the assumed hierarchy of scales This scaling is motivated as follows: First, the lattice box size should be large enough so that the Wilson line length L is the largest length scale and finite L effects and finite volume effects are suppressed. Second, we assume b T ∼ Λ −1 QCD to be a nonperturbative scale of the order of the characteristic size of transverse fluctuations of constituents in the hadron. Lastly, LaMET assumes large P z to approximate a lightlike correlator by boosting an equal-time correlator, so 1/P z should be the smallest length scale in the system other than the lattice spacing a. Since for power counting we take x as O(1) this also implies that b z ∼ 1/P z is small.
Note that for a matching formula like Eq. (3.2) to exist, C TMD ij must not depend on b T , since b T encodes infrared physics and is assumed to be a nonperturbative scale. This immediately implies the necessary condition that in perturbation theory, the b T dependence of both quasi-TMDPDF and TMDPDF must agree, up to power corrections, which will be our most stringent consistency test in the one-loop study in Sec. 4.
We can deduce further constraints on the matching relation from the Collins-Soper equation, and thus arrive at an improved schematic form for the matching relation. To do so we assume that ζ and xP z are independent variables, which could be achieved for example by considering a different hadronic momentumP z for the quasi-TMDPDF than the momentum P z used for the TMDPDF. In Eq. (3.2), the ζ dependence must then cancel between C TMD ij and f TMD j to yield a ζ-independent quasi-TMDPDFf TMD i , which implies that This clearly violates the requirement that C TMD ij must be independent of b T to carry out short distance matching. This mismatch occurs because the rapidity evolution governed by γ j ζ (µ, b T ) is also nonperturbative. To correct for this we therefore must consider the modified matching relatioñ where for brevity we suppress the power corrections which are the same as in Eq. (3.2). In Eq. (3.5), the ζ dependence cancels between the exponential and f TMD , at the cost that the kernel C TMD ij now depends on the auxiliary scaleζ =ζ(x,P z ). As indicated, this scale must be fixed in terms of x andP z , the relevant scales thatf TMD depends on, and hence does not technically add additional functional dependence to C TMD ij .
In order to interpret Eq. (3.5) as a true matching equation without any renormalization group evolution, it must be possible to make the exponential in Eq. (3.5) vanish, to yield a purely perturbative relation betweenf TMD and f TMD . Thus, a perturbative matching can only be possible if one can chooseζ(x,P z ) = ζ to cancel the Collins-Soper evolution to all orders in perturbation theory. From Eq. (2. 2) for f TMD j (y, b T , µ, ζ) we have ζ = (2yP z ) 2 (choosing y n = 0 without loss of generality). Accounting for dimensions we see that we must chooseP z ∝ P z , and since at tree level C TMD ij [x, y, . . .] = δ(1 − x/y) this fixes the constant of proportionality so thatζ = (2xP z ) 2 . With ζ ∝ y 2 andζ =ζ(x), the only possibility that allows for perturbative matching is to takeP z = P z and have C TMD ij [x, y, . . .] ∝ δ(1−x/y) to all orders. (We will confirm this explicitly at one loop in Sec. 4 below.) With this constraint the schematic relationship between quasi-TMDPDF and TMDPDF becomes multiplicative in x space, To derive a perturbative formula for C TMD ij then requires choosingP z = P z and ζ = ζ = (2xP z ) 2 . For this choice the effect of changing ζ in f TMD j is exactly balanced by a corresponding change of P z inf TMD i . Note that the lack of an integral over y in Eq. analogous to the fact that no such integral appears in the renormalization group equations for the TMDPDF, see Eq. (2.17).
We will demonstrate that a quasi-TMDPDF can be defined such that the use of Eq. (3.6) and matching with a short distance C TMD qq is fully satisfied at one loop for quark quasi-TMDPDF and TMDPDFs. However this is not automatic, and indeed we find that the most naive definition of a quasi-TMDPDF does not agree with Eq. (3.6). An all orders derivation of a formula like Eq. (3.6) would be required to completely address this relationship, and is left for future work.

Impact of finite-length Wilson lines
Before constructing quasi beam and soft functions, we discuss the impact of having infinitelylong Wilson lines in the definitions Eqs. (2.8) and (2.9) of the lightlike beam and soft functions. On the lattice, the finite lattice size prevents infinitely-long Wilson lines. Hence one has to truncate them at some finite length L, as illustrated in Fig. 3. Here, it is important to include transverse gauge links as required by the factorization theorem, which ensures gauge invariance of the soft and collinear matrix elements.
Naively, one might assume that any effect of L is suppressed as b T /L or b + /L for sufficiently large L. In practice, the finite L also regulates the analog of rapidity divergences on the lattice, and hence yields divergences as L → ∞. While this has the advantage of not having to implement a dedicated rapidity regulator in the lattice calculation, the drawback is that one cannot easily disentangle finite-L effects from rapidity divergences.
To explicitly illustrate that L < ∞ is a valid rapidity regulator, we consider the lightlike soft function S q (b T , , L), where we can easily compare to existing regulators in the literature. Intuitively, having a fixed finite length L for the collinear Wilson lines breaks boost invariance of the matrix elements, thereby distinguishing collinear and soft modes and regulating the rapidity divergences. To see this concretely, consider the Feynman rule for a Wilson line of size L stretching along the n direction, compared to its L → ∞ limit, In Sec. 2, an explicit example of rapidity-divergent integral was discussed, see Eq. (2.15). For finite L, the example integral changes to Here we see that possible divergences as either k ± → 0 are regulated by having finite L, and the leftover logarithmic divergence as either k ± → ∞ is taken care of by dimensional regularization.
In our construction of the quasi functions on lattice, we will replace the lightlike Wilson lines by spacelike Wilson lines, which affects the eikonal propagator, so the analog of Eq. (3.8) is Clearly, the exponentials regulate a possible divergence as k z → 0, and thus play a similar role as in the lightlike case. However, Eq. (3.9) contains a quadratic dependence on k z in the denominator, rather than the linear dependence on k + and k − in Eq. (3.8). Thus, we can also encounter linear divergences in L, as opposed to having only logarithmic divergences ln(L) in the lightlike case.

Example: Lightlike soft function at NLO
To give a concrete example of the effect of finite L, we consider in detail the lightlike soft function, defined in Eq. (2.9), at one loop. To account for the effect of finite lattice size, the Wilson lines along the n andn directions are truncated at Ln and Ln, respectively, and transverse gauge links are included, as shown in Fig. 3b. In Feynman gauge, there are four relevant diagrams, shown in Fig. 4, of which only (a) and (b) have rapidity divergences, while (c) and (d) do not.
Rapidity-divergent diagrams. Let us first discuss Fig. 4a, where a gluon connects the Wilson lines separated by the transverse displacement b T . Together with its mirror diagram, it is given by In the second line, one can see how keeping L < ∞ regulates divergences as the light-cone coordinates k ± approach zero, thereby regulating the whole integral. In the final result, the rapidity divergences are then reflected as a double logarithm in b T /L. The second rapidity-divergent diagram, Fig. 4b and its mirror diagram are independent of b T , and are given by Together, we obtain where we have defined b 0 = 2e −γ E . The rapidity logarithms ln(µL) are manifest, while the last term is an example of a finite-L contribution that vanishes in the limit b T L. Let us also compare this to the soft function using the δ regulator [113] [see also Eq. (B.27)] The two expression agree upon identifying 1/L 2 = e 2γ E δ + δ − , showing that for this function the two regulators are closely related.
Diagrams involving transverse gauge links. The Feynman rule for the transverse gauge links at offsets nL andnL are given by where, In Feynman gauge, this vertex can thus be neglected for L → ∞, and one would not consider Fig. 4c. For finite L, we instead obtain for Fig. 4c S (L) This result vanishes for L → ∞ (or b T L) as expected. The situation is more intricate for Fig. 4d. Again, for L → ∞ the diagram would not be considered. For finite L instead, the L dependence drops out and one obtains Here, the symmetry factor 1/2 is compensated by the mirror diagram of Fig. 4d. The relative minus sign of the momentum k in the vertices arises because k is incoming into one vertex and outgoing from the other. Due to this relative sign, the exponential factors cancel, yielding a nonvanishing result of the diagram. In particular, it is independent of L and thus does not vanish as L → ∞.
In order to obtain a result consistent with the known L → ∞ limit, where this diagram does not contribute, the transverse self-energy has to cancel with other diagrams to not give a physical contribution to the cross section. Indeed, when combining the unsubtracted beam function with the soft function and zero-bin subtraction into f TMD as in Eq. (2.6), we find that these transverse self-energies will exactly cancel.
Lastly, we remark that we have explicitly checked that the diagrams with the transverse gauge links are indeed necessary to ensure gauge invariance using a R ξ gauge, but the calculation is otherwise not instructive and hence not presented here.

Construction of the quasi beam function
Recall the definition Eq. (2.8) of the light-cone beam function, The Wilson lines W extend to light-cone infinity, where they are closed by W T in the transverse direction, see Fig. 1a.
Following the standard LaMET procedure, we define the quasi beam function analogous to the beam function by replacing the light-cone correlator with an equal-time correlator, which in particular includes replacing the n-collinear Wilson lines by Wilson lines along thê z direction. Due to the finite lattice size, they must also be truncated at a length L, where one needs to include transverse Wilson lines to ensure gauge invariance. The resulting Wilson line path is illustrated in Fig. 5a. The definition of the bare quasi beam function in position space thus reads Here, we also replaced γ − by the Dirac structure Γ, which can be chosen as either Γ = γ 0 or Γ = γ z , as booth can be boosted onto γ − . The Wilson lines of length L are defined by The transverse gauge links are given by Eq. (2.11). A crucial feature of Eq. (3.19) is that it is an equal-time correlation function, i.e. it neither depends on the time-dependent light-cone separation b + nor on the lightlike directions n andn as Eq. (3.18), which makes it computable on lattice. The physical picture underlying this specific Wilson line structure is that boosting a purely spatial separation along theẑ direction yields an almost lightlike separation. Concretely, for a Lorentz boost along theẑ direction with velocity v and γ = 1/ This is illustrated in Fig. 5b: The pure spatial separation (blue) is boosted along the orangedotted trajectory, thereby approaching the lightlike separation (red). Note that regardless of whether b z is positive or negative, it is always boosted onto the same lightlike axisn µ , as required for the n-collinear PDF. To boost onto n µ , one needs to reverse the boost, v > 0, as is appropriate for then-collinear PDF, since the corresponding proton is moving into the opposite direction. It is easy to check that by applying the Lorentz boost Eq. (3.21) to Eq. (3.19), one recovers the matrix element in Eq. (3.18). When evaluated in a large-momentum nucleon state, the quasi beam function defined in Eq. (3.19) is equivalent to the matrix element of an almost-lightlike correlator in a static nucleon state. According to LaMET [40,41], the quasi beam function is related to the beam function in Eq. (3.18) through a factorization formula which includes perturbative matching and nonperturbative power corrections determined by the large nucleon momentum. This has been proven rigorously for the collinear PDF [43,45,47]. For the TMDPDF the physical boost picture shown in Fig. 5b implies a relation for the bare operators. However we must test the extent to which this relation survives when rapidity regulators are implemented, which we will do in Sec. 4.2. For the TMDPDF such regulators are known to have a significant effect.

Construction of the quasi soft function
Recall the definition Eq. (2.9) of the bare TMD soft function, Note that this vacuum matrix element has no explicit time dependence, in contrast to the collinear matrix element Eq. (3.18). Time dependence only enters indirectly through the lightlike directions of the Wilson lines S n and Sn, which on its own prohibits a direct computation on lattice. To obtain a lattice-computable quasi soft function, it thus seems reasonable to follow the same logic as above and replace As before, the lattice computation also requires to truncate the Wilson lines at a length L, where they are joined by transverse gauge links. The most naive attempt of constructing a quasi version of the soft function Eq. (2.9) thus takes the form where the soft Wilson lines of finite length are given by The resulting Wilson line path is illustrated in Fig. 6a. Unfortunately, the physical boost argument that allowed us to relate spatial Wilson lines to lightlike Wilson lines in the quasi PDF [see Eq. (3.21)] does not apply to the quasi soft function. Since the soft function involves both light-cone directions n andn, it is necessary to simultaneously obtain them from boosting ±ẑ. However, this requires boosts of opposite signs, as illustrated in Fig. 6b. Note that if this were possible with a single boost, it would directly violate the boost argument forB, where it is essential that both positive and negative b z 's are boosted onto the same light-cone direction.
Despite the simple boost picture breaking down, one can still test whether the matching for the obtained quasi-TMDPDF in the form of Eq. (3.6) is possible, and we study this in the next section at NLO. Indeed, we find that for the naive quasi soft function the matching is spoiled by the structure of infrared b T dependence. In Sec. 4.5, we will suggest a modified quasi soft function that yields a quasi-TMDPDF which has the correct infrared b T dependence at one loop. Given the absence of an intuitive boost relation, a rigorous all orders proof for any such quasi-TMDPDF proposal will certainly be required before full confidence can be obtained.

One Loop Results
In this section we present explicit one-loop results for the TMDPDF, the quasi beam and naive quasi soft function, and their combination into the quasi-TMDPDF. Here, we work in the MS scheme, as opposed to considering renormalization schemes appropriate for lattice calculations as discussed in Sec. 3, since a pure MS calculation is fully sufficient to perturbatively test the matching relation. Furthermore, we limit ourselves to the quark PDF and neglect mixing with gluons for simplicity, which corresponds to considering a nonsinglet flavor combination. All results are calculated by evaluating the appropriate matrix elements for the (quasi) beam function with an on-shell external quark with momentum P µ = (P z , 0, 0, P z ).

Lightcone TMDPDF at one loop
The unrenormalized result for the TMDPDF at one loop is given by is the quark-to-quark splitting function and the subscript + denotes a plus distribution such that In appendix B, we show that Eq. (4.1) agrees with the vast majority of regulators used to define TMDPDFs.
As indicated, the divergence in the first line in Eq. (4.1) is of infrared origin and matches precisely the IR divergence in the collinear PDF, which is crucial for matching the TMDPDF onto the PDF for perturbative b T Λ −1 QCD , see Eq. (2.12). Likewise, it must be exactly reproduced by the quasi-TMDPDF for a matching relation to exist. The second line in Eq. (4.1) contains UV poles and constants, and the last line contains the b T dependence. Similar to the IR pole, the b T dependence must be identical in the quasi-TMDPDF in order for a perturbative matching for b T ∼ Λ −1 QCD to exist.

Quasi beam function
We first calculate the quasi beam function defined in Eq. (3.19) by evaluating the operator in a quark state with on-shell momentum P 2 = 0. Working in pure dimensional regularization and taking the physical limit L → ∞, P z → ∞, we can directly Fourier transform the result into x space. At one loop, there are four contributions, shown in Fig. 7. The calculation is quite lengthy and shown in detail in appendix C. The result is given bỹ As anticipated, it contains the same IR divergence as the TMDPDF, Eq. (4.1). Note the presence of a linear divergence in L/b T , which we interpret as the analog of a rapidity divergence. As discussed below Eq. (3.9), these divergences appear as power-law divergences. Thus, after regularization they yield a linear dependence on the regulator L, rather than a logarithmic dependence ln(L/b T ). This linear L/b T term will exactly cancel with a similar term in the quasi-soft factor when combiningB q and∆ q S as in Eq. (3.1). In order to directly match this quasi beam function onto the lightlike beam function, we require that the logarithmic b T dependence, arising from IR physics, must be equal between them. However, the b T dependence does not agree with any beam function known in the literature, see the results in appendix B which are summarized in Table 1 below. In particular, only in Collins' scheme with Wilson lines off the light-cone one has the correct double-logarithm − 1 2 ln 2 (b 2 T µ 2 /b 2 0 ), while in all schemes with Wilson lines on the light-cone this double logarithm is (at least partially) contained in the soft function. Even in Collins' scheme the single ln(b 2 T µ 2 /b 2 0 ) does not match up with the corresponding single logarithm in the quasi beam function. 4 Hence, for all the rapidity regulators used in the literature, which yield the same universal TMDPDF defined in Sec. 2, none are in agreement with the simple physical picture of relating beam function and quasi beam function. The Lorentz boost relation is spoiled by the presence of a rapidity regulator, which by construction is not boost invariant. Since it is well known that the choice of rapidity regulator can modify the logarithms of b T , one may still hope to find a regulator for the beam function which yields the same IR structure as the naive quasi beam function and thus yields a perturbative matching that agrees with the boost relation. However, the more important test is whether the quasi-TMDPDF can be matched to the TMDPDF, in which case the regulator dependence cancels. This requires considering the quasi soft function.

Naive quasi soft function
Next, we calculate the naive quasi soft function defined in Eq. (3.24). Working in Feynman gauge, there are six diagrams that contribute at NLO, shown in Fig. 8, where double lines represent Wilson lines and the labels denote the endpoints of the Wilson lines in position space. For later convenience, we distinguish diagrams where the gluon is exchanged between the +ẑ and −ẑ Wilson lines (upper row) and diagrams where the gluon is emitted between Wilson lines of the +ẑ sector alone (lower row). The latter is identical to the result for gluons exchanged within the −ẑ sector.
In Feynman gauge, the generic expression for a one-loop diagram in the MS scheme, parameterized by the spatial paths γ 1 and γ 2 of the Wilson lines, is Note that diagrams with the gluon attaching to the same line have an additional symmetry factor of 1/2. For example, for Fig. 8a one reads off the paths Together with its mirror diagram, this gives Note that this result also contains a divergence at = 1/2, which signals a power-law divergence ∝ µL in dimensional regularization, which is expected for the Wilson-line self energy. In pure dimensional regularization, these divergences are not visible when expanding at = 0, but on the lattice they explicitly arise and have to be canceled. For this reason, the beam function renormalization in Eq.
Summing Eqs. (4.5)-(4.7), we obtain the soft contribution from interactions between the two collinear sectors, where we have also taken the limit L b T for illustration. Interestingly, in this limit all dependence on L drops out, leaving only a pure logarithm in b T , the physical observable.
The three remaining diagrams in Figs. 8d-8f yield and their sum gives the contribution to the soft function from interactions within theẑ sector alone. The same result is obtained for the −ẑ sector, so we havẽ Note that this contains the same linear divergence in L/b T as the quasi beam function, Eq. (4.2). Interestingly, the logarithmic dependence on L cancels within each collinear sector of the quasi soft function. The full bare quasi soft function is then obtained by summing the contributions from Eqs. (4.8) and (4.12), We observe that the logarithmic b T dependence of the naive quasi soft function does not match that of any soft function in the literature, regardless of the employed rapidity regulator, see the comparison made below in Table 1. Here, this is not as surprising as for the comparisons for the beam function in Sec. 4.2, as there is no relation between soft and quasi soft function through a Lorentz boost even at the bare level, see Sec. 3.4.

Failure of the naive quasi-TMDPDF
We can now attempt to combine the quasi beam and naive quasi soft function to construct a quasi-TMDPDF in the form of Eq. (3.1). To fix the form of the soft factor∆ q S , we require that all divergences in L/b T cancel when combiningB q and∆ q S , analogous to the cancellation of rapidity divergences in the TMDPDF. Comparing the one-loop results in Eqs. (4.2) and (4.13), we deduce this to be∆ q S = 1/ Sq , sõ . (4.14) Note that this is similar to the δ regulator for the TMDPDF, where one has ∆ q S = 1/ √ S q .
Our result∆ q S = 1/ Sq for the quasi-TMDPDF is consistent with this, considering that we showed in Sec. 3.2.1 that the lightlike soft function using the finite-L regulator yields the same soft function as the δ regulator.
Here, we work purely in the MS scheme and in the physical limit, where the product Z Z uv is b z -independent, so for our particular one-loop study we can equivalently work with the longitudinal momentum space formulã Although our method of calculation is quite different, we note that our result in Eq. (4.16) fully agrees with the one loop calculation in Ref. [93] (up to trivial differences in our respective conventions for the MS scheme).
There is an important subtlety concerning the logarithmic dependence ln(2P z ) 2 in Eqs. (4.1) and (4.16), which arises from calculating matrix elements with an on-shell external quark of momentum P µ = (P z , 0, 0, P z ). In the ratio of the actual TMDPDF evaluated in a proton state, we have to replace this by P z → xP z , where xP z is the momentum of the struck parton (see also the discussion in Ref. [45]), so we obtaiñ (4.17) This is our final result for the quasi-TMDPDF which uses the natural quasi-beam function and the naive quasi-soft function.
To test whether a perturbative matching betweenf TMD q and f TMD q is possible, we need to UV renormalize both Eqs. (4.1) and (4.17) and study their difference: As expected, the explicit infrared poles in IR have canceled, as they arise entirely from the quasi beam function, which can be related to the beam function through a boost. However, the b T dependence off TMD q and f TMD q does not agree, leaving two uncanceled logarithms in Eq. (4.18). The second one multiplies a logarithm of ζ, and in fact is exactly the one-loop expansion of the Collins-Soper kernel in Eq. (2.21), as expected. This leaves The key problem with Eq. (4.21) is that it still contains a single infrared logarithm ln(b 2 T µ 2 /b 2 0 ) which is not associated with the Collins-Soper evolution, and thus cannot be eliminated in a similar fashion. Curiously, choosing ζ = (2xP z ) 2 /e would simultaneously cancel both logarithms in b T in Eq. (4.18), but the term involving a ln ζ is clearly related to the Collins-Soper kernel, whereas there is no clear relationship of the term ln(b 2 T µ 2 /b 2 0 ) with this evolution. Therefore we deem this choice with an extra e to be something that works at one loop by construction, but not a valid choice since it is very unlikely to continue to work at higher loop orders (unless there happens to be some unknown deep relationship). The presence of this extra ln(b 2 T µ 2 /b 2 0 ) thus indicates a failure of the naive quasi-TMDPDF to reproduce the same infrared physics at one-loop as required for the physical TMDPDF.
Our results can also be compared to those in Ref. [93], where the soft factor was not calculated separately, but immediately combined with the beam function to yield the quasi-TMDPDF. The final result obtained was a relation between the quasi-TMDPDF and TMDPDF at µ 2 = ζ = (2xP z ) 2 which was given as If we take our result in Eq. (4.18) and set µ 2 = ζ = (2xP z ) 2 , then we obtaiñ for the different regulators in the literature, as shown in Table 1. The dependence of the quasi constructionsB q ,∆ q S andf TMD q on L b is also shown in the lower part of the table. Here we only show the dependence on standalone factors of L b (providing references for the full expressions in the table caption). As discussed previously, the quasi functions do not match their lightlike counterparts in any of the regulators. In particular, the double-logarithm L 2 b only agrees with Collins' regulator. All the other regulators involve Wilson lines with light-like directions, and here the double logarithm is part of the soft function (it is split between these two in the exponential regulator). 5

Quasi TMDPDF using a bent soft function
The construction of the quasi beam function in Sec. 3.3 was motivated by the physical picture of boosting a spatial to a lightlike correlation function, while the (naive) quasi soft factor construction in Sec. 3.4 was simply the most straightforward attempt. This lead to a quasi-TMDPDF whose IR logarithms do not match those of the TMDPDF. However, there is significant freedom in constructing quasi functions on lattice, so we can consider alternate definitions with the goal of finding one which has the same infrared physics as the TMDPDF. When the quasi beam function and quasi soft function are combined, any dependence related to the method of regulating rapidity divergences (such as finite L) 5 Note that in none of these cases does one include the transverse self energy, which would add a single logarithm L b to Bq and −L b to ∆ q S = 1/ √ S q , see Eq. (3.16). Even after taking this into account, the mismatch persists. cancels. Since it was only the presence of rapidity regulators that causes problems for the physical boost argument for the beam function, one may infer that this issue is alleviated when considering the matching for the quasi-TMDPDF and TMDPDF. For this reason we will not try to adjust the definition of the quasi beam function here. However, the quasi soft function was not constructed based on a boost argument, and hence seems like the most likely culprit for the failure to match infrared logarithms. For the soft factors contained in the TMDPDF there are always two different spatial directions involved in the Wilson lines, while for our naive quasi soft factor there was only the z-direction. This motivates us to consider in this section a different "bent" quasi soft function which involves two spatial directions.
This need for this type of bent quasi soft function can also be motivated by studying the failure of the naive quasi-soft function to reproduce the IR physics needed for the TMDPDF in more detail. In particular, we can split the calculation of the naive quasi soft function into three distinct pieces, arising from gluon exchanges either within theẑ or −ẑ sector, or between them, as done in Sec. 4.3, (4.24) Physically, the first term is correctly boosted towards a n/n contribution by boosting with P z > 0, and likewise the second term is boosted towardsn/n for P z < 0. In practice, they are identical due to invariance under z ↔ −z. At one loop, the quasi-TMDPDF in Eq. (4.15) hence can be written as z/z exactly cancels the soft limit of the tadpole correction to the beam function. This can easily be verified by comparing Eqs. (4.12) and (C.51), from which one also sees that this subtraction cancels the divergence in L/b T . It remains to considerS (1) z/−z , given in Eq. (4.8). Its subtraction fromB q adds a single ln(b 2 T µ 2 /b 2 0 ), which is exactly the leftover logarithm found in the relation between the naive quasi-TMDPDF and TMDPDF in Eq. (4.18). In conclusion, it thus appears to be precisely the interaction between the +ẑ and −ẑ part of the naive soft function that is spoiling the matching of infrared logarithms.
With these motivations and observations we can define a bent quasi-soft function which gives a valid matching result between quasi-TMDPDF and TMDPDF at one-loop order. Crucially, it must still cancel the L/b T divergence in the quasi beam function and after combination with the quasi beam function produce the same logarithms in b T as the TMD-PDF. More concretely, we can demand the Wilson line structure in the z sector to match the soft-expanded (b z = 0) structure of the naive quasi beam function to ensure the cancellation of rapidity divergences. Given these restrictions we define the "bent" soft function wheren µ ⊥ is the transverse unit vector orthogonal to n µ ⊥ = b µ ⊥ /b T andẑ. Fig. 9 illustrates the Wilson line path in Eq. (4.26) and compares it to the path for naive quasi soft function defined in Eq. (3.24).
Above we deduced that the failure of a perturbative one-loop matching between quasi-TMDPDF and TMDPDF could be traced to soft diagrams mediating exchange between Wilson lines along the +ẑ and −ẑ directions. These diagrams precisely vanish for the bent soft function due ton ⊥ · n ⊥ =n ⊥ ·ẑ = 0, while all other diagrams are not affected by the new Wilson line paths. Hence the bent soft function at one loop yields As before, it is related to the soft subtration through∆ q S = 1/ S bent . This bent soft factor has precisely the infrared logarithms that are desired at one loop. 6 6 A soft factor with more than one transverse directions was also used in Ref. [92] where the goal was to express the physical factorization theorem directly in terms of quasi objects. They define a soft factor Using Eq. (4.15) to combine the bent quasi soft function from Eq. (4.27) together with the natural quasi beam function from Eq. (4.2) we obtain a new quasi-TMDPDF (4.29) Comparing this result to the TMDPDF at one loop yields where have again fixed ζ = (2xP z ) 2 as explained previously. Since there is no b T dependence on the RHS of Eq. (4.30), we see that all infrared logarithms of the TMDPDF are correctly reproduced by this quasi-TMDPDF construction at one loop. Thus this construction obeys the matching relation given in Eq. (3.6) with a one loop result for the matching coefficient that is given by (4.31) Here, we ignore possible mixing of quarks with gluons. Then since mixing of quark flavors can first arise at two loops, the one-loop coefficient is proportional to δ qq . This result provides a valid one-loop perturbative matching coefficient, which only depends on the hard scale of the struck parton, xP z . Assuming the validity of this quasi-TMDPDF construction beyond one loop, Eq. (4.31) can be used to match the lattice quasi-TMDPDF to the TMDPDF. To obtain the required input for this result one combines lattice calculations of the natural quasi beam function and bent quasi soft function to obtain a lattice quasi-TMDPDF, which is then converted into the MS scheme. Results for matching in more lattice friendly renormalization schemes should be straightforward to derive following a similar approach to the one used here (see e.g. [57,65]).

Results and Outlook
In this section, we briefly summarize the impact of our calculations in the previous sections for the matching between quasi-TMDPDF and TMDPDF, and what questions remain open for further study. Without relying on the existence of a quasi soft function that yields the correct infrared physics for a quasi-TMDPDF, we also discuss precisely what constraints on TMDPDFs can still be rigorously derived from lattice calculations.

Matching relation between quasi-TMDPDF and TMDPDF
The goal of this work was to establish a matching relation between the quasi-TMDPDF and TMDPDF analogous to the collinear PDF, where LaMET gives such a perturbative relation. However, the physical picture for the existence of such a matching relation is much more complicated than in the PDF case. For the beam function, the need for a non-trivial rapidity regulator on the TMDPDF side appears to spoil the simple boost correspondence between hadronic quasi and non-quasi matrix elements. We have confirmed that this is the case at one loop in Sec. 4.2 and Table 1 by making comparisons of the most natural quasi beam function with all modern TMDPDF definitions for the beam function. For the soft function, the vacuum matrix elements that appear necessarily involve two directions, and hence does not satisfy a simple boost relation to a quasi soft function even at the bare level. In Sec. 4.4 we computed the most naive quasi soft function at one loop, and showed that it yields a quasi-TMDPDF which does not have infrared logarithms that agree with those of the TMDPDF. Then in Sec. 4.5 we considered a modified bent quasi soft function at one loop, which when combined with the natural beam function gives our preferred quasi-TMDPDF definition. Its infrared logarithms at one loop properly agree with those of the TMDPDF, thus leading to a consistent matching formula at one loop.
The full utility of our preferred quasi-TMDPDF definition will depend on whether the correspondence between infrared logarithms continues to hold at higher orders in perturbation theory. For example, at two loops one can recouple the −ẑ and +ŷ sectors in the bent soft function, so these contributions will have to give contributions that match up correctly with corresponding contribution in the two loop soft function once its combined into the TMDPDF. Such calculations should be considered in the near future. More rigorously, one needs to show that our bent quasi soft function together with the natural quasi beam function yields a quasi-TMDPDF with the same infrared nonperturbative physics as the lightlike soft function, either nonperturbatively or at least to all orders in perturbation theory. This is obviously also an important avenue for future work.
With lack of further information, in order to proceed at the current time one can do one of two things, i) make the assumption that our preferred quasi-TMPDF holds nonperturbatively, or ii) assume that our bent quasi soft factor may not work to all orders, and see if interesting constraints can still be derived. We discuss these in turn.
In the case of i) we have the matching formula derived in Sec. 3.1, which for a nonsinglet (ns) quark flavor combination where there is no mixing reads: x, µ, P z ,ζ(x, P z ) = (2xP z ) 2 At this order, the kernel is diagonal in the quark flavors, as flavor mixing can first appear at two loops. Under the same assumptions one can proceed to carry out calculations for our preferred quasi-TMDPDF on the lattice, and with suitable renormalization and scheme changes, then use the matching relation in Eq. (5.2) to obtain the physical TMDPDF for the nonsinglet case.
In the case of ii) we assume that the matching between quasi-TMDPDF and TMDPDF is spoiled by a mismatch between the soft factors∆ q S and ∆ q S at higher orders. In this case, although we do not have a matching relation, we can still write down a formula relating the quasi-TMDPDF and TMDPDF. For the nonsinglet case without mixing it reads: Here, C TMD ns is by definition still a perturbative function. On the other hand, the function g S q (b T , µ) is nonperturbative and is associated with the failure of constructing a proper quasi soft function. It corrects the mismatch in the infrared physics. Since the quasi soft factor and soft factor which differ are flavor independent, there are at most two different g S i (b T , µ)'s, namely for quarks i = q and gluons i = g, and only g S q shows up for the non-singlet flavor case. Note that Eq. (5.3) is also satisfied if we use the naive quasi soft function construction. Summarizing our results obtained with the naive and bent quasi soft functions we have Interestingly, even with the less strong assumptions present in case ii) we can still extract non-trivial information about the TMDPDF by using ratios of distribution functions where the g S q factors cancel out. We discuss this further in the next section.

Ratios of TMDPDFs
While the presence of the soft sector seems to prohibit a straightforward computation of TMDPDFs on lattice, one can employ the fact that g S i (b T , µ) only differs for quarks and gluons but otherwise is flavor blind to try to construct ratios of quasi-TMDPDFs where g S i (b T , µ) cancels. To avoid possible mixing between quarks and gluons, here we only consider nonsinglet quasi-TMDPDFs such as ns = u−d, where Eq. (5.3) applies. The soft factor g S q cancels in any ratio of two (quasi) TMDPDFs with the same choice of b T and µ, so we havẽ For example, one can choose x 1 = x 2 and ζ 1 = ζ 2 to expose the Collins-Soper kernel γ q ζ as This allows one to determine the nonperturbative b T dependence of γ q ζ in terms of the quasi-TMDPDFsf TMD ns , and was proposed in Ref. [148]. Without this knowledge of γ q ζ , one is forced to choose x 1 P z 1 = x 2 P z 2 and ζ 1 = ζ 2 in Eq. (5.5) to cancel the Collins-Soper evolution kernel, .

(5.7)
Here, h 1 and h 2 for example refers to TMDPDF in different hadron states (which does not affect C TMD ns/h i ), or TMDPDFs of different spin structures (which can affect C TMD ns/h i ). The latter requires both spin structures to be either T -even or T -odd, as the soft function is not T -invariant, see e.g. Ref. [100]. Using the results given in Sec. 4, it is easy to confirm that the infrared logarithms cancel independently at one loop on the left and right hand sides of Eq. (5.7), confirming that it indeed is a matching equation at this order.
In Eqs. (5.5)-(5.7) we write ratios of quasi-TMDPDFs, which contains some choice for the quasi soft factor that satisfies Eq. (5.3). However, one can completely avoid the need for including a quasi soft factor in such ratios by employing Eq. (3.1), which yields . (5.8) In this ratio, the quasi soft factor cancels because it is independent of b z and P z 1,2 . Thus Eq. (5.8) removes the necessity to calculate a quasi soft matrix element on lattice. The leftover divergences as L → ∞ present in the individualB q functions also cancel between the numerator and denominator.
Note that our analysis of ratios of quasi-TMDPDFs here differs somewhat from that of Refs. [37][38][39]. In those references Lorentz invariance is used to directly relate ratios of b T -dependent spacelike TMDPDFs to the corresponding ratios of the physical TMDPDFs, analogous to our Eq. (5.7). However this is done by considering an adjustable spatial path for the beam functions, and taking the limit where this path approaches the light-cone. So far the focus has been on the case with b z → 0 corresponding to integrating over x. The relations they use do not require a non-trivial matching coefficient C TMD , but they do require a non-trivial adjustment of the path as the light-cone limit is taken. It would be interesting to consider a more detailed analysis of the difference between our approach and theirs.

Conclusion
In this paper we have studied the possibility to obtain quark TMDPDFs from lattice QCD using the LaMET approach. LaMET has been successfully applied to obtain collinear PDFs using a perturbative matching relation from quasi-PDFs, which are equal-time correlators evaluated in a highly-boosted hadron state. The construction of quasi-TMDPDFs is severely complicated by the presence of rapidity divergences that require a dedicated regulator and the need to combine the beam function B q , a collinear hadronic matrix element, with a soft factor ∆ q S defined through a soft vacuum matrix element. We have first discussed why the analog of rapidity divergences do not pose a problem for lattice calculations, as they are fully regulated by the finite length L of Wilson lines, both in the lightlike case of TMDPDFs and the equal-time case of quasi-TMDPDFs.
We then separately discussed constructions of the quasi beam functionB q and the quasi soft factor∆ q S . Since beam functions only depend on the flavor of the probed parton and the direction of the struck hadron, similar to normal PDFs, there is a natural definition of the quasi beam function as an equal-time matrix element following the LaMET procedure. With this definition the bare operator matrix element can be boosted onto the bare beam function operator matrix element. On the other hand, the soft function depends on the direction of both hadrons and thus can not be obtained from boosting a purely spatial matrix element.
We hence first studied a naive quasi soft function where one replaces the lightlike directions by n µ →ẑ µ andn µ → −ẑ µ . Both matrix elements separately suffer from rapidity divergences that are regulated by the length L, and we construct the naive quasi-TMDPDF f TMD q =B q∆ q S by demanding the cancellation of the L dependence. We also discussed implications of the Collins-Soper evolution on the matching relation between f TMD q and f TMD q , which (if it exists) requires one to fix the Collins-Soper scale ζ in the TMDPDF through the proton momentum P z in the quasi-TMDPDF, and for the relation not to involve a convolution in a momentum fraction.
We have carried out a detailed one-loop calculation to study whether a perturbative matching is feasible despite the failure of the physical boost picture. For this, we require that all logarithmic dependence on the transverse separation b T ∼ q −1 T ∼ Λ −1 QCD , assumed to be a nonperturbative scale, must agree between the lightlike functions and their quasi constructions. We find that this consistency test fails for the quasi beam function, the naive quasi soft function and the naive quasi-TMDPDF. For the quasi beam function, we interpret this to be due the need to regulate rapidity divergences, which necessarily breaks boost invariance and thereby invalidates the simple boost relation. Even after combining B q and∆ S q into the quasi-TMDPDFf TMD q , in which case all regularization dependence cancels, there is still a mismatch between TMDPDF and quasi-TMDPDF.
To fix this inconsistency we were motivated to consider a "bent" quasi soft function which involves an equal time operator with Wilson lines on two different spatial paths. At one loop, one can identify the diagrams that violate the boost relation in the soft function, which motivated the precise definition of a "bent" soft function. This leads to a quasi soft factor∆ q S which gives a quasi-TMDPDF which matches the b T dependence of the TMDPDF, establishing a perturbative matching relation at least up to one loop.
If our construction with a bent quasi soft function and the natural quasi beam function works beyond one loop order then it, for the first time, provides a method to fully access the complete physical TMDPDF from lattice QCD computations. Given the importance of such a construction, it is very important to further test this relationship beyond one loop, as we have emphasized repeatedly.
Even if this construction breaks down at higher order, one can still consider ratios of quasi-TMDPDFs, where any mismatch in the soft physics cancels. This enables one to study ratios of TMDPDFs. For example, this allows one to nonperturbatively determine the Collins-Soper kernel, which has been presented in detail in Ref. [148], or to study ratios of TMDPDFs using different hadron states or spin structures as described in Sec. 5. Such information is potentially useful to constrain TMDPDFs and thereby aid their extraction from experiment, particularly for TMDPDF spin structures or parameter ranges where only limited data is available.

A Notation and conventions
We briefly summarize the conventions used in this paper and compare them to different conventions used in the literature.
Lightcone coordinates. We use the SCET notation based on the two reference vectors n µ = (1, 0, 0, 1) ,n µ = (1, 0, 0, −1) , (A.1) which satisfy n 2 =n 2 = 0 and n ·n = 2. One can decompose any fourvector as where the plus and minus components are defined as and the transverse component is orthogonal to both reference vectors, n · k ⊥ =n · k ⊥ = 0. The Minkowski product of two fourvectors follows to be In particular, one has k 2 = k + k − − k 2 T . Another convention often used in the literature is k ± = (k 0 ±k z )/ √ 2, see e.g. Ref. [100]. In that notation, the Minkowski product is To translate results from this paper to that convention, replace k ± → √ 2k ∓ . Often, this leaves factor √ 2, that can be absorbed e.g. in integration measures.
Position space. We often write the collinear and soft matrix elements in position space, with the position denoted by b µ = (b 0 , b T , b z ). The Fourier transform to and from momentum space is defined as We use the same symbol for a function in f in both spaces, as we reserve the symbolf for the quasi-construction of f . In practice, we only perform the Fourier transform with respect to b + to obtain the momentum fraction x (b z for the quasi TMD), while we keep the transverse dependence in position space (often called impact parameter space), where the canonical logarithm is given by Renormalization scheme. Our results are expressed in the MS scheme. The associated renormalization scale µ ≡ µ MS is related to the MS scale µ 0 by Note that this differs from the convention used by Collins [100], where µ 2 = 4πµ 2 0 /Γ(1 − ), by terms of O( 2 ).

B Comparison of different schemes for TMD definitions
Here, we give more details on the various rapidity regulators employed in the literature that are compatible with the generic notation used in Sec. 2, where the TMDPDF was defined in Eq. (2.6) as and the Collins-Soper scale was given by Eq. (2.2), ζ = (xP − e −yn ) 2 = (xm P e y P −yn ) 2 .

(B.2)
We show how this formulation for the TMDPDF arises in the various regulators, including Wilson lines off the light-cone, the δ regulator, the η regulator, and the exponential regulator from Refs. [100,112,115,125] respectively. We also explicitly give one-loop results for the individual matrix element ingredients in the TMDPDF where available. In all cases, while the ingredients differ, the same universal result is obtained for f TMD q . We also give the correspondence to results for the analytic regulator of Ref. [109], where only the product of two TMDs are defined, but which also agree with this f TMD q .
At one-loop, the common result when evaluating the TMDPDF in an on-shell quark state of momentum P using pure dimensional regularization is Here, L b is defined in Eq. (A.6) and P qq (x) = (1 + x 2 )/(1 − x) is the regular part of the quark-quark splitting function. Absorbing the UV divergences in a multiplicative counter term one obtains the renormalized TMDPDF as The remaining 1/ IR pole here is of infrared origin and is the same collinear divergence that is present for the PDF. This correspondence between infrared divergences enables the TMDPDF to be matched on to the PDF for perturbative b T .

B.1 Wilson lines off the light-cone
In the modern definition by Collins [100], the lightlike Wilson lines are replaced by spacelike Wilson lines, This affects both the beam function as well as the soft factor, which is now a combination of soft matrix elements. From Eq. (13.106) in Ref. [100] we have for the n-collinear TMDPDF where the result in the last line was derived in Ref. [149]. Here, y P is the rapidity of the hadron (not the Wilson line direction), y B is the direction of the Wilson line as in Eq. (B.6) which acts as a rapidity regulator, and y n is a rapidity parameter that controls the split of soft radiation into the two TMDPDFs. The ζ scale is defined as in Eq. (B.2). To relate Eq. (B.7) to our notation in Eq. (B.6) one identifies τ = 1/(y B − y n ) and notes that so that the translation of functional forms is where B C q ≡f unsub q is the collinear matrix element as defined in Eq. (13.108) in Ref. [100] and S q C ≡S (0) is the soft function as defined in Eq. (13.39) therein. As usual, the dependence on m P and Λ QCD is implicit. Note that the two B q 's have the same number of arguments after using the fact that only the combination 1/τ − ln √ ζ appears in B q . Since ζ = (xm p e y P −yn ) 2 , one can obtain the Collins-Soper kernel as which is the analog of Eq. (2.18). This makes clear that the Collins-Soper kernel can be obtained as a differential equation involving f TMD , or B q , or ∆ q S independently. This fact is commonly exploited in the perturbative TMDPDF literature.

One-loop results
We are not aware of explicit individual one-loop results for B C q and S q C in the literature. Instead, in Refs. [100,150] the two functions are combined at the integrand level before integrating over the momentum of either real or virtual emissions. For example, in Ref. [100], where the soft factor is calculated in the context of the TMD fragmentation function, the limit y A,B → ±∞ is taken right away. The resulting divergent integral is then directly canceled against a divergent integral in the beam function, which for the TMDPDF then yields Eq. (B.5). The same order of operations was used for the calculation of the TMDPDF for double-parton scattering in Ref. [149].
We can calculate the soft function in the Collins scheme with Wilson lines along two spacelike vectors n A and n B . In Feynman gauge, purely virtual diagrams are scaleless and vanish, so we only need to take real emissions into account, as shown in Fig. 10. The two diagrams are given by Note that integral in in Eq. (B.12) is ill-defined due to a pinched pole singularity arising from the opposite signs of the i0 prescription in the two eikonal vertices, and we have evaluated it using the principal-value prescription. This is similar to the calculation in Ref. [151], where the soft function as defined in Ref. [94] was calculated in an axial gauge with v · A = 0 for a spacelike v. However, there the principal-value prescription is part of the axial gauge, while we employ it because a corresponding regularization in terms of Wilson lines for Feynman gauge is currently unknown, see also the discussion in appendix A of Ref. [151]. In the calculation of Ref. [126], the pinched pole singularities do not pose a problem due to employing timelike reference vectors n A and n B . Taking the mirror diagrams into account and taking the limit |y A − y B | 1, we obtain which after renormalization agrees with the result for timelike Wilson lines in Ref. [126]. Thus at one loop using Eq. (B.9) with τ = 1/(y B − y n ) we find 14) The one-loop result for f TMD q was calculated in Ref. [150] and agrees with Eq. (B.5) after adjusting for the difference in definitions of the MS scheme, see appendix A. Hence we can also deduce the bare beam function in this scheme, using Using Eq. (B.8), one can see that this result only depends on y P − y B .

B.2 δ regulator
The δ regulator was originally introduced in Ref. [112] and then modified in Refs. [140,152] to be applicable for NNLO calculations. Note that their notation for the light-cone momenta x ± are related to our notation by x ± ↔ x ∓ . For consistency of this paper, we convert results for the δ regulator in the literature to our notation.
Here, we follow the definition of the TMDPDF as presented in Ref. [113], which generalizes the definition in Ref. [112]. At one loop, the δ regulator amounts to modifying propagators in the full theory as are the momenta of the struck quark in the n-collinear andn-collinear proton, respectively. One can equivalently modify the eikonal propagators as (Additional integer m > 0 multiplicative factors, k ± + i m δ ± , appear at higher loop orders). The regulators are related by When working with ∆ ± as arising in the full theory, the soft function can be written as [113] ln S EIS However, the Q 2 dependence appearing in S here is artificial, as it is the δ ± regulators that are the fundamental Wilson line regulators. Using Eq. (B.19) together withp + p − = Q 2 , one can equivalently write this as where y n is an arbitrary parameter that cancels between the two logarithms. Eq. (B.21) allows one to split the soft function into a "left-moving" and "right-moving" component, 7 and it is this form of S EIS q that was used in NNLO calculations in Refs. [140,141,152]. The two factors in Eq. (B.23) are then combined with the n-collinear andn-collinear matrix elements, respectively. Since with the δ regulator the soft zero-bin subtractions on the two beam functions are identical to dividing by the original soft function, the soft subtraction amounts to dividing by √ S. For the n-collinear case, one identifies where B EIS q ≡ J n appears in Eq. (12) in Ref. [113]. Once again, only the combination arises in B q , such that the functional dependence of B q and J n in Eq. (B.24) agree.

One-loop results
The unsubtracted beam function has first been calculated in Ref. [112], where the δ regulator also acts as an IR regulator. The corresponding result with dimensional regularization as the IR regulator can be extracted from Ref. [141], Eq. (B.23) might seem to differ from Eq. (16) in Ref. [113], where the soft function is expressed as Since S only depends on the product of the last two arguments, one can identify α = e 2Y −2yn and use that α The one-loop soft function is given by [113] S EIS (1) According to Eqs. (B.24) and (B.1), the TMDPDF at NLO follows as

B.3 η regulator
The η regulator [115,123] modifies collinear and soft Wilson lines in momentum space as where the momentum operator P picks up the momentum of the gluon fields A, and η is the rapidity regulator with an associated rapidity scale ν. While the authors of Ref. [115] separately renormalize the beam and soft functions in η, one can trivially define a TMDPDF in the form of Eq. (B.1) by combining the bare beam and soft functions, in which case the dependence on η and ν cancels, where B q ≡ B CJNR q and S q ≡ S q CJNR are the beam and soft function as defined in Ref. [115]. For this regulator, the soft zero-bin subtractions vanish, so ∆ q S = S q CJNR only involves a √ S q in the numerator. For this regulator the correspondence with our notation is The choice of fixing y n = 0 arises because of the symmetric treatment of the two beam functions, but can be relaxed as in the other definitions if so desired. In Eq. (B.31), we combined bare beam and soft functions as usual for TMDPDFs. However, in practice the η regulator is often used to define rapidity-renormalized beam and soft functions by first expanding them individually about η → 0 and → 0 and then absorbing poles in η and through separate MS counterterms [115], When combined these counterterms give Using the fact the bare B CJNR q and S q CJNR are µ and ν independent yields the µ RGEs and the rapidity RGEs where the ν independence cancels on the right hand side, and we remind the reader that here ζ = (xP − ) 2 . Since B CJNR q only depends on the ratio xP − /ν, the Collins-Soper kernel can be extracted as This also shows that the rapidity anomalous dimension of Ref. [115] is related to the Collins-Soper kernel by γ q ν (µ, b T ) = 2γ q ζ (µ, b T ).

One-loop results
The bare quark-quark matrix element for the beam function is not explicitly given in the literature, but can easily be calculated analogous to the gluon beam function in Ref. [115]. We find The NLO soft function using this regulator can be found in Refs. [115,139], Combining the bare functions as in Eq. (B.31) gives the NLO TMDPDF as Alternatively, one can separately UV and rapidity renormalize beam and soft functions by absorbing the poles in η and in separate counterterms. The renormalized functions are given by

B.4 Exponential regulator
The exponential regulator [125] has been designed to connect the TMD soft function to the threshold soft function, which enabled it to be calculated up to O(α 3 s ) [142]. The regulator introduces a factor exp −k 0 τ e −γ E (B.44) into the phase space of real emissions. One then takes the τ → 0 limit, keeping only divergent terms. Identifying τ = ν −1 , one can connect this to the η regulator shown above.
Here, we will keep τ for clarity. Similar to the η regulator, in the exponential regulator one typically combines UV-and rapidity-renormalized beam and soft functions. A key difference is that the exponential regulator requires a zero-bin subtraction, which is equivalent to dividing by the soft function. Thus identifying one obtains the TMDPDF. Similar to our discussion for the η regulator this can be done either by first combining and then renormalizing the bare matrix elements, or by combining renormalized beam and soft functions.

One-loop results
The beam function is only given after UV renormalization in Ref. [153], The renormalized NLO soft function in this scheme is given by [125] S (1) Combining the renormalized functions gives the NLO TMDPDF as which reproduces Eq. (B.5) with ζ = Q 2 .

B.5 Analytic regulator
The analytic regulator as introduced in Ref. [109] for TMDs raises propagators including n andn-collinear momenta to a power, 8 n−collinear : where α and β are distinct parameters. One then has to consistently expand first in β → 0 and then in α → 0, or vice versa. Since this breaks the n ↔n symmetry, it regulates all rapidity divergences. However, the regulator does not act symmetrically in the n andncollinear beam function and thus yields two distinct bare quantities B BN q n and B BN qn . Following Ref. [109], we give all results for expanding first in β → 0 and then in α → 0, such that rapidity divergences are regulated by α alone.
A particular feature of this regulator is that the soft function is unity all orders, as all loop corrections are scaleless and vanish. This implies that rapidity divergences, i.e. poles in 1/α, cannot cancel between a single beam function and the soft factor, and hence one cannot define a rapidity-finite TMDPDF as in Eq. (B.1). Instead, divergences in α only cancel in the product, so one has to define the TMDPDFs through the limit is the so-called collinear anomaly coefficient [109], which exponentiates the rapidity logarithms. Note that the TMDPDFs on the right hand side of Eq. (B.52) do not have an explicit ζ dependence, as ζ a = ζ b = Q has already been chosen to fully exponentiate all rapidity logarithms.

One-loop results
The one-loop beam functions have been calculated in [109]. As discussed, the n andncollinear functions take different forms due to the regulator breaking the n ↔n symmetry.
The one-loop result for the matrix element in an on-shell quark state can be extracted from Ref. [109] in the form Only the product of the two beam functions is meaningful, as the poles in α and the ν dependence cancel, and one obtains the product of the two TMDPDFs, One can easily check that the same is obtained for f TMD

C Calculation of the quasi beam function
We calculate the spin-averaged quasi beam function by evaluating Eq. (3.19) with an onshell quark with momentum P µ = (P z , 0, 0, P z ), and carry out the Fourier transform to momentum space only for the z momentum, while the transverse component is kept in position space, where for compactness we suppress the dependence ofq n on and L. Here λ = 0 or λ = 3, and the path ordered exponential represents the Wilson line stretching along the path γ, as illustrated in Fig. 5. The explicit parameterization is given by three segments In Feynman gauge, there are four topologies contributing to the one-loop beam function, shown in Fig. 11. The sail and wave function diagram each have a mirror diagram. Since we work in pure dimensional regularization with on-shell quarks, the wave function renormalization diagram is scaleless and vanishes. In practice, it converts IR poles into UV poles, which is crucial for the renormalization of the beam function. The general strategy of the NLO calculation is as follows. First, we write down the integral fully in momentum space and introduce Feynman parameters. The resulting integral over the loop momentum q can then be reduced to a set of master integrals, defined as Here, K n (x) is the modified Bessel function of the second kind. We need only the following integrals over the Feynman parameter y, where µ 0 is the MS scale and µ is the MS scale, and x is the momentum fraction carried by the parton, not a Feynman parameter. Finally, we introduce plus distributions taking into account that prior to taking the physical limit, the quasi beam function has support x ∈ [−∞, ∞], rather than the physical support x ∈ [0, 1]. For the calculation, it is most convenient to define the plus distribution such that the integral over all x vanishes, For functions with physical support x ∈ [0, 1], this naturally reduces to the standard plus distribution. Any integrable expression h(x) can be written as a distribution using

C.1 Vertex correction
The vertex correction Fig. 11a is given bỹ (C.12) Note that the second term in square brackets vanishes for λ = 0, as the integrand is odd in k 0 . The Fourier transform w.r.t. b z can be carried out trivially, Introducing Feynman parameters, we obtaiñ (C.14) Using Eqs. (C.6) and (C.7), this can be expressed in terms of the master integrals, Next we use Eq. (C.11) to rewrite Eq. (C.15) in terms of a plus distribution. Note that the integral of the second line in Eq. (C.14) vanishes when integrated over all x and thus does not yield a boundary term. We obtaiñ Plugging in the master formulas and expanding in gives the final result where we singled out the terms with physical support x ∈ [0, 1]. These will be crucial to recover the collinear singularity of the standard PDF. The contribution with unphysical support x ∈ [−∞, ∞] is given by The first line in Eq. (C.18) is clearly suppressed for b T P z 1. The second line in Eq. (C.18) also vanishes in this limit, because the incomplete Gamma function behaves as 19) and the logarithmic enhancement for small arguments is compensated by the prefactor. The third line of Eq. (C.18) yields a power-suppressed contribution when convoluted with the TMDPDF, see appendix C.4.

(C.22)
For both λ = 0 and λ = z, the numerator can be simplified to P λ k µ − (P · k)g λµ + P µ k λ = P z k µ + k z g µ0 − k 0 g µz → P z (k µ − k 0 g µz ) , (C. 23) where the last step employs that the g µ0 piece vanishes when contracted with the spatial path γ (s) µ . Introducing Feynman parameters and dropping the terms linear in q 0 gives q (b) n (b µ ) = −g 2 C F µ 2 0 P z accordingly. Note that the second line in Eq. (C.26) only involves γ (s) z . All contributions from the transverse gauge link are thus fully captured in the first line, which becomes Taking the Fourier transform of this and pulling it inside a plus distribution yields q (b,1) n (x, b T , P z ) = 2ig 2 C F µ 2 0 P z I a 0 ±∞ + + 2ig 2 C F µ 2 0 δ(1 − x)I 2 . (C.28) The second contribution to Eq. (C.26) is more involved. We first simplify the exponentials by letting y → 1 − y and then q → q + yP , The first term in brackets has a singularity at q z = 0. However, when integrating over x, the singularity cancels with the second line, so it can be regulated by introducing the plus distribution prior to evaluating the q integration, For x = 1, this can be expressed using our master integrals, The δ(1 − x) term is fixed by integratingq (b,2) n (x, b T , P z ) over x, where we have singled out the terms with unphysical support x ∈ [−∞, ∞], (C.38)

C.3 Wilson line self energy
The general expression for the Wilson line self energy, Fig. 11c, in position space is given byq Diagrams (a) and (b) are independent of b z , so the Fourier transform w.r.t. b z is trivial and one can directly evaluate the line integrals similar to the soft function calculation shown in Sec. 4.3, obtaining Diagram (c) depends on b z , so care has to be taken with the Fourier transform. The position space result is The Fourier transform can be obtained using [45] db z 2π As shown in appendix C.4, when convoluted with the TMDPDF the phase term will be power suppressed. Note that if one were to expand Eq. (C.42) in b z L prior to Fourier transforming, one would instead obtaiñ q (c, 3) n (x, b T , P z ) b z L = db z 2π e ib z xP z α s C F 2π P z e −ib z P z 1 + ln so dropping the phase in Eq. (C.44) is equivalent to performing a small-distance expansion b z L in position space.
Diagram (d) gives rise to a factor 2 to lift the symmetry factor. Due to the finite separation b T , it is UV finite, so we can let → 0. In position space, we obtaiñ . (C.46) Here, it is quite difficult to take the Fourier transform while keeping th exact b z dependence. It is easier to use the first line of Eq. (C.39) and leave the k z integration until the end, q (c,4) n (x, b T , P z ) = ig 2 C F P z µ 2 0 db z 2π e i(x−1)P z b z L(L − b z ) While the second line is manifestly finite as k z → 0, the last line has an apparent singularity, which is as usual treated as a plus distribution, For comparison, the b z L limit of Eq. (C.46) is given bỹ .

(C.49)
This still has a nontrivial b z dependence, so in order to obtain the δ(1−x) term in Eq. (C. 48) we further need to expand in b z b T , giving (C.50) Hence dropping the phase term in Eq. (C.48) corresponds to a small-distance expansion b z L and b z b T .
Combined result. Combining Eqs. (C.40), (C.44) and (C.48), we obtain the full exact tadpole diagram as where the terms suppressed for b z b T L are

C.4 Power-suppressed contributions to the matching kernel
The unphysical contributions Eqs. (C.18), (C.38) and (C.52) contain fastly oscillating phases ∼ e −iP z L(1−x) , and one may thus expect that these give vanishing contributions as LP z → ∞. However, these phases vanish as x → 1, and furthermore can be associated with divergences in 1/(1 − x), so they can contribute nontrivially to the quasi beam function. Thus, in order to neglect them, one has to show that they do not contribute to the matched TMDPDF in the limit LP z → ∞. The unphysical terms cannot yield a multiplicative matching with the TMDPDF, as the support of quasi TMD and TMD do not match. Given our discussion in Sec. 3.1, this provides a direct indication for the fact that they are power suppressed, and we will show that this is consistent. Assuming that their contributions satisfy a convolution structure as in Eq. (3.2), then the unphysical part of the matching is given by ∆q(x, b T , P z ) = ∆q (a) n (x, b T , P z ) + ∆q (b) n (x, b T , P z ) + ∆q (c) n (x, b T , P z ) . (C. 53) In the following, we suppress all arguments except x and y as well as all superscripts and subscripts and the flavor indices for brevity. We also extend the integral in Eq. (C.53) to infinity by implicitly assuming that f TMD (y) = 0 for |y| ≥ 1. Changing the integration variable in Eq. (C.53) from y to x/y, we obtain where we can drop the plus prescription from now on, as the term in square bracket cancels an overall 1/(1 − y) divergence. Let us now consider the case where ∆q is given by an exponential phase factor, a regular function r(y) and potentially a divergence as y → 1, ∆q(y)