Renormalization and Matching for the Collins-Soper Kernel from Lattice QCD

The Collins-Soper kernel, which governs the energy evolution of transversemomentum dependent parton distribution functions (TMDPDFs), is required to accurately predict Drell-Yan like processes at small transverse momentum, and is a key ingredient for extracting TMDPDFs from experiment. Earlier we proposed a method to calculate this kernel from ratios of the so-called quasi-TMDPDFs determined with lattice QCD, which are defined as hadronic matrix elements of staple-shaped Euclidean Wilson line operators. Here we provide the one-loop renormalization of these operators in a regularizationindependent momentum subtraction (RI′/MOM) scheme, as well as the conversion factor from the RI′/MOM-renormalized quasi-TMDPDF to the MS scheme. We also propose a procedure for calculating the Collins-Soper kernel directly from position space correlators, which simplifies the lattice determination. ar X iv :1 91 0. 08 56 9v 1 [ he pph ] 1 8 O ct 2 01 9


Introduction
Transverse-momentum dependent parton distribution functions (TMDPDFs) describe the longitudinal and transverse momentum distribution of quarks and gluons in hadrons and nuclei, and thus are of vital interest to improving our understanding of hadronic and nuclear structure [1,2]. They are also crucial to predicting transverse momentum distributions in the Drell-Yan process, a key observable both for the Tevatron [3][4][5][6] and the LHC [7][8][9][10][11][12], as well as in semi-inclusive deep-inelastic scattering at low energies [13][14][15][16][17][18][19]. TMDPDFs measure the transverse momentum q T carried by the struck parton. For perturbative q T Λ QCD , they can be calculated in terms of collinear parton distribution functions, and the resulting matching formula is known to next-to-next-to-leading order (NNLO) [20][21][22][23][24][25][26][27]. In contrast, for nonperturbative q T Λ QCD , TMDPDFs become genuinely nonperturbative objects which so far have only been extracted from measurements by performing global fits to a variety of experimental data sets, see e.g. Refs. [28][29][30][31][32][33]. Since there are some issues associated to these extractions, in particular with reconciling low and high energy data, an independent determination from first principles is highly desirable. This has motivated studies with lattice QCD that have been carried out in Refs. [34][35][36][37][38], primarily for ratios of moments in the longitudinal momentum fraction.
The TMDPDF f i (x, b T , µ, ζ) for a parton of flavor i depends on the longitudinal momentum fraction x and the position space parameter b T , which is Fourier-conjugate to q T , as well as the renormalization scale µ and the Collins-Soper scale ζ [39,40]. The latter encodes the energy dependence of the TMDPDF, i.e. the momentum of the hadron or equivalently the hard scale of the scattering process, and the associated evolution is governed by the Collins-Soper equation [39,40] ζ d dζ The γ i ζ (µ, b T ) that appears here is referred to as either the Collins-Soper (CS) kernel or rapidity anomalous dimension for the TMDPDF. From consistency of the ζ and µ evolution equations, combined with information on the all order structure of the renormalization for Wilson line operators, it is known to have the all-order structure where Γ i cusp (α s ) is the cusp anomalous dimension and γ i ζ (α s ) is the noncusp anomalous dimension, both of which are known perturbatively in QCD at three loops, see Refs. [41][42][43] and Refs. [22][23][24][44][45][46][47][48][49][50], respectively. 1 As should be evident from eq. (1.2), the Collins-Soper kernel becomes genuinely nonperturbative when b −1 T Λ QCD , independent of the renormalization scale µ. Consequently, the scale evolution of TMDPDFs becomes nonperturbative itself, and relating TMDPDFs at different energies requires nonperturbative knowledge of γ i ζ (µ, b T ), even if one chooses a perturbative µ ≥ 1 GeV. When extracting TMDPDFs from global fits, it is thus also necessary to fit γ i ζ , which is typically achieved by splitting the kernel into a perturbative and nonperturbative piece, where µ b ≡ µ b (µ, b T ) is chosen to always be a perturbative scale, such that all nonperturbative physics is separated into the function g i (b T , µ b ). A common parameterization of g i (b T , µ b ) is to assume a Gaussian form, g i (b T , µ b ) = g i K b 2 T with constant g i K [32,33], which has also been motivated by a renormalon analysis [61], but other forms have also been employed [62,63]. In the literature, there is a considerable discrepancy between Refs. [33,62] and Ref. [32] on whether the nonperturbative part of γ i ζ is crucial to describe the measured data or not. This is perhaps not surprising, as Refs. [33,62] are based on Drell-Yan data at relatively large q T , where one expects nonperturbative effects to be suppressed, while they become much more important in the lower energy measurements included in Ref. [32].
The lack of precise knowledge of the nonperturbative part of γ i ζ (µ, b T ) from global fits motivates an independent determination from lattice QCD. Here, a key difficulty is that TMDPDFs are defined as lightcone correlation functions which depend on the Minkowski time, while first principles lattice QCD calculations are inherently restricted to the study of Euclidean time operators. Large-momentum effective theory (LaMET) was proposed to overcome this hurdle in a systematically improvable manner for collinear PDFs (and generalized PDFs) by relating so-called quasi-PDFs, defined as equal-time correlators, through a perturbative matching to the physical PDF [64,65]. For these collinear quasi-PDFs, significant progress has been made, in particular on their renormalization and matching onto PDFs , and the study of power corrections to the matching relation [95][96][97], and first lattice calculations of the x-dependence of PDFs and distribution amplitudes have been carried out in Refs. [83,95,[98][99][100][101][102][103][104][105][106][107][108][109][110][111][112][113][114][115]. Recent lattice calculations at the physical pion mass have shown encouraging results for a precise determination of PDFs using the LaMET, including in particular those of the European Twisted Mass Collaboration [106,109], and results reported by the Lattice Parton Physics Project Collaboration [107,111,113].
The application of LaMET to obtain TMDPDFs from lattice has only been studied very recently [116][117][118][119]. A key difference to collinear PDFs is the necessity to combine a hadronic matrix element with a soft vacuum matrix element in order to obtain a well-defined (quasi) TMDPDF. In Ref. [119] it was shown that this soft factor, which involves lightlike Wilson lines, can not be simply related to an equal-time correlation function computable on the lattice, and hence without a careful construction of the quasi-TMDPDF one can only generically expect to encounter a nonperturbative relation between TMDPDFs and quasi-TMDPDFs, rather than a relation that is determined by a perturbatively calculable short distance coefficient. 2 However, in certain ratios of TMDPDFs this soft factor, and physically related contributions in the TMD proton matrix elements, cancel out. Hence such ratios can be obtained from ratios of suitably defined quasi-TMDPDFs which can be obtained from lattice. In particular, Ref. [118] showed that the Collins-Soper kernel can be obtained from such a ratio using 3 where C ns is a perturbative matching coefficient given at one loop in Refs. [118,119],f ns is the nonsinglet (u−d) quasi-TMDPDF, and P z 1 = P z 2 are two different proton momenta that are used for the corresponding quasi-TMDPDF calculations.
In order to obtain γ i ζ in the MS scheme, in eq. (1.4) the quasi-TMDPDFf ns is assumed to be in the MS scheme as well. Consequently, a critical step in this approach is the renormalization of the bare quasi-TMDPDF on the lattice and its subsequent scheme conversion into the MS scheme. On the lattice, the renormalization should be performed in a scheme that is defined nonperturbatively to facilitate the removal of both linear and logarithmic divergences. In contrast the scheme conversion factor into MS is calculated order-by-order in continuum perturbation theory, and does not depend on the ultraviolet (UV) regulator since it is simply a difference of renormalized quantities. For the longitudinal quasi-PDFs, such nonperturbative renormalization, scheme conversions, and the associated matching to obtain the analog of C ns have been studied and implemented in Refs. [78,86,102,103,110] with the regularization-independent momentum subtraction (RI/MOM) schemes [120]. Such a calculation has also been carried out in Ref. [121] for staple-shaped Wilson line operators at vanishing longitudinal separation, which is connected to the calculations needed for determining TMDPDFs. In particular it corresponds to a special case of the quasi-TMDPDF operators studied here, which will involve stapleshaped Wilson lines but with an additional separation along the longitudinal direction. In this paper, we determine the scheme conversion coefficient between the RI /MOM scheme 4 and MS for quasi-TMDPDFs, including the longitudinal separation, and also calculate the corresponding one-loop matching coefficient C ns . This paper is structured as follows. In Sec. 2 we briefly review the definition of (quasi) TMDPDFs and how the Collins-Soper kernel γ q ζ can be extracted from lattice Refs. [118,119]. In Sec. 2.4 we also propose a new improved method for obtaining γ q ζ which reduces systematic uncertainties in the lattice analysis by directly exploiting the quasi-TMPDF correlators in longitudinal position space. We then proceed in Sec. 3 to discuss the general structure of the RI /MOM renormalization and scheme conversion from the RI /MOM scheme to the MS scheme, before giving details on our one-loop calculation of the required renormalization and conversion factors in Sec. 4. The impact of these results are numerically illustrated in Sec. 5, before concluding in Sec. 6. In appendix A we collect formulae for the master integrals used in Sec. 4.

Determination of the Collins-Soper kernel from lattice QCD
In this section we briefly review the definition of TMDPDFs and the construction of quasi-TMDPDFs computable on lattice, as well as how the Collins-Soper kernel can be determined from these, and refer to Refs. [118,119] for more details. We also show how to determine the Collins-Soper kernel directly in position space, which is better suited to a lattice calculation than the method proposed in Refs. [118] to obtain the kernel in momentum space.

Definition of TMDPDFs
We define the quark TMDPDF for a hadron moving close to the n µ = (1, 0, 0, 1) direction with momentum P µ as where we use the lightcone coordinates b ± = b 0 ∓ b z and b T are the Euclidean transverse coordinates. In eq. (2.1), the bare beam function B q is a hadronic matrix element encoding collinear radiation, and the bare soft factor ∆ q S is constructed from a soft vacuum matrix element, to be defined below. The TMDPDF gives the probability to obtain a quark with lightcone momentum p − = xP − and transverse momentum q T , which is Fourier-conjugate to the parameter b T . Z q uv is the UV renormalization constant, with being the UV regulator and µ the associated renormalization scale. Beam and soft functions individually suffer from so-called rapidity divergences [40,47,[122][123][124][125][126][127], which are regulated by an additional regulator denoted as τ , and these divergences give rise to the Collins-Soper scale ζ. However, the rapidity divergences cancel between beam and soft function as τ → 0, giving rise to a well-defined TMDPDF. For a detailed discussion of different rapidity regularization schemes, see e.g. Ref. [119].
The bare quark beam function is defined as where [...] τ denotes the rapidity regularization of the operator, h(P ) denotes the hadron state of momentum P µ , the quark fields are separated by b µ = b +nµ /2 + b µ T withn µ = (1, 0, 0, −1), and the Wilson lines are defined as 5 3) The bare quark soft function is defined as where as before [...] τ denotes the rapidity regularization, and the Wilson lines are given by The Wilson line paths of both beam and soft function are illustrated in Fig. 1. Finally, the soft factor ∆ q S entering eq. (2.1) is defined as

6)
5 Note that we have changed the sign of the strong coupling g compared to Refs. [118,119] to agree with the convention that the covariant derivative is given by D µ = ∂ µ + igA µ . This sign agrees with Ref. [132] which we use as a reference for Euclidean Feynman rules for our calculation.  where S q is the soft function defined in eq. (2.4) and S q 0 is a subtraction factor necessary to avoid double counting of soft physics in the beam and soft function. Its definition depends on the employed rapidity regulator τ , but as the notation indicates, it is typically closely related to S q itself. For example, in the scheme of Ref. [127] one has S q 0 = 1, while in the schemes of Ref. [48,125] one has S q 0 = S q . For more details, see Ref. [119].

Definition of quasi-TMDPDFs
The quasi-TMDPDF is defined analogous to eq. (2.1), but as an equal-time correlator rather than a lightcone correlation function, namelỹ HereB q is the quasi-beam function,∆ q S includes the quasi-soft function together with subtractions,Z q uv carries out UV renormalization in a lattice-friendly scheme, whereμ stands for any added scales introduced by this scheme choice, andZ q converts the result perturbatively to the MS scheme with scale µ. The quasi beam and soft functions will be constructed such that all Wilson line self energies proportional to L/a and b T /a, as well as divergences ∝ L/b T which correspond to rapidity divergences in the lightlike case [117,119], cancel betweenB q and∆ q S . Therefore, the remaining UV renormalizationZ q uv and the scheme conversionZ q only depend on b z , but not necessarily b T or L. 6 The bare quasi-beam function is defined as , and the UV regulator is denoted as a, following the notation for the finite lattice spacing that acts as a UV regulator in lattice calculations. Due to the finite lattice size, the longitudinal Wilson lines are truncated at a length L less than the size of the lattice, which also regulates the analog of rapidity divergences [117,119]. Compared to eq. (2.2), we also replaced γ − by the Dirac structure Γ, which can be chosen as Γ = γ 0 or Γ = γ z . (Technically, one can also use a combination, for example γ 0 + γ z .) The Wilson lines of finite length L are defined by while the transverse gauge links are identical to those in eq. (2.3).
For the quasi soft function, we use the bent soft function of Ref. [119], defined as Here,n µ ⊥ is the transverse unit vector orthogonal to n µ , then n µ ⊥ = (0, cos φ, sin φ, 0), and a valid choice for n µ ⊥ isn µ ⊥ = (0, − sin φ, cos φ, 0). The Wilson lines in eq. (2.10) are defined as The Wilson line paths of both quasi beam and quasi soft function are illustrated in Fig. 2 for the choice φ = 0. The quasi-soft factor is obtained from the bent soft function as whereS q 0 =S q bent is the subtraction factor which avoids double counting between quasi beam and soft functions. The overall length of the Wilson lines appearing in∆ q S must be chosen to ensure the cancellation of Wilson line self energies in eq. (2.7) [119], whereas implementing this with the specific choiceS q 0 =S q bent corresponds to a particular scheme.
Note that for the construction of the quasi-TMDPDF, different definitions of the quasi soft function could be employed as well. This yields different definitions of the quasi-TMDPDF, which will affect the (possibly nonperturbative) kernel relating quasi-TMDPDFs and TMDPDFs, see Ref. [119] for a more detailed discussion. With the bent soft function in eq. (2.4), this relation was shown to be short distance dominated and hence perturbative at one loop, which motivates its use here. Importantly, for the determination of the Collins-Soper kernel the soft factor always cancels, such that this precise definition does not matter.
The spacelike Wilson lines ofB q as given in eq. (2.8) and those ofS q as given in eq. (2.10) give rise to self energies that yield power law divergences proportional to e δm Ltot . Here, δm is a mass correction that absorbs divergences as a → 0, and the total lengths of the Wilson line structures are given by L B tot = L + |L − b z | + b T forB q and L S tot = 4L + 2b T forS q , respectively. After combining the quasi beam function with the square root of the quasi soft function, the Wilson line self-energies yield the overall power-law divergence which has to be absorbed byZ q uv (b z ,μ, a). To cancel this divergence on the lattice, the nonperturbative UV renormalization has to be applied before the Fourier transform, as shown in eq. (2.7), while in the lightlike case it is independent of b z and can be pulled out, see eq. (2.1). This distinction is important, implying in the ratio of TMDPDFs the UV renormalization factor Z q uv cancel out, whereas this is not possible for ratios of quasi-TMDPDFs.

Determination of the Collins-Soper kernel in momentum space
In this section, we briefly review the method proposed in Ref. [118] for calculating the Collins-Soper kernel from lattice QCD. As discussed in Ref. [118], and in more detail in Ref. [119], in general there is a mismatch between the infrared structure of the quasi-beam function and beam function due to the fact that the latter requires a dedicated rapidity regulator, whereas the former has rapidity divergences regulated by the finite length L. This spoils the simplest boost picture from LaMET (even when supplemented by short distance corrections), for relating these proton matrix elements. Nevertheless, when combined with the quasi-soft and soft functions, these divergences and the L dependence cancel, enabling the possibility of a matching equation between the quasi-TMDPDF and TMDPDF. However, even after these cancellations there can still be a mismatch in the remaining infrared structure of the quasi-soft and soft functions, leaving a relation of the form (2.14) Here, C ns is a perturbative kernel for the nonsinglet ns=u−d channel, g S q is a nonperturbative contribution which reflects the mismatch in soft physics, and γ q ζ is the standard Collins-Soper kernel, which allows one to relate the TMDPDF at the scale ζ to the quasi-TMDPDF at proton momentum P z . Corrections to this matching relation are suppressed for large L and P z , as indicated, and will be suppressed in the following. In Ref. [119] it was shown that the bent soft function yields To demonstrate that eq. (2.14) is a true matching equation requires an all-order proof that g S q (b T , µ) = 1, which has not been demonstrated. However the lack of this proof does not impact the determination of the anomalous dimension γ q ζ (µ, b T ), to which we now turn. Evaluating eq. (2.14) at two different proton momenta P z 1 = P z 2 but the same ζ, and taking the ratio of the results yields Here g S q and f ns have dropped out. In Ref. [118], this was solved for γ q ζ as .
Note that here we have canceled the quasi soft factor∆ q S (b T , a, L) in the ratio, as it is independent of b z . The advantage of doing so is that one needs to calculate one less nonperturbative function from lattice QCD. The price to pay is thatB ns still contains Wilson line self energies ∝ L/a, b T /a and divergences ∝ L/b T , which now only cancel in the ratio rather than in the numerator and denominator, respectively. To achieve the separate cancellation, we can simultaneously insert a b z -independent factorR B in both the numerator and denominator to separately cancel these leftover divergences, .
This factor has to be constructed such that it exactly removes all divergences that would normally be canceled by∆ q S (b T , a, L), i.e. all power-law divergences not yet absorbed bỹ Z q uv (b z ,μ, a). One trivial choice for this factor is thus to use the soft factor that was canceled before,R B =∆ q S , while another simple choice would beR B = [B ns (0, b T , a, P z R , L)] −1 , i.e. the quasi beam function at vanishing separation b z = 0 and some reference momentum P z R . In Sec. 3, we will construct a more refined expression by using the nonperturbative RI /MOM renormalization factor in a similar fashion, and the final definition for the com-binationZ q uvRB will be given in eq. (3.12). As stressed in Ref. [118], eqs. (2.17) and (2.18) are formally independent of x, P z 1 and P z 2 , up to power corrections as indicated in eq. (2.14), such that one can use any residual dependence of the lattice results on these parameters to assess systematic uncertainties.
The one-loop result for the matching coefficient C ns that enters eq. (2.17), whenf ns , Z q and γ q ζ are in the MS scheme, has been calculated in Refs. [118,119] and is given by This short distance coefficient can also be extracted from the results of Ref. [117]. Note that C ns is an even function of its second argument, C ns (µ, −xP z ) = C ns (µ, xP z ).

Determination of the Collins-Soper kernel in position space
A potential drawback of using eq. (2.17) and eq. (2.18) is that one has to Fourier transform the position-space correlatorB ns (b z , b T , a, P z , L). This can be a limiting factor due to the finite number of b z values available from lattice, which makes the numerical Fourier transform quite challenging. We hence propose in this section a related but modified formula which enables the matching to be performed directly in position space. This will turn the Fourier transform into a numerically more stable convolution.
To derive this relation, we need the Fourier transforms of the quasi TMDPDF, Note that for simplicity we distinguish the quasi-TMDPDF in position and momentum space only by their arguments, as it is always clear from context which one we refer to.
It will be convenient to work with the Fourier transform of the inverse of the kernel C ns , defined throughC Plugging eqs. (2.20) and (2.21) back into eq. (2.15), we get Next, we Fourier transform both sides from momentum fraction x to a dimensionless position y, by multiplying by e −ixy and integrating over x, obtaining This can trivially be solved for γ q ζ as .

(2.24)
As expected, Fourier transforming the product in eq. (2.15) yields a convolution in position space. In eq. (2.24),f ns is the renormalized nonsinglet quasi-TMDPDF in position-space as calculated on lattice. Using the expression eq. (2.7) forf ns , we obtain the final expression where we suppress the explicit limits L → ∞ and a → 0 for simplicity. As in eq. (2.18), we have inserted a factorR B that cancels all divergences in L/a, b T /a and L/b T separately in the numerator and denominator, which otherwise would only cancel in the ratio. Again, formally the dependence of the right hand side of eq. (2.25) on y, P z 1 and P z 2 cancels up to power corrections, such that one can use any residual dependence of the lattice results on these parameters to assess systematic uncertainties. To use the improved formula in eq. (2.25) one only needs the position-space proton matrix elementB ns (directly obtained on the lattice), its renormalization factorZ q uv (also obtained on the lattice) combined with the factorR B (discussed in Sec. 3), the MS-conversion factorZ q (which we calculate in Secs. 3 and 4 of this paper), and the Fourier-transformed matching kernelC ns (which we obtain below).
Note that we have chosen the definition eq. (2.21) of the position-space kernelC ns to be determined by the transform of the inverse of C ns in order to make eq. (2.25) particularly simple, with a numerator depending only on the momentum P z 1 and the denominator only on P z 2 . For comparison and completeness, we present in appendix B the corresponding derivation when using a position-space kernelC ns that is defined by the transform of C ns itself, in which case numerator and denominator would both depend on P z 1 and P z 2 . The Fourier transformC ns can be further simplified by employing that in the physical limit L, P z → ∞,f q has limited support x ∈ [0, 1] for quarks and x ∈ [−1, 0] for antiquarks. Hence, we can make different choices for the integration range in eq. (2.21) which lead to formally equivalent results when the resulting coefficientsC ns are employed in eq. (2.25).
To exploit this freedom we consider the two natural choices, defininḡ where the superscript inC D ns denotes the integration domain D. Physically,C Note that since C ns depends logarithmically on xP z /µ, its Fourier transform according to eq. (2.21) with unconstrained integration range will involve plus distributions which are complicated to implement numerically, so here we will refrain from advocating for using the unrestricted integration, and hence only present the simpler results obtained using eqs. (2.26) and (2.27).
Matching kernel in position space. Next we explicitly calculate the Fourier transform of C ns to position space as defined in eqs. (2.26) and (2.27). C ns was calculated at nextto-leading order (NLO) in the MS scheme in Refs. [118,119] and is given in eq. (2.19). Perturbatively inverting it gives the one-loop result where we abbreviated L z = ln[(2P z /µ) 2 ], and as before the superscript D on the three required functions f D i (y) denotes the integration range of the Fourier transform. For the case of integrating over D = [0, 1] the auxiliary integrals are (2.30) Here, n F n is a hypergeometric function. The results for n = 0 and n = 1 can be expressed using standard functions,

RI /MOM renormalization and matching
The determination of γ q ζ using either eq. (2.18) or (2.25) requires calculating the quasibeam functionB q from lattice, a renormalization of UV divergences withZ q uv , a definition ofR B to cancel remaining power-law divergences (one choice would be∆ q S ), and finallỹ Z q to convert to the MS scheme. Here, we specify in detail a preferred choice for how to construct these nonperturbative renormalization factors in the RI /MOM scheme, and how the conversion factorZ q can be calculated perturbatively.Z q is then calculated at one loop in Sec. 4.
Note that for b z = 0, the corresponding MS conversion kernel for the quasi beam functionB q has been calculated in Ref. [121], which is sufficient for the lattice studies of the x-moments of TMDPDFs carried out in Refs. [34][35][36][37][38], but does not suffice for the determination of the Collins-Soper kernel which requires the calculation for nonvanishing b z .
We first argue that in continuum theory, the staple-shaped quark Wilson line operator entering the quasi beam function can be renormalized multiplicatively in position space as Here, Wẑ and W T are Wilson lines as defined in eqs. (2.3) and (2.9), and for brevity we suppressed their explicit arguments, which are given in eq. (2.8). In the first line in eq. (3.1), we work with bare quark fields ψ 0 and Wilson lines built of bare gluon fields and bare couplings, while in the second line we work with renormalized fields and couplings, indicated by the subscript R. Z q,wf includes all the logarithmic UV divergences originating from the wave function renormalization and the quark-Wilson-line vertices. The exponential absorbs all linear power divergences from the self-energies of the spacelike Wilson lines, where L + |L − b z | + b T is the total length of the staple. In eq. (3.1), the dependence on the UV cutoff a is kept implicit on the right hand side. Eq. (3.1) resembles the multiplicative renormalization for the quasi-PDF operator, for which the renormalization in the RI /MOM scheme has been studied in Ref. [78,86,102,103,110]. For the staple-shaped operators discussed here, the multiplicative renormalization has also been used in Ref. [121], which carried out the RI /MOM renormalization for the special case b z = 0, i.e. vanishing longitudinal separation of the staple. However, as pointed out in Refs. [38,121], there will be mixing with other operators on a discretized lattice, and thus the renormalization on lattice requires an independent study [38,128]. In this work, we focus on the renormalization and scheme conversion in continuum theory, where the multiplicative renormalizability of the quasi beam function can be proven analogous to the proof for straight Wilson line operators by employing the auxiliary field formalism [81,83,92,129]. This auxiliary field formalism is also commonly used to derive Wilson line operators in the Soft Collinear Effective Theory, see the original work in Refs. [130,131]. By constructing three independent auxiliary "heavy quark" fields for each edge of the staple-shaped Wilson line, the nonlocal quark Wilson line operator can be reduced to the product of four composite operators in the effective theory that includes these auxiliary fields. BRST invariance implies that this effective theory is renormalizable through multiplicative counterterms to all orders in perturbation theory [81]. It follows that in continuum QCD, the staple-shaped quark Wilson line operator can indeed be renormalized multiplicatively as shown in eq. (3.1).
To implement the RI /MOM scheme for the quasi beam function, one first computes the amputated Green's function of the operator given in eq. (3.1), which is also referred to as the vertex function. Here and below b indicates dependence on b z and b T . In eq. (3.2), S 0 (p, a) is the bare quark propagator that can be calculated nonperturbatively on the lattice. Λ Γ 0 (b, a, p, L) is a linear combination of Dirac matrices that are allowed by the symmetries of space-time and the operator itself. For off-shell quarks, there will also be finite mixing with equation-of-motion operators that vanish in the on-shell limit.
In practice, one needs to choose a projection operator P to define the off-shell matrix element of the quasi-beam function from the amputated Green's function, The choice of P is not unique [86,110], but it must have overlap with Γ to project out all the UV divergences as a → 0. In Ref. [78,121], the choice is P = Γ, while in Refs. [86,110] both the choice P = / p, and a choice for P that effectively projects out the coefficient of Γ in the covariant decomposition of Λ Γ 0 , were considered. In the RI /MOM scheme, the renormalization constantZ Γ,P B of the bare operator O Γ 0 defined in eq. (3.1) is determined by requiring that at a specific momentum p µ R , the projection defined in eq. (3.3) reduces to its value at tree level in perturbation theory. Here, we actually need to define the RI /MOM condition for the quasi-TMDPDF, which also includes the soft factor. It reads Here, q Γ,P tree is the value of eq. (3.3) at tree-level in perturbation theory, which is nonzero only for particular choices of Γ and P, and each such pair (Γ, P) define a particularZ Γ,P q . The tree level soft factor is given by∆ q S tree = 1 and hence not explicitly given in eq. (3.4). Here the choice for the scales p R and b R T are part of the definition of the RI /MOM scheme. In eq. (3.4), the wave function renormalization factor Z wf arises to compensate for the renormalization of the bare quark fields in eq. (3.2). It is determined independently with the following condition on the quark propagator, where the 1/4 arises from the trace over Dirac indices. 7 The use of eq. (3.5) in eq. (3.4) defines the RI scheme, while in the closely related RI scheme Z wf is defined by imposing vector current conservation using Ward identities [120]. From eqs. (2.7) and (3.4), it follows that where we have splitZ Γ,P q into a pieceZ Γ,P B arising from the RI /MOM prescription applied to the quasi beam function only, and the quasi soft factor∆ q S (b R T , a, L). TheZ Γ,P B is given by the RI /MOM conditioñ . (3.7) From eq. (3.6) we can also identify the RI /MOM renormalization scaleμ = (b R T , p R ), which contains a choice for both the momentum p µ R and the transverse separation b R T to be used when defining the renormalization constant. This is unusual for a RI /MOM scheme, where one would normally only specify p µ R , but not b R T . The reason to also specify b T = b R T here is that b T itself can become a nonperturbative scale, and hence must not enter the perturbative scheme conversion factorZ q (b z , µ,μ). In contrast, b R T can always be chosen to be a perturbative scale, similar to p µ R , thus ensuring that this scheme conversion factor to MS remains perturbatively calculable.
UsingZ Γ,P q , the bare quasi-TMDPDF can be renormalized in position space as The RI /MOM-renormalized quasi-TMDPDF obtained from eq. (3.8) is independent of the UV regulator, and therefore can be matched perturbatively onto the MS renormalized quasi-TMDPDF, which is given bỹ . (3.10) Note that in eq. (3.8), all divergences in L/a, b T /a and L/b T cancel amongB q and∆ q S , rather than being absorbed byZ Γ,P q . However, for the determination of the Collins-Soper kernel γ q ζ as suggested in eqs. (2.18) and (2.25), it was advantageous to cancel out the soft factor in the ratio so it does not have to be calculated on the lattice. In this case these power law divergences also only cancel in the ratio. Such power law divergences can be problematic since it is generally unwise to attempt to extract a signal only after canceling out large contributions. This can be remedied by constructing the factorR B to precisely cancel these divergences. A convenient choice in the RI /MOM scheme is . (3.11) Hence the combination that enters eqs. (2.18) and (2.25) is given bỹ In this expression, the quasi soft factor∆ q S has canceled betweenZ Γ,P q andR B , as desired so that it does not need to be calculated on the lattice. The factorZ Γ,P B (0, b T , p R , a, L) in eq. (3.12) cancels all divergences in L/a, b T /a and L/b T that are present inB q (b z , b T , a, P z , L), while the remaining fraction in eq. (3.12) removes all leftover UV divergences, in particular those proportional to b z /a. Thus this result fulfills all requirements of eqs. (2.18) and (2.25).
As a result, both the numerator and denominator in eqs. (2.18) and (2.25) have well defined continuum (a → 0) and L → ∞ limits before one calculates their ratio. Nevertheless, if one loosens the requirement for convergence by taking the continuum (a → 0) and L → ∞ limits after calculating the ratio, then it is important to note that the combinatioñ Z Γ,P B (0, b T , p R , a)/Z Γ,P B (0, b R T , p R , a) will cancel out in the ratios in eqs. (2.18) and (2.25). If this is done, then the UV and L/b T divergences it accounts for will still cancel out between the numerator and denominator in those limits.

One-loop results
In this section, we provide details on the one-loop calculation of the quasi-beam function in an off-shell state, as well as the resulting RI /MOM renormalization factorZ q uv and the RI /MOM to MS conversion factorZ q for the quasi-TMDPDF. Throughout this section, we work in Euclidean space with d = 4 − 2 dimensions, as our aim is to calculateZ q with a lattice friendly definition of the Lorentz indices.

Quasi-beam function with an off-shell regulator
For completeness, we first give the Feynman rules in Euclidean space, following the notation of Ref. [132]. For a covariant gauge with gauge parameter ξ, the gluon propagator reads The Feynman rules for quark propagators and the QCD vertex are given by where the sign of the strong coupling constant g is such that the covariant derivative is given by D µ = ∂ µ + igA µ . In eqs. (4.1) and (4.2), k E is a Euclidean momentum such that k 2 The γ µ E in eq. (4.2) are Dirac matrices in Euclidean space, which are related to the Dirac matrices γ µ M in Minkowski space by γ 0 E = γ 0 M and γ i E = iγ i M , and obey In the remainder of this section, we suppress the explicit subscript "E", as we will always work in Euclidean space.
We consider the matrix element of the quasi TMD beam function operator in eq. (3.1) with an off-shell quark state |q s (p) of Euclidean momentum p 2 > 0, amputated to remove the spinors, The full set of possible projection operators is Note that only P 1 = γ ρ with ρ = λ yields a nonvanishing tree-level result and thus a valid renormalization. However, it is also interesting to study the mixing between different Dirac structures, and hence we also consider all other projectors eq. (4.4) yielding nonvanishing one-loop results. For our continuum analysis, this includes only the axial and vector projection operators, so we consider two different projections ofΛ λ ξ to define the off-shell matrix element of the quasi-beam function: Here, the subscript "a" refers to axial. We split all results into a piece corresponding to Feynman gauge (ξ = 1) plus a correction for ξ = 1, and similarly for the axial projection. Here, ξ = 0 corresponds to the Landau gauge most relevant for lattice. The tree level results are given bỹ At one loop, there are four topologies contributing toq ρλ ξ andq ρλ a,ξ , as shown in Fig. 3. To evaluate these, we introduce two master integrals, Explicit results for these are collected in appendix A. In eq. (4.8), µ 0 is the MS renormalization scale, which is related to the MS scale by In the following, we derive results for the different diagram topologies in terms of these master integrals, keeping the Dirac indices ρ and λ as well as the gauge parameter ξ generic. Note that only the sail diagram is nonvanishing for the axial projectionq ρλ a,ξ , and hence for the other diagrams we only discussq ρλ ξ .

Vertex diagram
The vertex diagram shown in Fig. 3a is given bỹ After evaluating the Dirac trace, the integrals can be expressed in terms of the master integrals defined in eq. (4.8) as Note that all poles explicitly cancel between the different master integrals, as infrared poles are regulated by the offshellness p 2 > 0 and UV poles are regulated by b 2 > 0.

Sail diagram
The sail topology of Fig. 3b and its mirror diagram are given bỹ For compactness, we parameterize the Wilson lines by a path γ(s), such that where γ µ (s) = dγ µ (s)/ds and γ(s) is composed of three straight segments given by For brevity, here we suppress the explicit dependence of the γ i (s) on b µ and L. After evaluating the Dirac trace in eq. (4.13), the Feynman gauge piece can be expressed as Note that the terms proportional to δ µλ and δ µρ cancel each other in the case ρ = λ, and that only the term proportional to k µ yields a 1/ pole, which can easily be extracted since k · γ (s) involves a total derivative in s. We find where we have made the UV pole in 1/ explicit. In the covariant-gauge piece in eq. (4.13), the derivatives of the path always combine to k·γ i (s), such that the ds integration only involves a total derivative, i.e. one only encounters This gives a simple result in terms of master integrals, This result contains a UV pole inducing a logarithmic contribution, given by Axial projection. The sail diagram is the only diagram contributing for the axial projector P = 1 2 γ ρ γ 5 . It is obtained similar to eq. (4.13) as The gauge-dependent piece is easily seen to vanish, such that only the Feynman piece needs to be considered. The relevant traces are given by where the antisymmetric tensor is normalized such that 0123 = 1. Inserting into eq. (4.21), we obtaiñ Here, the I 1,1 are the usual master integrals defined in eq. (4.8). Note that both I 1,1 and I ν 1,1 are IR and UV finite, and hence eq. (4.24) does not contain any poles in . In particular, this implies that there is no ambiguity in defining γ 5 in d = 4 − 2 dimensions for this calculation. The γ i (s) in eq. (4.24) are the three line segments defined in eq. (4.15).

Tadpole diagram
The Wilson line self energy, Fig. 3c, is given bỹ where as in Sec. 4.1.2 γ(s) is the Wilson line path and we included a symmetry factor 1/2. The Feynman piece can be obtained from Ref. [119], . (4.26) The covariant piece only involves an integral over a total derivative in s, and is given by where the two pieces are given bỹ The individual pieces can be found in eqs. (4.11), (4.17) and (4.26) forq ρλ (1) , and in eqs. (4.12), (4.19) and (4.27) for ∆q ρλ (1) , respectively. We have checked that the poles in agree with those reported in Ref. [119,121], and verified numerically that after dropping these poles our result at b z = 0 agrees with Ref. [121]. Note that our results are significantly more involved than those in Ref. [121] because we keep b z = 0, which is necessary for the quasi-beam function that is needed as input for the calculation of γ q ζ (µ, b T ). For the axial projection, there is only one nonvanishing contribution, such that where q ρλ (1) a s is given in eq. (4.24).

RI /MOM renormalization factor and conversion to MS
Having calculated the full one-loop result for the off-shell amputated Green's functionq ρλ , we can now proceed to calculate the RI /MOM renormalization and the conversion to the MS scheme. This also requires the one-loop wave function renormalization to account for the external state in the amputated Green's function. In the RI /MOM scheme it is given by [120] The RI /MOM renormalization of the quasi-TMDPDF also requires us to include the oneloop soft factor. Using eq. (2.12), it can be written as The required one-loop result for the bent soft function can be obtained from Ref. [119] as S q (1) bent (b T , µ, , L) = 12 + 12 ln The RI /MOM to MS conversion kernel follows from eqs. (3.6), (3.7) and (3.10) as where all the 1/ poles are canceled byZ MS q (µ, ), so only the O( 0 ) terms are extracted from the terms in the square brackets.
Last but not least, one may note that∆ q S (b R T , , L) also formally cancels out in the ratios in eqs. (2.18) and (2.25) due to its b z -independence, and therefore equivalently we can drop it in eq. (4.34) and obtain the conversion factor that matches the RI /MOM-renormalized quasi-beam function to the MS scheme, However, thisZ B will suffer from L/b T divergence that makes its numerical value much larger than one, indicating that the perturbation series does not converge. In constrast,Z q in eq. (4.34), which includes the correction from∆ q S , is free from such divergences and has good perturbative convergence, as we will demonstrate numerically in the following section.

Numerical results
In this section, we numerically illustrate the importance of the perturbative matching from the RI /MOM to the MS scheme. We assume a lattice with spacing a = 0.06 fm and size L lat = 32 a, and set the length of the Wilson line to L = 10 a. The MS renormalization scale is chosen as µ = 3 GeV, with α s (µ) = 0.2492 obtained using three-loop running from α s (m Z ) = 0.118. We always work in Landau gauge with ξ = 0. To show the effect of canceling linear divergences in L/b T , we will consider both the conversion factorZ q for the quasi-TMDPDF andZ B for the quasi-beam function alone.
We first consider the Euclidean momentum 2π L lat ≈ (3.9, 1.9, 1.9, 1.0) GeV , In Fig. 4, we showZ q in the left panel andZ B in the right panel. The b T dependence is shown for fixed b z = 0 (top row) and b z = 3a (middle row), while the bottom row shows the b z dependence for fixed b T = 5a. In each plot, we show real and imaginary parts for the ρ = λ = 0 Dirac structure in solid red and dotted green, respectively, as well as the real and imaginary parts for ρ = λ = 3 in dashed blue and dashed orange, respectively. The imaginary parts in the first two rows are amplified by a factor of ten to increase their visibility. The off-diagonal Dirac structures with ρ = λ are very small and not shown here. In all cases, we find a very small imaginary part of both Z B and Z q , and that the two choices λ = ρ = 0 and λ = ρ = 3 are very similar. Hence in the following, we restrict our discussion to the real part and the choice ρ = λ = 0 only. ForZ B (right panels of Fig. 4), the presence of the L/b R T divergence is clearly visible, and leads to large values for this factor. Since this coefficient is 1 at lowest order, clearly perturbation theory is not converging forZ B , as anticipated.
ForZ q (left panels of Fig. 4), we generically observe corrections close to Z q = 1, indicating that the O(α s ) corrections are rather moderate and of the expected size of a NLO  correction. However, there is a significant dependence on both b R T and b z . In particular, one can observe a mild logarithmic dependence on b R T as b R T → 0. Since b R T is a free parameter in the renormalization procedure, one can choose it freely to yield small matching corrections, as long as b R T is perturbative. The results in Fig. 4 indicate that b R T = O(a) is a good choice. There is also a significant b z dependence, arising from the fact that in the RI /MOM scheme -25 - one fully absorbs the b z -dependence at p = p R and b T = b R T into the UV renormalization, and this b z dependence must therefore be corrected perturbatively through the conversion to the MS scheme.
For larger b z the correction fromZ q becomes numerically significant, as can seen from the bottom left panel of Fig. 4. The impact of this large b z region is suppressed by the large parton momentum when the quasi-TMDPDF is Fourier transformed into the x-space, namely through the oscillation caused by the Fourier exponents involving (xP z b z ) in eq. (2.18). For the position space version in eq. (2.25) the analogous suppression of the large b z region occurs from the falling and oscillating behavior ofC ns with (P z b z ). The derivation of eq. (2.14), and thus both the momentum and position space formulae for γ q ζ , assume the hierarchy b z L, and hence dominance of the integral away from the large b z ∼ L region. Studying numerically the dependence of eqs. (2.18) and (2.25) on the upper limit of b z used in the integrations will help us understand how well this hierarchy is obeyed.
Finally, we study the p z R dependence for fixed b T = 5 a and b z = 3 a. For a given p z R there is still considerable freedom for the other parameters of p R , and we choose the Euclidean momentum where p 0 R is a function of p z R such that p 2 R is fixed. The largest value of p z R yielding a real solution for p 0 R is then given by p z R = √ 79(2π/L lat ) ≈ 5.7 GeV. Fig. 5 shows the resulting scheme conversion factors, withZ q on the left andZ B on the right. As before, in both cases the imaginary part (blue dashed) is very small, and the real part (red) is close to unity forZ q indicating that perturbation theory is working as expected, while it has significant deviation from unity forZ B (where we see that perturbation theory is breaking down). ForZ q it can also be observed that there is a relatively mild dependence on p z R .

Conclusion
In this paper we have elaborated on the method to determine the Collins-Soper kernel using ratios of quasi-TMDPDFs. Originally, in Ref. [118] a method was proposed which used ratios of properly matched and renormalized quasi-TMDPDFs in momentum space. This requires a Fourier transformation of spatial correlations obtained from the lattice to momentum space, which can be numerically challenging. Here we have extended this proposal to demonstrate how to carry out the matching for renormalized ratios directly in position space. This trades the Fourier transformation for a convolution with a position space matching coefficient, which we expect will improve the numerical stability of the method. The required position space matching coefficientC ns (µ, y, P z ) was obtained here at O(α s ).
In addition, we have calculated a renormalization scheme conversion factor that is needed for the lattice calculation. Renormalization on the lattice must necessarily be done nonperturbatively to properly handle power law divergences from spatial Wilson line self energies. Here we calculated the one-loop renormalization factor for the transversemomentum dependent quasi-TMDPDFs in the regularization-independent momentum subtraction RI /MOM scheme with b z = 0, and used this result to obtain the one-loop conversion factorZ q (b z , µ,μ) that converts from the RI /MOM scheme to the MS scheme. This conversion factor is necessary to obtain results for the Collins-Soper kernel γ q ζ (µ, b T ) in the desired MS scheme. Our results are thus key to determining the Collins-Soper kernel from lattice QCD using ratios of quasi-TMDPDFs as proposed in Ref. [118], and elaborated on in Ref. [119]. These results will also be used in the lattice study of nonperturbative renormalization of the quasi-beam functions [128].
Together the results obtained here provide important ingredients to enable a first nonperturbative determination of the Collins-Soper kernel from lattice QCD.
x → 1 can be extracted using the asymptotic limit of K n (z → 0), where the integral kernels are defined as B Alternative determination of γ ζ in position space Here, we present a slightly modified method compared to that presented in Sec. 2.4 for determining γ ζ in position space. There, the Fourier transform of the inverse of the kernel C ns was employed. Here, we directly Fourier transform the kernel C ns , As before, we Fourier transform with respect to b z i P z i to absorb superfluous factors of P z i . Plugging eqs. (2.20) and (B.1) into eq. (2.15), we get P z 2 db z 1 db z 2 e ix(b z 1 P z 1 +b z 2 P z 2 )C ns (µ, b z 2 P z 2 , P z 2 )f ns (b z 1 , b T , µ, P z 1 ) (B.2) Next, we Fourier transform both sides from x to y by integrating over x against e −ixy , obtaining db z 1C ns (µ, y − b z 1 P z 1 , P z 2 )f ns (b z 1 , b T , µ, P z 1 ) (B.3) This can trivially be solved for γ q ζ as γ q ζ (µ, b T ) = 1 ln(P z 1 /P z 2 ) ln db zC ns (µ, y − b z P z 1 , P z 2 )f ns (b z , b T , µ, P z 1 ) db zC ns (µ, y − b z P z 2 , P z 1 )f ns (b z , b T , µ, P z 2 ) .

(B.4)
Using the expression eq. (2.7) forf ns and inserting a factorR B to separately cancel divergences in b T /a, L/a and L/b T in numerator and denominator, we obtain the final expression γ q ζ (µ, b T ) = 1 ln(P z 1 /P z 2 )
The key difference to eq. (2.25) is that in eq. (B.5), both numerator and denominator depend on P z 1 and P z 2 , sinceC ns depends on both momenta. In contrast, in eq. (2.25) the numerator only depends on P z 1 and the denominator only depends on P z 2 , which makes the bookkeeping simpler for an analysis that separately determines the numerator and denominator before taking ratios.