Towards Quasi-Transverse Momentum Dependent PDFs Computable on the Lattice

TL;DR

This paper explores the construction of quasi-TMDPDFs for lattice QCD, analyzing their infrared properties and proposing a modified soft matrix element that better approximates true TMDPDFs.

Contribution

It introduces a new 'bent' quasi soft matrix element that passes the infrared consistency test at one loop, advancing lattice calculations of TMDPDFs.

Findings

Standard quasi-TMDPDF definitions fail the infrared test at one loop.

Ratios of impact parameter quasi-TMDPDFs can be matched to TMDPDF ratios.

The 'bent' quasi soft matrix element passes the infrared test at one loop.

Abstract

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.

Figures

Tables

Table 1. Table 1 : Dependence of unsubtracted beam function B q ≡ B q ( unsub ) subscript 𝐵 𝑞 superscript subscript 𝐵 𝑞 unsub B_{q}\equiv B_{q}^{\rm(unsub)} , soft subtraction Δ S q superscript subscript Δ 𝑆 𝑞 \Delta_{S}^{q} and TMDPDF f q TMD superscript subscript 𝑓 𝑞 TMD f_{q}^{\mathrm{TMD}} (upper part) and their quasi constructions B ~ q subscript ~ 𝐵 𝑞 \tilde{B}_{q} , Δ ~ S q superscript subscript ~ Δ 𝑆 𝑞 \tilde{\Delta}_{S}^{q} , f ~ q TMD superscript subscript ~ 𝑓 𝑞 TMD \tilde{f}_{q}^{\mathrm{TMD}} (lower part) on the logarithm L b = ln ⁡ ( b T 2 ​ μ 2 / b 0 2 ) subscript 𝐿 𝑏 superscript subscript 𝑏 𝑇 2 superscript 𝜇 2 superscript subscript 𝑏 0 2 L_{b}=\ln(b_{T}^{2}\mu^{2}/b_{0}^{2}) in various rapidity regularization schemes. Results are shown for terms from the one loop matrix elements, pulling out an overall α s / ( 2 ​ π ) subscript 𝛼 𝑠 2 𝜋 \alpha_{s}/(2\pi) . Only pure L b subscript 𝐿 𝑏 L_{b} terms are shown. The full functional form for B q subscript 𝐵 𝑞 B_{q} , Δ S q superscript subscript Δ 𝑆 𝑞 \Delta_{S}^{q} , and f q TMD superscript subscript 𝑓 𝑞 TMD f_{q}^{\mathrm{TMD}} can be found in appendix  B for all regulators. The corresponding results are given in Eqs. ( 53 ), ( 4.3 ) and ( 67 ) for the naive quasi-TMDPDF, and in Eqs. ( 53 ), ( 78 ) and ( 4.5 ) for the quasi-TMDPDF using the bent soft function.

Regulator	Beam function $B_q$	Soft factor ${\Delta }_S^q$	TMDPDF $f_q^{\mathrm{TMD}}=B_q{\Delta }_S^q$
Collins Collins:1350496	$-\frac{1}{2}{L}_b^2,\frac{5}{2}{L}_b$	$-L_b$	$-\frac{1}{2}{L}_b^2,\frac{3}{2}{L}_b$
$\delta$ regulator GarciaEchevarria:2011rb ; Echevarria:2012js	$\frac{3}{2}{L}_b$	$-\frac{1}{2}{L}_b^2$	$-\frac{1}{2}{L}_b^2,\frac{3}{2}{L}_b$
$\eta$ regulator Chiu:2011qc ; Chiu:2012ir	$\frac{3}{2}{L}_b$	$-\frac{1}{2}{L}_b^2$	$-\frac{1}{2}{L}_b^2,\frac{3}{2}{L}_b$
Exp. regulator Li:2016axz	$-L_b^2,\frac{3}{2}{L}_b$	$\frac{1}{2}{L}_b^2$	$-\frac{1}{2}{L}_b^2,\frac{3}{2}{L}_b$
	quasi ${\tilde{B}}_q$	quasi ${\tilde{\Delta }}_S^q$	quasi ${\tilde{f}}_q^{\mathrm{TMD}}={\tilde{B}}_q{\tilde{\Delta }}_S^q$
Finite $L$, naive ${\tilde{\Delta }}_S^q$	$-\frac{1}{2}{L}_b^2,\frac{9}{2}{L}_b$	$-2L_b$	$-\frac{1}{2}{L}_b^2,\frac{5}{2}{L}_b$
Finite $L$, bent ${\tilde{\Delta }}_S^q$	$-\frac{1}{2}{L}_b^2,\frac{9}{2}{L}_b$	$-3L_b$	$-\frac{1}{2}{L}_b^2,\frac{3}{2}{L}_b$

Equations

$$f_i(x,\mu )$$

$$\displaystyle W\Big{(}b^{+}\frac{{\bar{n}}}{2},0\Big{)}$$

$${\tilde{f}}_i(x,P^z,\tilde{\mu })$$

$$W_z(z,0)$$

$$\displaystyle\tilde{f}_{i}(x,P^{z},\tilde{\mu})=\int_{-1}^{1}\frac{\mathrm{d}y}{|y|}C_{ij}\biggl{(}\frac{x}{y},\frac{\tilde{\mu}}{\mu},\frac{\mu}{|y|P^{z}}\biggr{)}f_{j}(y,\mu)+\mathcal{O}\biggl{(}{\Lambda_{\mathrm{QCD}}^{2}\over P_{z}^{2}},{M^{2}\over P_{z}^{2}}\biggr{)}\,,$$

$$\frac{\mathrm{d}\sigma }{\mathrm{d}Q\mathrm{d}Y{\mathrm{d}}^2{\vec{q}}_T}$$

$$\displaystyle\hskip 28.45274pt\times\biggl{[}1+\mathcal{O}\biggl{(}\frac{q_{T}^{2}}{Q^{2}},\frac{\Lambda_{\mathrm{QCD}}^{2}}{Q^{2}}\biggr{)}\biggr{]}$$

$$=\sum _{i,j}{H}_{ij}(Q,\mu )\int {\mathrm{d}}^2{\vec{b}}_T{e}^{\mathrm{i}{\vec{b}}_T\cdot {\vec{q}}_T}{f}_i^{\mathrm{TMD}}(x_a,{\vec{b}}_T,\mu ,{\zeta }_a)f_j^{\mathrm{TMD}}(x_b,{\vec{b}}_T,\mu ,{\zeta }_b)$$

$$\displaystyle\hskip 28.45274pt\times\biggl{[}1+\mathcal{O}\biggl{(}\frac{q_{T}^{2}}{Q^{2}},\frac{\Lambda_{\mathrm{QCD}}^{2}}{Q^{2}}\biggr{)}\biggl{]}\,.$$

$${\zeta }_a=(x_a{P}_a^{-}{)}^2{e}^{-2y_n},{\zeta }_b=(x_b{P}_b^{+}{)}^2{e}^{2y_n},{\zeta }_a{\zeta }_b=Q^4.$$

$$B_i(x,{\vec{b}}_T,\mu ,\zeta /{\nu }^2)$$

$$=\underset{\begin{array}{c}ϵ\to 0\\ \tau \to 0\end{array}}{\lim }{Z}_B^i(b_T,\mu ,\nu ,ϵ,\tau ,x{P}^{-})\frac{B_i^{(\mathrm{unsub})}(x,{\vec{b}}_T,ϵ,\tau ,x{P}^{-})}{S_i^0(b_T,ϵ,\tau )},$$

$$S^i(b_T,\mu ,\nu )$$

$$f_i^{\mathrm{TMD}}(x,{\vec{b}}_T,\mu ,\zeta )$$

$$\displaystyle\equiv\lim_{\begin{subarray}{c}\epsilon\to 0\\ \tau\to 0\end{subarray}}Z_{\rm uv}^{i}(\mu,\zeta,\epsilon)\,B_{i}^{\rm(unsub)}\bigl{(}x,\vec{b}_{T},\epsilon,\tau,xP^{-}\bigr{)}\Delta_{S}^{i}(b_{T},\epsilon,\tau)\,.$$

$$Z_{\mathrm{uv}}(\mu ,\zeta ,ϵ)=Z_B(b_T,\mu ,\nu ,ϵ,\tau ,x{P}^{-})\sqrt{Z_S(b_T,\mu ,\nu ,ϵ,\tau )}.$$

$$B_q(x,{\vec{b}}_T,ϵ,\tau ,x{P}^{-})$$

$$\displaystyle\hskip 156.49014pt\times W^{\dagger}(0)q(0)\Bigr{]}_{\tau}\Bigl{|}h(P)\Bigr{>},$$

$$S^q(b_T,ϵ,\tau )$$

$$\displaystyle\hskip 56.9055pt\times S^{\dagger}_{\bar{n}}(\vec{0}_{T})S_{n}(\vec{0}_{T})S_{T}^{\dagger}\bigl{(}-\infty n;\vec{b}_{T},\vec{0}_{T}\bigr{)}\bigr{]}_{\tau}\bigl{|}0\bigr{>}\,.$$

$$W(x^{\mu })$$

$$S_n(x^{\mu })$$

$$W_T(x^{\mu };{\vec{b}}_T,{\vec{0}}_T)$$

$$\displaystyle f_{i}^{\mathrm{TMD}}(x,\vec{b}_{T},\mu,\zeta)=\sum_{j}\int_{x}^{1}\frac{\mathrm{d}y}{y}C_{ij}\Bigl{(}\frac{x}{y},\vec{b}_{T},\mu,\zeta\Bigr{)}f_{j}(y,\mu)+\mathcal{O}(b_{T}\Lambda_{\mathrm{QCD}})\,.$$

$$B_n:$$

$$B_{\bar{n}}:$$

$$S:$$

$$p_n^2\sim p_{\bar{n}}^2\sim p_s^2\sim q_T^2,$$

$$I_{\mathrm{div}}=\int \mathrm{d}{k}^{+}\mathrm{d}{k}^{-}\frac{f(k^{+}{k}^{-})}{(k^{+}{k}^{-}{)}^{1+ϵ}}=\frac{1}{2}\int \frac{\mathrm{d}(k^{-}/k^{+})}{k^{-}/k^{+}}\int \mathrm{d}(k^{+}{k}^{-})\frac{f(k^{+}{k}^{-})}{(k^{+}{k}^{-}{)}^{1+ϵ}}.$$

$$\displaystyle I_{\rm div}=\int\mathrm{d}k^{+}\mathrm{d}k^{-}\frac{f(k^{+}k^{-})}{(k^{+}k^{-})^{1+\epsilon}}\to\int\mathrm{d}k^{+}\mathrm{d}k^{-}\frac{f(k^{+}k^{-})}{(k^{+}k^{-})^{1+\epsilon}}\biggl{|}\frac{k^{-}-k^{+}}{\nu}\biggr{|}^{-\eta}\,.$$

$$\mu \frac{\mathrm{d}}{\mathrm{d}\mu }{f}_i^{\mathrm{TMD}}(x,{\vec{b}}_T,\mu ,\zeta )$$

$$\zeta \frac{\mathrm{d}}{\mathrm{d}\zeta }{f}_i^{\mathrm{TMD}}(x,{\vec{b}}_T,\mu ,\zeta )$$

$$\mu \frac{\mathrm{d}}{\mathrm{d}\mu }{\gamma }_{\zeta }^i(\mu ,b_T)$$

$${\gamma }_{\mu }^i(\mu ,\zeta )$$

$${\gamma }_{\zeta }^i(\mu ,b_T)$$

$${\gamma }_{\mu }^i(\mu ,\zeta )$$

$${\gamma }_{\zeta }^i(\mu ,b_T)$$

Full text

††preprint:

MIT–CTP 5039

1 Introduction

Transverse momentum dependent (TMD) distributions are important ingredients for describing high-energy scatterings at low transverse momentum, such as for Drell-Yan production at the Tevatron Affolder:1999jh ; Abbott:1999yd ; Abazov:2007ac ; Abazov:2010kn or LHC Aad:2011gj ; Chatrchyan:2011wt ; Aad:2014xaa ; Khachatryan:2015oaa ; Aad:2015auj ; Khachatryan:2016nbe which is a benchmark observable of the Standard Model (SM). For Higgs bosons produced in proton collisions, the transverse momentum spectrum is also one of the primary observables describing the production kinematics and is of great interest in LHC measurements Aad:2014lwa ; Aad:2014tca ; Aad:2016lvc ; Aaboud:2017oem ; Aaboud:2018xdt ; Khachatryan:2015rxa ; Khachatryan:2015yvw ; Khachatryan:2016vnn ; Sirunyan:2018kta . TMDs also play an important role in measuring semi-inclusive deep-inelastic scattering (SIDIS) at low energies Ashman:1991cj ; Derrick:1995xg ; Adloff:1996dy ; Aaron:2008ad ; Airapetian:2012ki ; Adolph:2013stb ; Aghasyan:2017ctw , and in improving our understanding of hadron structure Boer:2011fh ; Accardi:2012qut .

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}\gg\Lambda_{\mathrm{QCD}}$, where they can be obtained in terms of collinear parton distribution functions (PDFs). In contrast, nonperturbative TMDs with small $q_{T}\sim\Lambda_{\mathrm{QCD}}$ have only been extracted from measurement by performing global fits to a variety of experimental data sets, see e.g. Refs. Landry:1999an ; Landry:2002ix ; Konychev:2005iy ; DAlesio:2014mrz ; Bacchetta:2017gcc ; Scimemi:2017etj . 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. Musch:2010ka ; Musch:2011er ; Engelhardt:2015xja ; Yoon:2016dyh ; Yoon:2017qzo , 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. Ji:2013dva ; Ji:2014gla 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 $\overline{\mathrm{MS}}$ scheme as

[TABLE]

where $i$ is the flavor index and the square bracket with subscript $\mu$ denotes that the operator is renormalized at scale $\mu$. For later use we introduce the lightlike reference vectors $n^{\mu}=(1,0,0,1)$ and ${\bar{n}}^{\mu}=(1,0,0,-1)$. The light-cone coordinates are defined as $b^{\pm}=b^{0}\mp b^{3}$, and $p(P)$ denotes a proton state with momentum $P^{\mu}=(E,0,0,P^{z})$. The path-ordered light-cone collinear Wilson line is

[TABLE]

To calculate the PDF $f_{i}(x,\mu)$ using LaMET, one starts from a quasi-PDF, which is defined from equal-time correlation functions,

[TABLE]

where the operator is renormalized with a lattice friendly renormalization scheme and $\tilde{\mu}$ denotes a corresponding scale. The spacelike Wilson line is

[TABLE]

and $\Gamma=\gamma^{0}$ or $\Gamma=\gamma^{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 Hatta:2013gta . Unlike the PDF $f_{i}(x,\mu)$ that is boost invariant, the quasi-PDF $\tilde{f}_{i}(x,P^{z},\tilde{\mu})$ depends nontrivially on the nucleon momentum $P^{z}$. When $P^{z}$ is large compared to $\Lambda_{\mathrm{QCD}}$ as well as the nucleon mass $M$, the quasi-PDF satisfies the following factorization theorem Ji:2013dva ; Ji:2014gla ; Ma:2014jla ; Ma:2017pxb ; Izubuchi:2018srq ,

[TABLE]

where $\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}/P_{z}^{2},M^{2}/P_{z}^{2})$ terms are power corrections. Here $f_{j}(y,\mu)$ for $-1<y<0$ corresponds to the anti-quark PDF. The $C_{ij}$ are perturbative matching coefficients which come from a hard region of momentum space, see Refs. Ma:2014jla ; Izubuchi:2018srq for further details.

Recently, significant progress has been made on various aspects of the LaMET procedure, including the renormalization and matching Xiong:2013bka ; Ma:2014jla ; Ma:2014jga ; Ji:2015jwa ; Ji:2015qla ; Xiong:2015nua ; Li:2016amo ; Ishikawa:2016znu ; Chen:2016fxx ; Carlson:2017gpk ; Briceno:2017cpo ; Xiong:2017jtn ; Constantinou:2017sej ; Rossi:2017muf ; Ji:2017rah ; Ji:2017oey ; Ishikawa:2017faj ; Green:2017xeu ; Wang:2017qyg ; Chen:2017mie ; Stewart:2017tvs ; Wang:2017eel ; Spanoudes:2018zya ; Izubuchi:2018srq ; Xu:2018mpf ; Rossi:2018zkn ; Zhang:2018diq ; Li:2018tpe ; Liu:2018tox of the quasi-PDF, the power corrections Chen:2016utp ; Radyushkin:2017ffo ; Braun:2018brg , as well as the lattice calculation of the $x$-dependence of PDFs and distribution amplitudes Lin:2014zya ; Alexandrou:2015rja ; Chen:2016utp ; Alexandrou:2016jqi ; Zhang:2017bzy ; Alexandrou:2017huk ; Chen:2017mzz ; Green:2017xeu ; Chen:2017lnm ; Chen:2017gck ; Alexandrou:2018pbm ; Chen:2018xof ; Chen:2018fwa ; Alexandrou:2018eet ; Liu:2018uuj ; Lin:2018qky ; Fan:2018dxu ; Liu:2018hxv . Notably, the most recent lattice results at physical pion mass Alexandrou:2018pbm ; Chen:2018xof ; Alexandrou:2018eet ; Lin:2018qky ; Liu:2018hxv and large nucleon momenta Chen:2018xof ; Lin:2018qky ; Liu:2018hxv have shown encouraging signs that the LaMET approach can lead to a precise determination of the PDFs.

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. Engelhardt:2015xja ; Yoon:2016dyh ; Yoon:2017qzo . 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 collinear111Note that the second use of the word “collinear” here is in the context of factorized collinear and soft fields as defined for example in SCET, not to distinguish between collinear and TMD factorization. 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. Ji:2014hxa , 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. Ji:2018hvs , 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. Ji:2018hvs , 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}\sim\Lambda_{\mathrm{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. (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 $\alpha_{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 $\overline{\rm 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.

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

2.1 TMD Factorization and TMDPDFs

TMD factorization was originally derived by Collins, Soper and Sterman (CSS) in Refs. Collins:1981uk ; Collins:1981va ; Collins:1984kg . Refs. Collins:1985ue ; Collins:1988ig ; Collins:1989gx ; Collins:1350496 ; Diehl:2015bca showed the cancellation of potentially factorization-violating Glauber modes, and the factorization was further elaborated on and extended in Refs. Catani:2000vq ; deFlorian:2001zd ; Catani:2010pd ; Collins:1350496 . It has also been considered in the framework of Soft-Collinear Effective Theory (SCET) Bauer:2000ew ; Bauer:2000yr ; Bauer:2001ct ; Bauer:2001yt by various authors Becher:2010tm ; Becher:2011xn ; Becher:2012yn ; GarciaEchevarria:2011rb ; Echevarria:2012js ; Echevarria:2014rua ; Chiu:2012ir . For a relation of the different approaches to each other, see e.g. Ref. Collins:2017oxh , and for a historical review on TMDPDFs we refer the reader to Collins:1350496 . In appendix A we give a summary of our notation for light-cone coordinates. In appendix B we give a comparison between different TMDPDF constructions from the literature, and provide explicit relations to the notation we use here.

For simplicity we consider TMD factorization in the context of the production of a color-singlet final state $F$ in the scattering of two unpolarized energetic protons moving along the $n^{\mu}=(1,0,0,1)$ and ${\bar{n}}^{\mu}=(1,0,0,-1)$ directions. Examples of such cross sections include Drell-Yan, $W$, and Higgs production. Here we only measure the total four momentum of $F$ through its invariant mass $Q$, rapidity $Y$ and transverse momentum $\vec{q}_{T}$, so that we only have unpolarized TMDPDFs. In the limit $q_{T}\equiv|\vec{q}_{T}|\ll Q$ the cross section can be factorized as

[TABLE]

The factorized cross section receives power corrections in $q_{T}^{2}/Q^{2}$ and $\Lambda_{\mathrm{QCD}}^{2}/Q^{2}$, and first studies of their perturbative and nonperturbative structure have been performed in Refs. Balitsky:2017flc ; Balitsky:2017gis ; Ebert:2018gsn . Importantly, Eq. (6) remains valid for nonperturbative $q_{T}\sim\Lambda_{\mathrm{QCD}}$. The factorization is most elegantly written in impact parameter space, with $\vec{b}_{T}$ being Fourier conjugate to the measured transverse momentum $\vec{q}_{T}$. In Eq. (6) $i,j$ are parton indices for quark flavors or gluons,222For gluon-induced processes, $H$ and $B$ also carry helicity indices even for unpolarized TMDPDFs which we suppress, and $H_{ij}$ is the hard function containing virtual corrections to the underlying hard process $ab\to 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_{\rm cm}$ and $x_{b}=Qe^{-Y}/E_{\rm cm}$ where $E_{\rm cm}$ is the center of mass energy of the $pp$ collision. Finally, $\zeta_{a}$ and $\zeta_{b}$ are related to the momenta of the struck partons and given by

[TABLE]

Here $P_{a}^{\mu}=P_{a}^{-}\frac{n^{\mu}}{2}+\frac{m_{P}^{2}}{P_{a}^{-}}\frac{{\bar{n}}^{\mu}}{2}$ and $P_{b}^{\mu}=P_{b}^{+}\frac{{\bar{n}}^{\mu}}{2}+\frac{m_{P}^{2}}{P_{b}^{+}}\frac{n^{\mu}}{2}$ are the momenta of the protons, such that $P_{a}^{-}P_{b}^{+}=E_{\rm cm}^{2}$, and $y_{n}$ is a parameter that encodes scheme dependence. Note that whenever possible we neglect target mass corrections from $m_{P}^{2}\ll Q^{2}<E_{\rm cm}^{2}$.

The small transverse momentum $q_{T}\ll Q$ of the final state is generated from soft and collinear radiation off the incoming protons. In Eq. (6), 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 $n$- and ${\bar{n}}$-collinear protons, and to distinguish them from the final TMDPDF we will follow the language of Ref. Stewart:2009yx and refer to $B_{i}$ and $B_{j}$ as beam functions. The soft function $S^{i}\equiv S_{ij}$ is a vacuum matrix element that encodes soft exchange between the incoming partons. It only differs for gluons, $S^{g}\equiv S_{gg}$, and quarks, $S^{q}\equiv S_{q\bar{q}}$, 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$ and ${\bar{n}}$.

The different ingredients in Eq. (6) depend on both the renormalization scale $\mu$ and an additional rapidity renormalization scale $\nu$.333Some regulators and schemes have two distinct scales $\nu_{n},\nu_{\bar{n}}$ where $\nu_{n}$ appears in the $n$-collinear beam function $B_{i}$, $\nu_{\bar{n}}$ appears in the ${\bar{n}}$-collinear beam function $B_{j}$, and both scales appear in $S_{ij}$. We suppress this possibility in our review. The latter scale arises because the matrix elements not only suffer from UV divergences, but also from so-called rapidity divergences that require a dedicated regulator Soper:1979fq ; Collins:1981uk ; Collins:1992tv ; Collins:2008ht ; Becher:2010tm ; GarciaEchevarria:2011rb ; Chiu:2011qc ; Chiu:2012ir , and we will discuss their physical origin in more detail in Sec. 2.2. Crucially, the regularization of these divergences requires an explicit rapidity regulator, which we generically denote as $\tau$. One can distinguish two classes of schemes by whether introducing the rapidity regulator affects the hard function $H_{ij}$ or not. The hard function describes virtual corrections to the underlying Born process and thus carries the full process dependence (the TMDPDFs are only sensitive to the hard initial state, namely $q\bar{q}$ or $gg$), so it is simplest to consider regulators that do not affect $H_{ij}$. In the following, we will only consider such regulators, and fix the hard function to be in the $\overline{\mathrm{MS}}$ scheme, which then immediately fixes the scheme for the product of the TMDPDFs. This applies to most modern regulators, but not to the formulations of Refs. Soper:1979fq ; Collins:1981uk ; Ji:2004wu where $H_{ij}$ and $f_{i}^{\mathrm{TMD}}$ depend on an additional parameter $\rho$.

For schemes considered here we can generically define UV and rapidity-renormalized beam and soft functions as

[TABLE]

Note that for simplicity we use the same notation for bare and renormalized soft functions $S^{i}$, since it will always be clear from the arguments or context which one we are discussing. For the soft function, this is straightforward. For the beam function, we note that the renormalized beam function depends on $\zeta\propto(xP^{-})^{2}$, where the proportionality is scheme-dependent, i.e. depends on the precise definition of $\tau$. Secondly, we note that $B_{i}$ is defined to only describe collinear radiation. In practical calculations, one often encounters an overlap with the soft function when the collinear radiation becomes soft. To avoid double counting with the soft function, one has to subtract out this overlap from the unsubtracted beam function $B_{i}^{\rm(unsub)}$, and the subtraction factor is denoted as $S_{i}^{0}$ in Eq. (2.1). (In SCET, this is referred to as soft zero-bin subtraction Manohar:2006nz .) Since $S_{i}^{0}$ describes soft physics, it cannot depend on the large momentum $xP^{-}$. As indicated by the notation for $S_{i}^{0}(b_{T},\epsilon,\tau)$, this factor is often equivalent or closely related to the soft function $S^{i}(b_{T},\epsilon,\tau)$ itself, but this is not always the case. In particular, one can find rapidity renormalization schemes which (a) have no zero-bin subtraction, $S_{i}^{0}=1$ Chiu:2012ir , (b) where the zero-bin is equal to the soft function, $S_{i}^{0}=S^{i}$ GarciaEchevarria:2011rb ; Li:2016axz , and (c) where the subtraction is given by a combination of soft functions in different non-lightlike directions Collins:1350496 .

In practice, one often combines the beam and soft functions to yield two separate TMDPDFs $f^{\mathrm{TMD}}_{i}$ and $f^{\mathrm{TMD}}_{j}$, as used in the factorized cross section in Eq. (6). This can be achieved either by combining the renormalized beam and soft functions,

[TABLE]

or by combining the bare functions and performing the UV renormalization afterwards,

[TABLE]

In Eq. (10), the rapidity renormalization scale $\nu$ cancels between $B_{i}$ and $S^{i}$, leaving only a dependence on $\zeta$ in the TMDPDF. Likewise, in Eq. (2.1) the dependence on the regulator $\tau$ cancels between $B_{i}$ and $S^{i}$. Note that Eq. (2.1) has been written in terms of the unsubtracted beam function $B_{i}^{\rm(unsub)}$, as is common in the literature, and for ease of notation we have combined the zero-bin subtraction $S_{i}^{0}$ and the soft function $S^{i}$ to define the soft factor $\Delta_{S}^{i}=\sqrt{S^{i}}/S_{i}^{0}$. The UV renormalization factor $Z_{\rm uv}$ in Eq. (2.1) is trivially related to those of $B_{i}$ and $S^{i}$ in Eqs. (2.1) and (9),

[TABLE]

As before, the $\zeta$ scale arises as a rapidity-regulator dependent function with the general form $\zeta\propto(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 Soper:1979fq ; Collins:1981uk . Since the same non-lightlike axial gauge has to be used in the calculation of the hard function, this regulator does not fall in the class of rapidity regulators considered here.444Note that in their original work, the hard function was absorbed in the TMDPDFs. The separation into process-dependent form factors and process-independent TMDPDFs was first noted in Ref. Catani:2000vq . In Ref. Ji:2004wu rapidity divergences are regulated by taking Wilson lines off the light cone. Since the same Wilson lines enter the hard factor, this again falls outside the class of regulators we consider, as is explicitly visible by a dependence of their hard function and TMDPDFs on an extra parameter $\rho$. This is distinct from Collin’s new regulator Collins:1350496 where again Wilson lines are taken off the lightcone, but the lightlike limit is taken such that the hard function is independent of the regulator. The regulators used in the SCET-based approaches are the analytic regulator acting on eikonal propagators Beneke:2003pa ; Chiu:2007yn ; Becher:2011dz , the $\eta$-regulator inserted into Wilson lines on the light-cone Chiu:2011qc ; Chiu:2012ir , the $\delta$-regulator which adds mass-like terms to eikonal propagators Chiu:2009yx ; GarciaEchevarria:2011rb , and the exponential regulator inserted into the phase space Li:2016axz . The definition of the TMDPDF in terms of bare beam and soft functions, Eq. (2.1), is used in Refs. Collins:1350496 ; Ji:2004wu ; GarciaEchevarria:2011rb ; Echevarria:2012js , where renormalized beam and soft functions were not defined. Both the bare and renormalized forms can be used in Refs. Chiu:2012ir ; Li:2016axz . In the analytic regulator approach of Refs. Becher:2010tm ; Becher:2012yn , 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.

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

[TABLE]

where $b^{\mu}=b^{+}{\bar{n}}^{\mu}/2+b_{T}^{\mu}$. Here the bracket $[\cdots]_{\tau}$ denotes that the operator inside is considered by implementing the rapidity regulator $\tau$. For the matrix element in Eq. (2.1), we note that diagrams that have no fields contracted with the states $|h(P)\rangle$ are excluded. 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,

[TABLE]

as well as finite-length gauge links with transverse paths,

[TABLE]

The Wilson paths for matrix elements in Eqs. (2.1) and (2.1) are shown in the $(b^{+},b^{-},\vec{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 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. Ji:2002aa ; Belitsky:2002sm ; Idilbi:2010im ; GarciaEchevarria:2011md . These gauge links are also important for constructing the analogues of Eqs. (2.1) and (2.1) on a finite-size lattice.

Note that the inclusion of a rapidity regulator can in principle spoil the gauge invariance of the matrix elements in Eqs. (2.1) and (2.1). Gauge invariance trivially holds for regulators only affecting the Wilson line paths, as for example Collins’ regulator Collins:1350496 , the exponential regulator Li:2016ctv , or the finite-length Wilson lines to be introduced in Sec. 3 for lattice calculations. Gauge invariance has been explicitly shown to hold in the $\tau\to 0$ limit for the $\eta$ regulator Chiu:2012ir , the analytic phase space regulator of Ref. Becher:2011dz , while it is known to be violated in for TMDs with the analytic regulator used in Ref. Becher:2010tm . Gauge invariance is also known to be more tricky for the $\Delta$ regulator, where the limit $\Delta\to 0$ is also required Echevarria:2015usa ; Echevarria:2015byo , and individual beam function matrix elements may only be gauge invariant after including 0-bin subtractions Chiu:2009yx .

Lastly, note that extracting TMDPDFs from lattice QCD (or experiment) is only necessary for nonperturbative $q_{T}\sim b_{T}^{-1}\sim\Lambda_{\mathrm{QCD}}$. For perturbative values, one can instead perform an operator-product expansions to match the TMDPDF, or equivalently the beam function, onto the collinear PDF Collins:1981uw ; Collins:1984kg ,

[TABLE]

Here, $C_{ij}$ are perturbative matching kernels that are known to NNLO Catani:2011kr ; Catani:2012qa ; Gehrmann:2014yya ; Luebbert:2016itl ; Echevarria:2015byo ; Echevarria:2016scs , and even to N3LO for the soft contribution Li:2016ctv . The nonperturbative input is given solely by the standard longitudinal PDFs. Throughout this work, we will limit our discussion to the case $b_{T}\sim\Lambda_{\mathrm{QCD}}^{-1}$, where Eq. (17) cannot be applied.

2.2 Rapidity Divergences in TMDs

Quantum corrections to the beam and soft function defined in Eqs. (2.1) and (2.1) 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 Collins:1981uk ; Collins:1992tv ; Manohar:2006nz ; Collins:2008ht ; Becher:2010tm ; GarciaEchevarria:2011rb ; Chiu:2011qc ; Chiu:2012ir . Rapidity divergences arise because $n$- and ${\bar{n}}$-collinear beam functions and the soft function, or equivalently the collinear and soft Wilson lines, are defined to describe modes with momenta scaling as

[TABLE]

where we follow the usual SCET conventions, and use the light-cone notation $p=(p^{+},p^{-},p_{T})$, see appendix A, and $\lambda\sim q_{T}/Q\ll 1$. These modes indicate approximations that are used to derive the corresponding operators and their Feynman rules. This is illustrated in Fig. 2, where the orange dots denote the dominant region for the $n$- and ${\bar{n}}$-collinear modes, and the green dots denotes that for the soft modes. Hard modes mediating the underlying hard process are shown in blue. Although these dots indicate the dominant momentum region, in matrix elements the corresponding fields are still integrated over all momenta. Since collinear and soft momenta have the same virtuality

[TABLE]

which in Fig. 2 is shown by lying on the same hyperbola, the soft and collinear momentum regions are only distinguished by their rapidity $y=1/2\ln(k^{-}/k^{+})$. In calculations using regulators such as dimensional regularization, which only regulate virtualities, one can therefore encounter additional rapidity divergences which arise in soft and collinear matrix elements when integrating over all $y$, and have to be resolved using a dedicated regulator. This is indicated schematically 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 Vladimirov:2016dll ; Vladimirov:2017ksc ,

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,

[TABLE]

Here the integrand only depends on the product $k^{+}k^{-}$. Singularities as $k^{+}k^{-}\to 0$ or $k^{+}k^{-}\to\infty$ 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. (20) require an additional regulator. A large variety of regulators has been suggested in the literature, see Refs. Collins:1981uk ; Collins:1350496 ; Ji:2004wu ; Beneke:2003pa ; Chiu:2007yn ; Becher:2011dz ; Chiu:2011qc ; Chiu:2012ir ; Chiu:2009yx ; GarciaEchevarria:2011rb ; Li:2016axz ; Ebert:2018gsn . The key idea in all rapidity regulators is to regulate the behavior as $y=1/2\ln(k^{-}/k^{+})\to\pm\infty$, which in the example of Eq. (20) amounts to lifting the dependence on $k^{+}k^{-}$ only. For example, in the $\eta$ regulator of Refs. Chiu:2011qc ; Chiu:2012ir one regulates the integral Eq. (20) by inserting a factor $|2k^{z}/\nu|^{-\eta}$, such that

[TABLE]

Divergences as $k^{-}/k^{+}\to 0$ or $k^{-}/k^{+}\to\infty$ are thus made manifest as poles in $\eta$. Similar to dimensional regularization, one introduces a new dimensionful scale $\nu$ to keep the regulator scaleless, and a parameter $\eta$ to be taken to zero at the end of the calculation, analogous to the limit $\epsilon\to 0$. In this way, renormalized beam and soft functions acquire an additional scale dependence, namely $\nu$, or equivalently the bare functions entering Eq. (2.1) depend on $\eta$. The chosen rapidity regulator determines the precise form of this dependence on $\nu$ and $\eta$, respectively. As discussed, this dependence cancels in the combination $f_{i}^{\mathrm{TMD}}=B_{i}\sqrt{S^{i}}$, but the cancellation leaves a residual dependence on the scale $\zeta$. Intuitively, the definitions of $\zeta_{a}$ and $\zeta_{b}$ in Eq. (7) reflects how the hyperbola in Fig. 2 is split between soft and collinear modes, see e.g. Ref. Echevarria:2012js for more details.

2.3 Evolution of TMDPDFs

The cross section in Eq. (6) must be independent of the unphysical scale $\mu$ and of the precise choice of $\zeta_{a}$ and $\zeta_{b}$ as long as $\zeta_{a}\zeta_{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. Chiu:2012ir ; Ebert:2016gcn and in appendix B.3. The RGEs for $f_{i}^{\mathrm{TMD}}$ read

[TABLE]

where $\Gamma_{\rm cusp}^{i}[\alpha_{s}]$ is the cusp anomalous dimension. The second equation in Eq. (2.3) is known as the Collins-Soper equation Collins:1981va ; Collins:1981uk . The subscripts $\mu$ and $\zeta$ on the anomalous dimensions $\gamma_{\mu}^{i}$ and $\gamma_{\zeta}^{i}$ denote the scale evolution they govern, and the superscript $i=q,g$ distinguishes quarks from gluons. These anomalous dimensions are defined as

[TABLE]

Here, we used that the $\nu$-dependence in $f^{\mathrm{TMD}}=B\sqrt{S}$ cancels and that $B$ only depends on the combination $\nu^{2}/\zeta$, and thus $\gamma_{\zeta}^{i}$ can be equivalently obtained from either $f^{\mathrm{TMD}}$ or $B$ or $S$. Since the definition of $S$, Eq. (2.1), is independent of the hadron state entering the TMDPDF and the light quark flavor, this immediately implies that $\gamma_{\zeta}^{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

[TABLE]

where $\gamma^{i}[\alpha_{s}]$ denotes the noncusp piece.

Note that $\gamma_{\zeta}^{i}(\mu,b_{T})$ has an intrinsically nonperturbative component for $b_{T}\sim\Lambda_{\mathrm{QCD}}^{-1}$, independently of the scale $\mu$, as is clear from Eq. (2.3). However, once it has been nonperturbatively determined at a scale $\mu_{0}\gg\Lambda_{\mathrm{QCD}}$, it can be perturbatively evolved to any other scale $\mu\gg\Lambda_{\mathrm{QCD}}$ using

[TABLE]

The boundary term $\gamma_{\zeta}^{i}(\mu_{0},b_{T})$ is known to three loops for perturbative $b_{T}$ and $\mu_{0}$ Li:2016axz ; Li:2016ctv ; Vladimirov:2016dll .

Combining the solutions to Eq. (2.3), the TMDPDF can be evolved from arbitrary initial scales $(\mu_{0},\zeta_{0})$ to the desired final scales $(\mu,\zeta)$ through

[TABLE]

Here, we have chosen to first evolve $\zeta_{0}\to\zeta$ at fixed $\mu_{0}$, and then $\mu_{0}\to\mu$. Since the evolution must be path independent, other choices are also possible, see also Refs. Chiu:2012ir ; Scimemi:2018xaf for a discussion of this path independence. Also note that due to the $b_{T}$ dependence of $\gamma_{\zeta}$, the TMD evolution is severely more complicated in momentum ($\vec{q}_{T}$) space Ebert:2016gcn , which is part of the reason why the factorization is commonly written in $\vec{b}_{T}$ space.

Eq. (26) is crucial to relate the TMDPDF at reference scales $(\mu_{0},\zeta_{0})$, where they are either measured or determined from lattice QCD, to the phenomenological scales $(\mu,\zeta)$.

The boundary term $f_{i}^{\mathrm{TMD}}(x,\vec{b}_{T},\mu_{0},\zeta_{0})$ in Eq. (26) is by definition nonperturbative. For perturbative $b_{T}\ll\Lambda_{\mathrm{QCD}}^{-1}$, one can can match it onto the collinear PDF, see Eq. (17), thereby reducing all the nonperturbative input to these more well-studied PDFs. In this case, the RG-evolved TMDPDF Eq. (26) serves to resum large logarithms $\ln(\mu b_{T})$ and $\ln(\zeta b_{T}^{2})$, which otherwise spoil the perturbative convergence of the matching kernel $C_{ij}$ in Eq. (17).

In this paper, we are instead interested in obtaining nonperturbative TMDPDFs valid also for the case $b_{T}\sim\Lambda_{\rm QCD}^{-1}$, where the boundary term $f_{i}^{\mathrm{TMD}}(x,\vec{b}_{T},\mu_{0},\zeta_{0})$ in Eq. (26) is obtained from lattice QCD (or equivalently is measured from experiment). In this case, the values $\mu_{0}$ and $\zeta_{0}$ are fixed by the lattice calculation (or measurement), and the role of Eq. (26) is to evolve $f_{i}^{\mathrm{TMD}}(x,\vec{b}_{T},\mu_{0},\zeta_{0})$ to the scales $(\mu,\zeta)$ required in the phenomenological application. This is completely analogous to the determination of collinear PDFs, which are extracted at a reference scale $\mu_{0}$ and then evolved via the DGLAP evolution. Similar to the DGLAP evolution, the $\mu$ evolution of the TMDPDF encoded in the first exponential in Eq. (26) is perturbative as long as $\mu_{0}$ and $\mu$ are perturbative. In contrast, the $\zeta$ evolution intrinsically contains a nonperturbative component $\gamma_{\zeta}^{i}(\mu,b_{T})$ for $b_{T}\sim\Lambda_{\mathrm{QCD}}^{-1}$, even if the TMDPDF is extracted at perturbative $(\mu_{0},\zeta_{0})$. Thus one needs to determine both $f_{i}^{\mathrm{TMD}}(x,\vec{b}_{T},\mu_{0},\zeta_{0})$ and $\gamma_{\zeta}(\mu,b_{T})$ nonperturbatively. In particular, the nonperturbative determination of $\gamma_{\zeta}(\mu,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}$ Scimemi:2016ffw . A dedicated discussion of the extraction of $\gamma_{\zeta}(\mu,b_{T})$ from lattice QCD, exploiting some of the results discussed in this paper, has been given in Ref. Ebert:2018gzl .

3 Towards Constructing quasi-TMDPDFs

The goal of this work is to define a quasi-TMDPDF $\tilde{f}^{\mathrm{TMD}}$ which involves matrix elements that are calculable with lattice QCD, and which can be matched onto the TMDPDF $f^{\mathrm{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}\sim\Lambda_{\mathrm{QCD}}^{-1}$. 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” $\tilde{B}_{i}$ and “quasi soft function” $\tilde{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.1). 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 power law divergences from Wilson line self energies that have to be subtracted. Since the self energies are also proportional to the length of the Wilson lines, and $L$ dependence cancels when combining $\tilde{B}_{i}$ and $\tilde{S}^{i}$, this gives a $b^{z}$-dependent 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 $\overline{\mathrm{MS}}$ scheme analogous to Eq. (2.1) as

[TABLE]

Here, $\tilde{B}_{i}\equiv\tilde{B}_{i}^{\rm(unsub)}$ and $\tilde{\Delta}_{S}^{i}$ 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.1). We will always consider the unsubtracted beam function and absorb the zero-bin subtraction factor into $\tilde{\Delta}_{S}^{i}$, so for simplicity we drop the superscript “(unsub)”. $\tilde{Z}_{\rm uv}$ is the lattice renormalization constant, and $\tilde{Z}^{\prime}$ converts from the lattice renormalization scheme to the $\overline{\mathrm{MS}}$ scheme. These schemes are typically distinct, and hence we distinguish the lattice renormalization scale $\tilde{\mu}$ from the $\overline{\mathrm{MS}}$ scale $\mu$. The finite lattice spacing $a$ takes the role of $\epsilon$ as an UV regulator. On the lattice, one also has to truncate the Wilson lines that enter $\tilde{B}_{i}$ and $\tilde{\Delta}_{S}^{i}$ 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 $\tau$ dependence is replaced by additional dependence on $L$. The $L$ dependence associated with both the rapidity regularization and finite length of Wilson lines cancels between $\tilde{B}_{i}$ and $\tilde{\Delta}_{S}^{i}$, and hence $\tilde{f}_{i}^{\mathrm{TMD}}$ does not depend on $L$. (In practice, there will be a finite $L$ dependence from power corrections that vanishes in the limit $L\to\infty$, and we suppress these in our notation.) Finally, $\tilde{f}_{i}^{\mathrm{TMD}}$ 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 $\zeta$ scale.

3.1 Constraints on the matching relation

The main focus of this paper is to give constructions of $\tilde{B}_{i}$ and $\tilde{\Delta}_{S}^{i}$ such that $\tilde{f}_{i}^{\mathrm{TMD}}$ can be perturbatively matched onto the TMDPDF $f_{i}^{\mathrm{TMD}}$. In this section we 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

[TABLE]

where $i,j=q,g$ are the parton flavors and $f^{\mathrm{TMD}}_{j}$ for negative $y$ corresponds to the antiquark distribution, $f^{\mathrm{TMD}}_{j}(-y,\ldots)=f^{\mathrm{TMD}}_{\bar{j}}(y,\ldots)$. The precise dependence of the matching kernel $C_{ij}^{\mathrm{TMD}}$ on its arguments is not obvious a priori, see for example Eq. (5) for the exact form for collinear PDFs. In Eq. (3.1) we have kept things generic and allowed for dependence on $x$ and $y$ to match the parton momenta. Since we have already converted the quasi-TMDPDF in Eq. (3) into the $\overline{\mathrm{MS}}$ scheme, $C^{\mathrm{TMD}}_{ij}$ only depends on the $\overline{\mathrm{MS}}$ scale $\mu$ and not on the lattice renormalization scale $\tilde{\mu}$. The kernel also depends on the finite proton momentum $P^{z}$ and the $\zeta$ scale that enters the TMDPDF $f_{j}^{\mathrm{TMD}}$.

Given the physical picture underlying LaMET, we can also write down the expected power corrections to the matching Eq. (3.1), obtained from the assumed hierarchy of scales

[TABLE]

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}\sim\Lambda_{\mathrm{QCD}}^{-1}$ 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 $\mathcal{O}(1)$ this also implies that $b^{z}\sim 1/P^{z}$ is small.

Note that for a matching formula like Eq. (3.1) to exist, $C^{\mathrm{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 $\zeta$ and $xP^{z}$ are independent variables, which could be achieved for example by considering a different hadronic momentum $\tilde{P}^{z}$ for the quasi-TMDPDF than the momentum $P^{z}$ used for the TMDPDF. In Eq. (3.1), the $\zeta$ dependence must then cancel between $C^{\mathrm{TMD}}_{ij}$ and $f_{j}^{\mathrm{TMD}}$ to yield a $\zeta$-independent quasi-TMDPDF $\tilde{f}_{i}^{\mathrm{TMD}}$, which implies that

[TABLE]

This clearly violates the requirement that $C^{\mathrm{TMD}}_{ij}$ must be independent of $b_{T}$ to carry out short distance matching. This mismatch occurs because the rapidity evolution governed by $\gamma_{\zeta}^{j}(\mu,b_{T})$ is also nonperturbative. To correct for this we therefore must consider the modified matching relation

[TABLE]

where for brevity we suppress the power corrections which are the same as in Eq. (3.1). In Eq. (3.1), the $\zeta$ dependence cancels between the exponential and $f^{\mathrm{TMD}}$, at the cost that the kernel $C_{ij}^{\mathrm{TMD}}$ now depends on the auxiliary scale $\tilde{\zeta}=\tilde{\zeta}(x,\tilde{P}^{z})$. As indicated, this scale must be fixed in terms of $x$ and $\tilde{P}^{z}$, the relevant scales that $\tilde{f}^{\mathrm{TMD}}$ depends on, and hence does not technically add additional functional dependence to $C_{ij}^{\mathrm{TMD}}$.

In order to interpret Eq. (3.1) as a true matching equation without any renormalization group evolution, it must be possible to make the exponential in Eq. (3.1) vanish, to yield a purely perturbative relation between $\tilde{f}^{\mathrm{TMD}}$ and $f^{\mathrm{TMD}}$. Thus, a perturbative matching can only be possible if one can choose $\tilde{\zeta}(x,\tilde{P}^{z})=\zeta$ to cancel the Collins-Soper evolution to all orders in perturbation theory. From Eq. (7) for $f_{j}^{\mathrm{TMD}}(y,\vec{b}_{T},\mu,\zeta)$ we have $\zeta=(2yP^{z})^{2}$ (choosing $y_{n}=0$ without loss of generality). Accounting for dimensions we see that we must choose $\tilde{P}^{z}\propto P^{z}$, and since at tree level $C_{ij}^{\mathrm{TMD}}[x,y,\ldots]=\delta(1-x/y)$ this fixes the constant of proportionality so that $\tilde{\zeta}=(2xP^{z})^{2}$. With $\zeta\propto y^{2}$ and $\tilde{\zeta}=\tilde{\zeta}(x)$, the only possibility that allows for perturbative matching is to take $\tilde{P}^{z}=P^{z}$ and have $C_{ij}^{\mathrm{TMD}}[x,y,\ldots]\propto\delta(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,

[TABLE]

To derive a perturbative formula for $C^{\mathrm{TMD}}_{ij}$ then requires choosing $\tilde{P}^{z}=P^{z}$ and $\zeta=\tilde{\zeta}=(2xP^{z})^{2}$. For this choice the effect of changing $\zeta$ in $f^{\mathrm{TMD}}_{j}$ is exactly balanced by a corresponding change of $P^{z}$ in $\tilde{f}_{i}^{\mathrm{TMD}}$. Note that the lack of an integral over $y$ in Eq. (3.1) is analogous to the fact that no such integral appears in the renormalization group equations for the TMDPDF, see Eq. (2.3).

We will demonstrate that a quasi-TMDPDF can be defined such that the use of Eq. (3.1) and matching with a short distance $C^{\mathrm{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.1). An all orders derivation of a formula like Eq. (3.1) would be required to completely address this relationship, and is left for future work.

3.2 Impact of finite-length Wilson lines

Before constructing quasi beam and soft functions, we discuss the impact of having infinitely-long Wilson lines in the definitions Eqs. (2.1) and (2.1) 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\to\infty$. 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 show that finite $L$ is a sufficient rapidity regulator, we first note that at leading power in $q_{T}^{2}/Q^{2}$, all emissions arise from Wilson lines and thus rapidity divergences can be regulated by modifying Wilson lines alone, as is done in most regulators in the literature. (This does not hold at subleading power, where regulating Wilson lines alone is insufficient, see Ebert:2018gsn .) Intuitively, by modifying the Wilson lines in such a way that boost invariance is broken, for example by adjusting its geometry (taking it off the light cone or restricting its length) or explicitly regulating the momentum flowing into the Wilson line, one distinguishes collinear and soft modes and thus regulates the rapidity divergences. To concretely show how this is achieved for finite $L$, consider the one-gluon Feynman rule for a Wilson line of size $L$ stretching along the $n$ direction, compared to its $L\to\infty$ limit,

[TABLE]

For $n\cdot k\to 0$, the limit where the emission of momentum $k^{\mu}$ becomes collinear to the Wilson line, i.e. it approaches infinite rapidity, is clearly regulated by exponential phase, whereas it is unregulated for $L\to\infty$. This pattern continues for multiple gluon emission since the eikonal denominators are always in one-to-one correspondence with the regulating factors involving $L$ in the numerator. The opposite limit is regulated analogously by an ${\bar{n}}$-collinear Wilson of finite length. For example, in Sec. 2 a generic example of rapidity-divergent integral was discussed, Eq. (20). For finite $L$, the example integral changes to

[TABLE]

Here we see that possible divergences as either $k^{\pm}\to 0$, corresponding to the rapidity $y=\frac{1}{2}\ln(k^{-}/k^{+})\to\pm\infty$, are regulated by having finite $L$, and the leftover logarithmic divergence as either $k^{\pm}\to\infty$ 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. (34) is

[TABLE]

Clearly, the exponentials regulate a possible divergence as $k_{z}\to 0$, and thus play a similar role as in the lightlike case. However, Eq. (35) contains a quadratic dependence on $k_{z}$ in the denominator, rather than the linear dependence on $k^{+}$ and $k^{-}$ in Eq. (34). Thus, we can also encounter linear divergences in $L$, as opposed to having only logarithmic divergences $\ln(L)$ in the lightlike case.

3.2.1 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.1), at one loop. To account for the effect of finite lattice size, the Wilson lines along the $n$ and ${\bar{n}}$ directions are truncated at $Ln$ and $L{\bar{n}}$, respectively, and transverse gauge links are included, as shown in Fig. 3(b). 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. 4(a), where a gluon connects the Wilson lines separated by the transverse displacement $\vec{b}_{T}$. Together with its mirror diagram, it is given by

[TABLE]

In the second line, one can see how keeping $L<\infty$ regulates divergences as the light-cone coordinates $k^{\pm}$ 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. 4(b) and its mirror diagram are independent of $b_{T}$, and are given by

[TABLE]

Together, we obtain

[TABLE]

where we have defined $b_{0}=2e^{-\gamma_{E}}$. The rapidity logarithms $\ln(\mu L)$ are manifest, while the last term is an example of a finite-$L$ contribution that vanishes in the limit $b_{T}\ll L$. Let us also compare this to the soft function using the $\delta$ regulator Echevarria:2012js [see also Eq. (124)]

[TABLE]

The two expression agree upon identifying $1/L^{2}=e^{2\gamma_{E}}\delta^{+}\delta^{-}$, 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$ and ${\bar{n}}L$ are given by

[TABLE]

where, ${\vec{n}}_{\perp}=\vec{b}_{T}/|\vec{b}_{T}|$. In Feynman gauge, this vertex can thus be neglected for $L\to\infty$, and one would not consider Fig. 4(c). For finite $L$, we instead obtain for Fig. 4(c)

[TABLE]

This result vanishes for $L\to\infty$ (or $b_{T}\ll L$) as expected.

The situation is more intricate for Fig. 4(d). Again, for $L\to\infty$ the diagram would not be considered. For finite $L$ instead, the $L$ dependence drops out and one obtains

[TABLE]

Here, the symmetry factor $1/2$ is compensated by the mirror diagram of Fig. 4(d). 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\to\infty$.

In order to obtain a result consistent with the known $L\to\infty$ 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^{\mathrm{TMD}}$ as in Eq. (2.1), 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_{\xi}$ gauge, but the calculation is otherwise not instructive and hence not presented here.

Full result.

The full result for the soft function with Wilson lines of finite length $L$ is given by adding Eqs. (3.2.1), (37), (41) and (3.2.1),

[TABLE]

3.3 Construction of the quasi beam function

Recall the definition Eq. (2.1) of the light-cone beam function,

[TABLE]

where $b^{\mu}=b^{+}{\bar{n}}^{\mu}/2+b_{\perp}^{\mu}$. The Wilson lines $W$ extend to light-cone infinity, where they are closed by $W_{T}$ in the transverse direction, see Fig. 1(a).

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 the $\hat{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. 5(a). The definition of the bare quasi beam function in position space thus reads

[TABLE]

where $b^{\mu}=(0,\vec{b}_{T},b^{z})$. Here, we also replaced $\gamma^{-}$ by the Dirac structure $\Gamma$, which can be chosen as either $\Gamma=\gamma^{0}$ or $\Gamma=\gamma^{z}$, as booth can be boosted onto $\gamma^{-}$. The Wilson lines of length $L$ are defined by

[TABLE]

The transverse gauge links are given by Eq. (16). Again in Eq. (45), diagrams that have no fields contracted with the states $|h(P)\rangle$ are excluded.

A crucial feature of Eq. (45) 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$ and ${\bar{n}}$ as Eq. (3.3), which makes it computable on lattice. The physical picture underlying this specific Wilson line structure is that boosting a purely spatial separation along the $\hat{z}$ direction yields an almost lightlike separation. Concretely, for a Lorentz boost along the $\hat{z}$ direction with velocity $v$ and $\gamma=1/\sqrt{1-v^{2}}$, one obtains

[TABLE]

This is illustrated in Fig. 5(b): The pure spatial separation (blue) is boosted along the orange-dotted 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 axis ${\bar{n}}^{\mu}$, as required for the $n$-collinear PDF. To boost onto $n^{\mu}$, one needs to reverse the boost, $v>0$, as is appropriate for the ${\bar{n}}$-collinear PDF, since the corresponding proton is moving into the opposite direction. It is easy to check that by applying the Lorentz boost Eq. (47) to Eq. (45), one recovers the matrix element in Eq. (3.3).

When evaluated in a large-momentum nucleon state, the quasi beam function defined in Eq. (45) is equivalent to the matrix element of an almost-lightlike correlator in a static nucleon state. According to LaMET Ji:2013dva ; Ji:2014gla , the quasi beam function is related to the beam function in Eq. (3.3) 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 Ma:2014jla ; Ma:2014jga ; Izubuchi:2018srq . For the TMDPDF the physical boost picture shown in Fig. 5(b) 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, and therefore we do not expect a simple short distance matching relation between the quasi-beam function and beam function alone. This expectation is also consistent with the known importance of the soft region in the TMD factorization theorem.

3.4 Construction of the quasi soft function

Recall the definition Eq. (2.1) of the bare TMD soft function,

[TABLE]

Note that this vacuum matrix element has no explicit time dependence, in contrast to the collinear matrix element Eq. (3.3). Time dependence only enters indirectly through the lightlike directions of the Wilson lines $S_{n}$ and $S_{{\bar{n}}}$, 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

[TABLE]

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.1) thus takes the form

[TABLE]

where the soft Wilson lines of finite length are given by

[TABLE]

The resulting Wilson line path is illustrated in Fig. 6(a).

Unfortunately, the physical boost argument that allowed us to relate spatial Wilson lines to lightlike Wilson lines in the quasi-PDF [see Eq. (47)] does not apply to the quasi soft function. Since the soft function involves both light-cone directions $n$ and ${\bar{n}}$, it is necessary to simultaneously obtain them from boosting $\pm\hat{z}$. However, this requires boosts of opposite signs, as illustrated in Fig. 6(b). Note that if this were possible with a single boost, it would directly violate the boost argument for $\tilde{B}$, 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.1) 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.

4 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 $\overline{\mathrm{MS}}$ scheme, as opposed to considering renormalization schemes appropriate for lattice calculations as discussed in Sec. 3, since a pure $\overline{\mathrm{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 non-singlet 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^{\mu}=(P^{z},0,0,P^{z})$.

4.1 Lightcone TMDPDF at one loop

The unrenormalized result for the TMDPDF at one loop is given by

[TABLE]

where $P_{qq}(x)=(1+x^{2})/(1-x)$ is the quark-to-quark splitting function and the subscript + denotes a plus distribution such that $\int_{0}^{1}\mathrm{d}x\,[f(x)]_{+}^{1}=0$. 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}\ll\Lambda_{\mathrm{QCD}}^{-1}$, see Eq. (17). 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}\sim\Lambda_{\mathrm{QCD}}^{-1}$ to exist.

4.2 Quasi beam function

We first calculate the quasi beam function defined in Eq. (45) 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\to\infty,P^{z}\to\infty$, 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 by

[TABLE]

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. (35), 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 combining $\tilde{B}_{q}$ and $\tilde{\Delta}_{S}^{q}$ as in Eq. (3).

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 $-\frac{1}{2}\ln^{2}(b_{T}^{2}\mu^{2}/b_{0}^{2})$, 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_{T}^{2}\mu^{2}/b_{0}^{2})$ does not match up with the corresponding single logarithm in the quasi beam function.555We have also checked that this problem is not simply due to the contribution from the transverse Wilson line self energy diagram. 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.

4.3 Naive quasi soft function

Next, we calculate the naive quasi soft function defined in Eq. (3.4). 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 $+\hat{z}$ and $-\hat{z}$ Wilson lines (upper row) and diagrams where the gluon is emitted between Wilson lines of the $+\hat{z}$ sector alone (lower row). The latter is identical to the result for gluons exchanged within the $-\hat{z}$ sector.

In Feynman gauge, the generic expression for a one-loop diagram in the $\overline{\mathrm{MS}}$ scheme, parameterized by the spatial paths $\gamma_{1}$ and $\gamma_{2}$ of the Wilson lines, is

[TABLE]

Note that diagrams with the gluon attaching to the same line have an additional symmetry factor of $1/2$.

For example, for Fig. 8(a) one reads off the paths

[TABLE]

Together with its mirror diagram, this gives

[TABLE]

Note that this result also contains a divergence at $\epsilon=1/2$, which signals a power-law divergence $\propto\mu 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 $\epsilon=0$, but on the lattice they explicitly arise and have to be canceled. For this reason, the beam function renormalization in Eq. (3) contains a $b^{z}$-dependent counter term to absorb the divergence associated with the self energy of the Wilson-line segment of length $L-b^{z}$ [see Eq. (C.3)], while the other self-energies will cancel against the soft factor.

The diagrams in Figs. 8(b) and 8(c) yield

[TABLE]

Summing Eqs. (4.3)–(58), we obtain the soft contribution from interactions between the two collinear sectors,

[TABLE]

where we have also taken the limit $L\gg 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. 8(d)–8(f) yield

[TABLE]

and their sum gives the contribution to the soft function from interactions within the $\hat{z}$ sector alone. The same result is obtained for the $-\hat{z}$ sector, so we have

[TABLE]

Note that this contains the same linear divergence in $L/b_{T}$ as the quasi beam function, Eq. (53). 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.3) and (4.3),

[TABLE]

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.

4.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). To fix the form of the soft factor $\tilde{\Delta}_{S}^{q}$, we require that all divergences in $L/b_{T}$ cancel when combining $\tilde{B}_{q}$ and $\tilde{\Delta}_{S}^{q}$, analogous to the cancellation of rapidity divergences in the TMDPDF. Comparing the one-loop results in Eqs. (53) and (4.3), we deduce this to be $\tilde{\Delta}_{S}^{q}=1/\sqrt{\tilde{S}^{q}}$, so

[TABLE]

Note that this is similar to the $\delta$ regulator for the TMDPDF, where one has $\Delta_{S}^{q}=1/\sqrt{S^{q}}$. Our result $\tilde{\Delta}_{S}^{q}=1/\sqrt{\tilde{S}^{q}}$ 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 $\delta$ regulator.

Here, we work purely in the $\overline{\mathrm{MS}}$ scheme and in the physical limit, where the product $\tilde{Z}^{\prime}\tilde{Z}_{\rm uv}$ is $b^{z}$-independent, so for our particular one-loop study we can equivalently work with the longitudinal momentum space formula

[TABLE]

Combining Eqs. (53) and (4.3) according to Eq. (66), we obtain the NLO result for the quasi-TMDPDF evaluated in an on-shell quark state,

[TABLE]

Although our method of calculation is quite different, we note that our result in Eq. (67) fully agrees with the one loop calculation in Ref. Ji:2018hvs (up to trivial differences in our respective conventions for the $\overline{\rm MS}$ scheme).

There is an important subtlety concerning the logarithmic dependence $\ln(2P^{z})^{2}$ in Eqs. (4.1) and (67), which arises from calculating matrix elements with an on-shell external quark of momentum $P^{\mu}=(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}\to xP^{z}$, where $xP^{z}$ is the momentum of the struck parton (see also the discussion in Ref. Izubuchi:2018srq ), so we obtain

[TABLE]

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 between $\tilde{f}_{q}^{\mathrm{TMD}}$ and $f_{q}^{\mathrm{TMD}}$ is possible, we need to UV renormalize both Eqs. (4.1) and (4.4) and study their difference:

[TABLE]

As expected, the explicit infrared poles in $\epsilon_{\mathrm{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 of $\tilde{f}_{q}^{\mathrm{TMD}}$ and $f_{q}^{\mathrm{TMD}}$ does not agree, leaving two uncanceled logarithms in Eq. (69). The second one multiplies a logarithm of $\zeta$, and in fact is exactly the one-loop expansion of the Collins-Soper kernel in Eq. (26),

[TABLE]

This confirms our argument in Sec. 3.1 that the Collins-Soper equation prohibits a perturbative matching, unless $\zeta$ is fixed in terms of $P^{z}$ such that the $\zeta$-dependence of $\tilde{f}_{q}^{\mathrm{TMD}}(x,\vec{b}_{T},\tilde{\mu},P^{z})$ and $f_{q}^{\mathrm{TMD}}[x,\vec{b}_{T},\mu,\zeta(P^{z})]$ exactly cancel. From Eq. (69), this is fulfilled with

[TABLE]

as expected. This leaves

[TABLE]

The key problem with Eq. (4.4) is that it still contains a single infrared logarithm $\ln(b_{T}^{2}\mu^{2}/b_{0}^{2})$ which is not associated with the Collins-Soper evolution, and thus cannot be eliminated in a similar fashion. Curiously, choosing $\zeta=(2xP^{z})^{2}/e$ would simultaneously cancel both logarithms in $b_{T}$ in Eq. (69), but the term involving a $\ln\zeta$ is clearly related to the Collins-Soper kernel, whereas there is no clear relationship of the term $\ln(b_{T}^{2}\mu^{2}/b_{0}^{2})$ 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_{T}^{2}\mu^{2}/b_{0}^{2})$ 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. Ji:2018hvs , 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 $\mu^{2}=\zeta=(2xP^{z})^{2}$ which was given as

[TABLE]

If we take our result in Eq. (69) and set $\mu^{2}=\zeta=(2xP^{z})^{2}$, then we obtain

[TABLE]

This agrees with expanding Eq. (4.4) to ${\cal O}(\alpha_{s})$, where the single infrared logarithm is generated by the exponential term. (There is a trivial mismatch from the constant $\pi^{2}/12$ term which arises because Ref. Ji:2018hvs uses a different definition of the $\overline{\mathrm{MS}}$ scheme, see appendix A.) In Ref. Ji:2018hvs Eq. (4.4) was interpreted as being a valid matching formula between quasi-TMDPDF and TMDPDF. However, for nonperturbative $b_{T}$ the exponential in Eq. (4.4) becomes a nonperturbative function and cannot be included in a short distance matching coefficient, in agreement with our conclusions.

We conclude this section with an overview of the dependence of the one-loop coefficients of beam function $B_{q}$, soft subtraction $\Delta_{S}^{q}$ and TMDPDF $f_{q}^{\mathrm{TMD}}$ on the logarithm $L_{b}=\ln(b_{T}^{2}\mu^{2}/b_{0}^{2})$ for the different regulators in the literature, as shown in Table 1. The dependence of the quasi constructions $\tilde{B}_{q}$, $\tilde{\Delta}_{S}^{q}$ and $\tilde{f}_{q}^{\mathrm{TMD}}$ 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_{b}^{2}$ 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).666Note that in none of these cases does one include the transverse self energy, which would add a single logarithm $L_{b}$ to $B_{q}$ and $-L_{b}$ to $\Delta_{S}^{q}=1/\sqrt{S^{q}}$, see Eq. (3.2.1). Even after taking this into account, the mismatch persists.

4.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$) 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 $\hat{z}$ or $-\hat{z}$ sector, or between them, as done in Sec. 4.3,

[TABLE]

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 towards ${\bar{n}}/{\bar{n}}$ for $P^{z}<0$. In practice, they are identical due to invariance under $z\leftrightarrow-z$. At one loop, the quasi-TMDPDF in Eq. (66) hence can be written as

[TABLE]

Next, recall the Wilson line structures of quasi beam and quasi soft functions in Eqs. (45) and (3.4), see also Figs. 5(a) and 6(a). Taking the soft limit of the quasi beam function, $b^{z}=0$, clearly gives the same Wilson lines as half of the quasi soft function, and hence subtracting $\tilde{S}_{z/z}^{(1)}$ exactly cancels the soft limit of the tadpole correction to the beam function. This can easily be verified by comparing Eqs. (4.3) and (203), from which one also sees that this subtraction cancels the divergence in $L/b_{T}$. It remains to consider $\tilde{S}_{z/-z}^{(1)}$, given in Eq. (4.3). Its subtraction from $\tilde{B}_{q}$ adds a single $\ln(b_{T}^{2}\mu^{2}/b_{0}^{2})$, which is exactly the leftover logarithm found in the relation between the naive quasi-TMDPDF and TMDPDF in Eq. (69). In conclusion, it thus appears to be precisely the interaction between the $+\hat{z}$ and $-\hat{z}$ 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 TMDPDF. 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 as

[TABLE]

where ${\bar{n}}_{\perp}^{\mu}$ is the transverse unit vector orthogonal to $n_{\perp}^{\mu}=b_{\perp}^{\mu}/b_{T}$ and $\hat{z}$. Fig. 9 illustrates the Wilson line path in Eq. (4.5) and compares it to the path for naive quasi soft function defined in Eq. (3.4).

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 $+\hat{z}$ and $-\hat{z}$ directions. These diagrams precisely vanish for the bent soft function due to ${\bar{n}}_{\perp}\cdot n_{\perp}={\bar{n}}_{\perp}\cdot\hat{z}=0$, while all other diagrams are not affected by the new Wilson line paths. Hence the bent soft function at one loop yields

[TABLE]

As before, it is related to the soft subtration through $\tilde{\Delta}_{S}^{q}=1/\sqrt{\tilde{S}_{\rm bent}}$. This bent soft factor has precisely the infrared logarithms that are desired at one loop.777 A soft factor with more than one transverse directions was also used in Ref. Ji:2014hxa where the goal was to express the physical factorization theorem directly in terms of quasi objects. They define a soft factor

$\displaystyle\tilde{\Delta}_{S}^{q\,{\rm JSXY}}(b_{T})=\sqrt{\frac{\tilde{S}_{n_{x},n_{y}}(b_{T})}{\tilde{S}_{n_{x},n_{z}}(b_{T})\tilde{S}_{n_{z},n_{y}}(b_{T})}}\,,$

(79)

where $\tilde{S}_{n_{1},n_{2}}(b_{T})$ is the same as the result obtained from our naive quasi soft matrix element, Eq. (3.4), by replacing $\hat{z}\to n_{1}$ and $-\hat{z}\to n_{2}$ but using infinitely-long Wilson lines ($L\to\infty$). Taking $b_{T}^{\mu}=b_{T}n_{x}^{\mu}$ and using instead finite Wilson lines, we have checked that the resulting combination of terms in $\tilde{\Delta}_{S}^{q}(b_{T},L)$ gives the same result at one-loop as that of our bent soft function.

Using Eq. (66) to combine the bent quasi soft function from Eq. (78) together with the natural quasi beam function from Eq. (53) we obtain a new quasi-TMDPDF

[TABLE]

Comparing this result to the TMDPDF at one loop yields

[TABLE]

where have again fixed $\zeta=(2xP^{z})^{2}$ as explained previously. Since there is no $b_{T}$ dependence on the RHS of Eq. (4.5), 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.1) with a one loop result for the matching coefficient that is given by

[TABLE]

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 $\delta_{qq^{\prime}}$. 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. (82) 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 $\overline{\rm 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. Constantinou:2017sej ; Stewart:2017tvs ).

5 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.

5.1 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 $-\hat{z}$ and $+\hat{y}$ 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 3our 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 non-singlet (${\rm ns}$) quark flavor combination where there is no mixing reads:

[TABLE]

If the infrared structure of $\tilde{f}_{{\rm ns}}^{\mathrm{TMD}}(x,\vec{b}_{T},\mu,P^{z})$ and $f^{\mathrm{TMD}}_{{\rm ns}}(x,\vec{b}_{T},\mu,\zeta=(2xP^{z})^{2})$ match to all orders, then our one loop calculation in Sec. 4.5 provides a valid result for $C^{\mathrm{TMD}}_{{\rm ns}}$. For the matching from the $\overline{\mathrm{MS}}$ renormalized quasi-TMDPDF it gives

[TABLE]

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.1) 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 $\tilde{\Delta}_{S}^{q}$ and $\Delta_{S}^{q}$ 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:

[TABLE]

Here, $C_{\rm ns}^{\mathrm{TMD}}$ is by definition still a perturbative function. On the other hand, the function $g^{S}_{q}(b_{T},\mu)$ 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},\mu)$’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.1) 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

[TABLE]

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_{q}^{S}$ factors cancel out. We discuss this further in the next section.

5.2 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_{i}^{S}(b_{T},\mu)$ only differs for quarks and gluons but otherwise is flavor blind to try to construct ratios of quasi-TMDPDFs where $g_{i}^{S}(b_{T},\mu)$ cancels. To avoid possible mixing between quarks and gluons, here we only consider nonsinglet quasi-TMDPDFs such as ${\rm ns}=u{-}d$, where Eq. (5.1) applies. The soft factor $g_{q}^{S}$ cancels in any ratio of two (quasi) TMDPDFs with the same choice of $b_{T}$ and $\mu$, so we have

[TABLE]

For example, one can choose $x_{1}=x_{2}$ and $\zeta_{1}=\zeta_{2}$ to expose the Collins-Soper kernel $\gamma_{\zeta}^{q}$ as

[TABLE]

This allows one to determine the nonperturbative $b_{T}$ dependence of $\gamma_{\zeta}^{q}$ from the ratio of the quasi-TMDPDFs $\tilde{f}^{\mathrm{TMD}}_{\rm ns}$ computed in lattice QCD, and was proposed in Ref. Ebert:2018gzl .888In the context of lattice calculations, it is also worth noting that $\gamma_{\zeta}^{q}$ depends on the quark masses in the discretized action used in the computation, however, due to its independence of the hadronic state it is independent of the valence quark masses used to construct correlation functions. We thank Phiala Shanahan and Michael Wagman for discussions on this point.

Without this knowledge of $\gamma_{\zeta}^{q}$, one is forced to choose $x_{1}P_{1}^{z}=x_{2}P_{2}^{z}$ and $\zeta_{1}=\zeta_{2}$ in Eq. (5.2) to cancel the Collins-Soper evolution kernel,

[TABLE]

Here, $h_{1}$ and $h_{2}$ for example refers to TMDPDF in different hadron states (which does not affect $C^{\mathrm{TMD}}_{{\rm ns}/h_{i}}$), or TMDPDFs of different spin structures (which can affect $C^{\mathrm{TMD}}_{{\rm 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. Collins:1350496 . 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. (89), confirming that it indeed is a matching equation at this order.

In Eqs. (5.2)–(89) we write ratios of quasi-TMDPDFs, which contains some choice for the quasi soft factor that satisfies Eq. (5.1). However, one can completely avoid the need for including a quasi soft factor in such ratios by employing Eq. (3), which yields

[TABLE]

In this ratio, the quasi soft factor cancels because it is independent of $b^{z}$ and $P_{1,2}^{z}$. Thus Eq. (90) removes the necessity to calculate a quasi soft matrix element on lattice. The leftover divergences as $L\to\infty$ present in the individual $\tilde{B}_{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. Engelhardt:2015xja ; Yoon:2016dyh ; Yoon:2017qzo . 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. (89). 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}\to 0$ corresponding to integrating over $x$. The relations they use do not require a non-trivial matching coefficient $C^{\mathrm{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.

6 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 $\Delta_{S}^{q}$ 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 function $\tilde{B}_{q}$ and the quasi soft factor $\tilde{\Delta}_{S}^{q}$. 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^{\mu}\to\hat{z}^{\mu}$ and ${\bar{n}}^{\mu}\to-\hat{z}^{\mu}$. Both matrix elements separately suffer from rapidity divergences that are regulated by the length $L$, and we construct the naive quasi-TMDPDF $\tilde{f}_{q}^{\mathrm{TMD}}=\tilde{B}_{q}\tilde{\Delta}_{S}^{q}$ by demanding the cancellation of the $L$ dependence. We also discussed implications of the Collins-Soper evolution on the matching relation between $f_{q}^{\mathrm{TMD}}$ and $\tilde{f}_{q}^{\mathrm{TMD}}$, which (if it exists) requires one to fix the Collins-Soper scale $\zeta$ 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}\sim q_{T}^{-1}\sim\Lambda_{\mathrm{QCD}}^{-1}$, 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 $\tilde{B}_{q}$ and $\tilde{\Delta}_{q}^{S}$ into the quasi-TMDPDF $\tilde{f}_{q}^{\mathrm{TMD}}$, 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 $\tilde{\Delta}_{S}^{q}$ 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 within LaMET works beyond one loop order then it, for the first time, provides a definite 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. Ebert:2018gzl , 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.

Acknowledgements.

We thank W. Detmold, M. Diehl, M. Engelhardt, P. Shanahan, and M. Wagman for discussions. This work was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, from DE-SC0011090 and within the framework of the TMD Topical Collaboration. I.S. was also supported in part by the Simons Foundation through the Investigator grant 327942. M.E. was also supported by the Alexander von Humboldt Foundation through a Feodor Lynen Research Fellowship. M.E. and I.S. also thank KIT for hospitality while portions of this work were completed.

Appendix 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

[TABLE]

which satisfy $n^{2}={\bar{n}}^{2}=0$ and $n\cdot{\bar{n}}=2$. One can decompose any fourvector as

[TABLE]

where the plus and minus components are defined as

[TABLE]

and the transverse component is orthogonal to both reference vectors, $n\cdot k_{\perp}={\bar{n}}\cdot k_{\perp}=0$. The Minkowski product of two fourvectors follows to be

[TABLE]

In particular, one has $k^{2}=k^{+}k^{-}-k_{T}^{2}$.

Another convention often used in the literature is $k^{\pm}=(k^{0}\pm k^{z})/\sqrt{2}$, see e.g. Ref. Collins:1350496 . In that notation, the Minkowski product is x\cdot k=\bigl{(}x^{+}k^{-}+x^{-}k^{+}\bigr{)}-\vec{x}_{T}\cdot\vec{k}_{T}, and $k^{2}=2k^{+}k^{-}-k_{T}^{2}$. To translate results from this paper to that convention, replace $k^{\pm}\to\sqrt{2}k^{\mp}$. Often, this leaves factor $\sqrt{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^{\mu}=(b^{0},\vec{b}_{T},b^{z})$. The Fourier transform to and from momentum space is defined as

[TABLE]

We use the same symbol for a function in $f$ in both spaces, as we reserve the symbol $\tilde{f}$ 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

[TABLE]

Renormalization scheme.

Our results are expressed in the $\overline{\mathrm{MS}}$ scheme. The associated renormalization scale $\mu\equiv\mu_{\overline{\mathrm{MS}}}$ is related to the MS scale $\mu_{0}$ by

[TABLE]

Note that this differs from the convention used by Collins Collins:1350496 , where $\mu^{2}=4\pi\mu_{0}^{2}/\Gamma(1-\epsilon)$, by terms of $\mathcal{O}(\epsilon^{2})$.

Appendix 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.1) as

[TABLE]

and the Collins-Soper scale was given by Eq. (7),

[TABLE]

We show how this formulation for the TMDPDF arises in the various regulators, including Wilson lines off the light-cone, the $\delta$ regulator, the $\eta$ regulator, and the exponential regulator from Refs. Collins:1350496 ; GarciaEchevarria:2011rb ; Chiu:2012ir ; Li:2016axz 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_{q}^{\mathrm{TMD}}$. We also give the correspondence to results for the analytic regulator of Ref. Becher:2010tm , where only the product of two TMDs are defined, but which also agree with this $f_{q}^{\mathrm{TMD}}$.

At one-loop, the common result when evaluating the TMDPDF in an on-shell quark state of momentum $P$ using pure dimensional regularization is

[TABLE]

Here, $L_{b}$ is defined in Eq. (96) 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

[TABLE]

one obtains the renormalized TMDPDF as

[TABLE]

The remaining $1/\epsilon_{\mathrm{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 Collins:1350496 , the lightlike Wilson lines are replaced by spacelike Wilson lines,

[TABLE]

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. Collins:1350496 we have for the $n$-collinear TMDPDF

[TABLE]

where the result in the last line was derived in Ref. Buffing:2017mqm . 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.1) 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 $\zeta$ scale is defined as in Eq. (99). To relate Eq. (B.1) to our notation in Eq. (B.1) one identifies $\tau=1/(y_{B}-y_{n})$ and notes that

[TABLE]

so that the translation of functional forms is

[TABLE]

where $B_{q}^{C}\equiv\tilde{f}^{\rm unsub}_{q}$ is the collinear matrix element as defined in Eq. (13.108) in Ref. Collins:1350496 and $S^{q}_{C}\equiv\tilde{S}_{(0)}$ is the soft function as defined in Eq. (13.39) therein. As usual, the dependence on $m_{P}$ and $\Lambda_{\mathrm{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/\tau-\ln\sqrt{\zeta}$ appears in $B_{q}$.

Since $\zeta=(xm_{p}e^{y_{P}-y_{n}})^{2}$, one can obtain the Collins-Soper kernel as

[TABLE]

which is the analog of Eq. (23). This makes clear that the Collins-Soper kernel can be obtained as a differential equation involving $f^{\mathrm{TMD}}$, or $B_{q}$, or $\Delta_{S}^{q}$ 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_{q}^{C}$ and $S_{C}^{q}$ in the literature. Instead, in Refs. Aybat:2011zv ; Collins:1350496 the two functions are combined at the integrand level before integrating over the momentum of either real or virtual emissions. For example, in Ref. Collins:1350496 , where the soft factor is calculated in the context of the TMD fragmentation function, the limit $y_{A,B}\to\pm\infty$ 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. (102). The same order of operations was used for the calculation of the TMDPDF for double-parton scattering in Ref. Buffing:2017mqm .

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

[TABLE]

Note that integral in in Eq. (B.1) is ill-defined due to a pinched pole singularity arising from the opposite signs of the $\mathrm{i}0$ prescription in the two eikonal vertices, and we have evaluated it using the principal-value prescription. This is similar to the calculation in Ref. Bacchetta:2008xw , where the soft function as defined in Ref. Collins:1981uk was calculated in an axial gauge with $v\cdot\mathcal{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. Bacchetta:2008xw . In the calculation of Ref. Ji:2004wu , 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}|\gg 1$, we obtain

[TABLE]

which after renormalization agrees with the result for timelike Wilson lines in Ref. Ji:2004wu . Thus at one loop using Eq. (B.1) with $\tau=1/(y_{B}-y_{n})$ we find

[TABLE]

The one-loop result for $f_{q}^{\mathrm{TMD}}$ was calculated in Ref. Aybat:2011zv and agrees with Eq. (102) after adjusting for the difference in definitions of the $\overline{\mathrm{MS}}$ scheme, see appendix A. Hence we can also deduce the bare beam function in this scheme, using

[TABLE]

Using Eq. (105), one can see that this result only depends on $y_{P}-y_{B}$.

B.2 $\delta$ regulator

The $\delta$ regulator was originally introduced in Ref. GarciaEchevarria:2011rb and then modified in Refs. Echevarria:2015usa ; Echevarria:2015byo to be applicable for NNLO calculations. Note that their notation for the light-cone momenta $x^{\pm}$ are related to our notation by $x^{\pm}\leftrightarrow x^{\mp}$. For consistency of this paper, we convert results for the $\delta$ regulator in the literature to our notation.

Here, we follow the definition of the TMDPDF as presented in Ref. Echevarria:2012js , which generalizes the definition in Ref. GarciaEchevarria:2011rb . At one loop, the $\delta$ regulator amounts to modifying propagators in the full theory as

[TABLE]

where

[TABLE]

are the momenta of the struck quark in the $n$-collinear and ${\bar{n}}$-collinear proton, respectively. One can equivalently modify the eikonal propagators as

[TABLE]

(Additional integer $m>0$ multiplicative factors, $k^{\pm}+\mathrm{i}\,m\,\delta^{\pm}$, appear at higher loop orders). The regulators are related by

[TABLE]

When working with $\Delta^{\pm}$ as arising in the full theory, the soft function can be written as Echevarria:2012js

[TABLE]

However, the $Q^{2}$ dependence appearing in $S$ here is artificial, as it is the $\delta^{\pm}$ regulators that are the fundamental Wilson line regulators. Using Eq. (116) together with $\bar{p}^{+}p^{-}=Q^{2}$, one can equivalently write this as

[TABLE]

where $y_{n}$ is an arbitrary parameter that cancels between the two logarithms. Eq. (118) allows one to split the soft function into a “left-moving” and “right-moving” component,999 Eq. (120) might seem to differ from Eq. (16) in Ref. Echevarria:2012js , where the soft function is expressed as

\displaystyle S\biggl{(}b_{T},\frac{\Delta^{-}}{p^{+}},\frac{\Delta^{+}}{\bar{p}^{-}}\biggr{)}=\sqrt{S\biggl{(}b_{T},\frac{\Delta^{-}}{\bar{p}^{+}},\alpha\frac{\Delta^{-}}{p^{-}}\biggr{)}}\sqrt{S\biggl{(}b_{T},\frac{1}{\alpha}\frac{\Delta^{+}}{\bar{p}^{+}},\frac{\Delta^{+}}{p^{-}}\biggr{)}}\,.

(119)

Since $S$ only depends on the product of the last two arguments, one can identify $\alpha=e^{2Y-2y_{n}}$ and use that $\alpha(\Delta^{-})^{2}/Q^{2}=\alpha(\delta^{-}e^{-Y})^{2}=(\delta^{-}e^{-y_{n}})^{2}$ and $(\Delta^{+})^{2}/(\alpha Q^{2})=(\delta^{+}e^{Y})^{2}/\alpha=(\delta^{+}e^{y_{n}})^{2}$ to show that Eq. (119) is equivalent to our Eq. (120).

[TABLE]

and it is this form of $S^{q}_{\rm EIS}$ that was used in NNLO calculations in Refs. Echevarria:2015usa ; Echevarria:2015byo ; Echevarria:2016scs . The two factors in Eq. (120) are then combined with the $n$-collinear and ${\bar{n}}$-collinear matrix elements, respectively. Since with the $\delta$ 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 $\sqrt{S}$. For the $n$-collinear case, one identifies

[TABLE]

where $B_{q}^{\rm EIS}\equiv J_{n}$ appears in Eq. (12) in Ref. Echevarria:2012js . Once again, only the combination

[TABLE]

arises in $B_{q}$, such that the functional dependence of $B_{q}$ and $J_{n}$ in Eq. (B.2) agree.

One-loop results.

The unsubtracted beam function has first been calculated in Ref. GarciaEchevarria:2011rb , where the $\delta$ regulator also acts as an IR regulator. The corresponding result with dimensional regularization as the IR regulator can be extracted from Ref. Echevarria:2016scs ,

[TABLE]

The one-loop soft function is given by Echevarria:2012js

[TABLE]

According to Eqs. (B.2) and (B), the TMDPDF at NLO follows as

[TABLE]

which reproduces the result in Eq. (100) with $\zeta$ as given in Eq. (99).

B.3 $\eta$ regulator

The $\eta$ regulator Chiu:2011qc ; Chiu:2012ir modifies collinear and soft Wilson lines in momentum space as

[TABLE]

where the momentum operator $\mathcal{P}$ picks up the momentum of the gluon fields $A$, and $\eta$ is the rapidity regulator with an associated rapidity scale $\nu$. While the authors of Ref. Chiu:2012ir separately renormalize the beam and soft functions in $\eta$, one can trivially define a TMDPDF in the form of Eq. (B) by combining the bare beam and soft functions, in which case the dependence on $\eta$ and $\nu$ cancels,

[TABLE]

where $B_{q}\equiv B_{q}^{\rm CJNR}$ and $S^{q}\equiv S^{q}_{\rm CJNR}$ are the beam and soft function as defined in Ref. Chiu:2012ir . For this regulator, the soft zero-bin subtractions vanish, so $\Delta_{S}^{q}=\sqrt{S^{q}_{\rm CJNR}}$ only involves a $\sqrt{S^{q}}$ in the numerator. For this regulator the correspondence with our notation is

[TABLE]

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. (128), we combined bare beam and soft functions as usual for TMDPDFs. However, in practice the $\eta$ regulator is often used to define rapidity-renormalized beam and soft functions by first expanding them individually about $\eta\to 0$ and $\epsilon\to 0$ and then absorbing poles in $\eta$ and $\epsilon$ through separate $\overline{\mathrm{MS}}$ counterterms Chiu:2012ir ,

[TABLE]

When combined these counterterms give

[TABLE]

Using the fact the bare $B_{q}^{\rm CJNR}$ and $S^{q}_{\rm CJNR}$ are $\mu$ and $\nu$ independent yields the $\mu$ RGEs

[TABLE]

and the rapidity RGEs

[TABLE]

Using Eqs. (128), (B.3) and (131) gives the renormalized TMDPDF in terms of the renormalized beam and soft functions as

[TABLE]

where the $\nu$ independence cancels on the right hand side, and we remind the reader that here $\zeta=(xP^{-})^{2}$. Since $B_{q}^{\rm CJNR}$ only depends on the ratio $xP^{-}/\nu$, the Collins-Soper kernel can be extracted as

[TABLE]

This also shows that the rapidity anomalous dimension of Ref. Chiu:2012ir is related to the Collins-Soper kernel by $\gamma_{\nu}^{q}(\mu,b_{T})=2\gamma_{\zeta}^{q}(\mu,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. Chiu:2012ir . We find

[TABLE]

The NLO soft function using this regulator can be found in Refs. Chiu:2012ir ; Luebbert:2016itl ,

[TABLE]

Combining the bare functions as in Eq. (128) gives the NLO TMDPDF as

[TABLE]

which yields Eq. (100) with $\zeta=(xP^{-})^{2}$. Here, the poles in $\eta$ manifestly cancel.

Alternatively, one can separately UV and rapidity renormalize beam and soft functions by absorbing the poles in $\eta$ and $\epsilon$ in separate counterterms. The renormalized functions are given by

[TABLE]

where the remaining $1/\epsilon$ pole in Eq. (B.3) is of infrared origin. Combining these as in Eq. (134) yields the UV-renormalized TMDPDF Eq. (102) with $\zeta=(xP^{-})^{2}$.

B.4 Exponential regulator

The exponential regulator Li:2016axz has been designed to connect the TMD soft function to the threshold soft function, which enabled it to be calculated up to ${\cal O}(\alpha_{s}^{3})$ Li:2016ctv . The regulator introduces a factor

[TABLE]

into the phase space of real emissions. One then takes the $\tau\to 0$ limit, keeping only divergent terms. Identifying $\tau=\nu^{-1}$, one can connect this to the $\eta$ regulator shown above. Here, we will keep $\tau$ for clarity.

Similar to the $\eta$ 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

[TABLE]

one obtains the TMDPDF. Similar to our discussion for the $\eta$ 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. HuaXing ,

[TABLE]

The renormalized NLO soft function in this scheme is given by Li:2016axz

[TABLE]

Combining the renormalized functions gives the NLO TMDPDF as

[TABLE]

which reproduces Eq. (102) with $\zeta=Q^{2}$.

B.5 Analytic regulator

The analytic regulator as introduced in Ref. Becher:2010tm for TMDs raises propagators including $n$ and ${\bar{n}}$-collinear momenta to a power,101010In practice, one can also introduce a factor $(\nu/n\cdot k)^{\alpha}$ for each unresolved final-state parton with momentum $k$, as suggested in Ref. Becher:2011dz and employed in the NNLO calculation in Ref. Gehrmann:2014yya .

[TABLE]

where $\alpha$ and $\beta$ are distinct parameters. One then has to consistently expand first in $\beta\to 0$ and then in $\alpha\to 0$, or vice versa. Since this breaks the $n\leftrightarrow{\bar{n}}$ symmetry, it regulates all rapidity divergences. However, the regulator does not act symmetrically in the $n$ and ${\bar{n}}$-collinear beam function and thus yields two distinct bare quantities $\mathcal{B}_{q\,n}^{\rm BN}$ and $\mathcal{B}_{q\,{\bar{n}}}^{\rm BN}$. Following Ref. Becher:2010tm , we give all results for expanding first in $\beta\to 0$ and then in $\alpha\to 0$, such that rapidity divergences are regulated by $\alpha$ alone.

A particular feature of this regulator is that the soft function is unity all orders,

[TABLE]

as all loop corrections are scaleless and vanish. This implies that rapidity divergences, i.e. poles in $1/\alpha$, cannot cancel between a single beam function and the soft factor, and hence one cannot define a rapidity-finite TMDPDF as in Eq. (B). Instead, divergences in $\alpha$ only cancel in the product, so one has to define the TMDPDFs through the limit

[TABLE]

Here, $F_{q\bar{q}}(\mu,b_{T})=-\gamma^{q}_{\zeta}(\mu,b_{T})$ is the so-called collinear anomaly coefficient Becher:2010tm , which exponentiates the rapidity logarithms. Note that the TMDPDFs on the right hand side of Eq. (B.5) do not have an explicit $\zeta$ dependence, as $\zeta_{a}=\zeta_{b}=Q$ has already been chosen to fully exponentiate all rapidity logarithms.

One-loop results.

The one-loop beam functions have been calculated in Becher:2010tm . As discussed, the $n$ and ${\bar{n}}$-collinear functions take different forms due to the regulator breaking the $n\leftrightarrow{\bar{n}}$ symmetry. The one-loop result for the matrix element in an on-shell quark state can be extracted from Ref. Becher:2010tm in the form

[TABLE]

Only the product of the two beam functions is meaningful, as the poles in $\alpha$ and the $\nu$ dependence cancel, and one obtains the product of the two TMDPDFs,

[TABLE]

One can easily check that the same is obtained for $f_{q}^{\mathrm{TMD}}(x_{1},\vec{b}_{T},\epsilon,\zeta_{a})f_{q}^{\mathrm{TMD}}(x_{2},\vec{b}_{T},\epsilon,\zeta_{b})$ using the result in Eq. (100) together with $\zeta_{a}\zeta_{b}=Q^{4}$.

Appendix C Calculation of the quasi beam function

We calculate the spin-averaged quasi beam function by evaluating Eq. (45) with an on-shell quark with momentum $P^{\mu}=(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,

[TABLE]

where for compactness we suppress the dependence of $\tilde{q}_{n}$ on $\epsilon$ and $L$. Here $\lambda=0$ or $\lambda=3$, and the path ordered exponential represents the Wilson line stretching along the path $\gamma$, as illustrated in Fig. 5. The explicit parameterization is given by three segments

[TABLE]

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

[TABLE]

Here, $K_{n}(x)$ is the modified Bessel function of the second kind. We need only the following integrals over the Feynman parameter $y$,

[TABLE]

where $\mu_{0}$ is the MS scale and $\mu$ is the $\overline{\mathrm{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\in[-\infty,\infty]$, rather than the physical support $x\in[0,1]$. For the calculation, it is most convenient to define the plus distribution such that the integral over all $x$ vanishes,

[TABLE]

For functions with physical support $x\in[0,1]$, this naturally reduces to the standard plus distribution. Any integrable expression $h(x)$ can be written as a distribution using

[TABLE]

C.1 Vertex correction

The vertex correction Fig. 11(a) is given by

[TABLE]

Note that the second term in square brackets vanishes for $\lambda=0$, as the integrand is odd in $k^{0}$. The Fourier transform w.r.t. $b^{z}$ can be carried out trivially,

[TABLE]

Introducing Feynman parameters, we obtain

[TABLE]

Using Eqs. (C) and (159), this can be expressed in terms of the master integrals,

[TABLE]

Next we use Eq. (163) to rewrite Eq. (167) in terms of a plus distribution. Note that the integral of the second line in Eq. (C.1) vanishes when integrated over all $x$ and thus does not yield a boundary term. We obtain

[TABLE]

Plugging in the master formulas and expanding in $\epsilon$ gives the final result

[TABLE]

where we singled out the terms with physical support $x\in[0,1]$. These will be crucial to recover the collinear singularity of the standard PDF. The contribution with unphysical support $x\in[-\infty,\infty]$ is given by

[TABLE]

The first line in Eq. (C.1) is clearly suppressed for $b_{T}P^{z}\gg 1$. The second line in Eq. (C.1) also vanishes in this limit, because the incomplete Gamma function behaves as

[TABLE]

and the logarithmic enhancement for small arguments is compensated by the prefactor. The third line of Eq. (C.1) yields a power-suppressed contribution when convoluted with the TMDPDF, see appendix C.4.

C.2 Sail diagrams

The two sail diagrams are shown in Fig. 12. The Feynman rule for the Wilson line following from Eq. (153) is given by

[TABLE]

where $k$ is the momentum flowing into the Wilson line and $\gamma$ the Wilson line path, given in Eq. (154).

The two diagrams prior to Fourier transforming are given by

[TABLE]

After averaging over spin, they can be combined into

[TABLE]

For both $\lambda=0$ and $\lambda=z$, the numerator can be simplified to

[TABLE]

where the last step employs that the $g^{\mu 0}$ piece vanishes when contracted with the spatial path $\gamma^{\prime}(s)_{\mu}$. Introducing Feynman parameters and dropping the terms linear in $q^{0}$ gives

[TABLE]

Next, we use the relation

[TABLE]

to replace the $q^{\mu}$ piece in the numerator,

[TABLE]

Since the first line in Eq. (C.2) is much simpler than the second line, we have split $\tilde{q}_{n}^{(b)}$ accordingly. Note that the second line in Eq. (C.2) only involves $\gamma^{\prime}(s)^{z}$. All contributions from the transverse gauge link are thus fully captured in the first line, which becomes

[TABLE]

Taking the Fourier transform of this and pulling it inside a plus distribution yields

[TABLE]

The second contribution to Eq. (C.2) is more involved. We first simplify the exponentials by letting $y\to 1-y$ and then $q\to q+yP$,

[TABLE]

To proceed, we need the explicit parameterizations for the Wilson line path $\gamma$, Eq. (154). Evaluating all line integrals and taking the Fourier transform w.r.t. $b^{z}$, Eq. (181) yields

[TABLE]

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,

[TABLE]

For $x\neq 1$, this can be expressed using our master integrals,

[TABLE]

The $\delta(1-x)$ term is fixed by integrating $\tilde{q}^{(b,2)}_{n}(x,\vec{b}_{T},P^{z})$ over $x$,

[TABLE]

Here, the finite $L$ terms are crucial to regulate the integral at $q^{z}=0$, which makes the integral IR finite. Note that due to the factor $1/q^{z}$, the above integral is UV finite, so we can evaluate the integral in $d=4$ dimensions, which after some manipulations yields

[TABLE]

We have not been able to obtain a closed form for this integral. However, it is clearly convergent for $v\to 0$ and $v\to\infty$, so we can take the $LP^{z}\to\infty$ limit first. To cleanly extract the logarithmic dependence on $P^{z}$, take

[TABLE]

The full expression for the sail diagram, Eq. (C.2), is obtained by adding Eqs. (180) and (183),

[TABLE]

Using the approximate result Eq. (C.2), but keeping all master integrals exact, we obtain

[TABLE]

where we have singled out the terms with unphysical support $x\in[-\infty,\infty]$,

[TABLE]

C.3 Wilson line self energy

The general expression for the Wilson line self energy, Fig. 11(c), in position space is given by

[TABLE]

where a symmetry factor $1/2$ is included and $\gamma$ is the path of the Wilson line. Since $\gamma$ can be split into three straight lines, see Fig. 5, there are four distinct contributions in Feynman gauge, shown in Fig. 13.

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

[TABLE]

Diagram (c)

depends on $b^{z}$, so care has to be taken with the Fourier transform. The position space result is

[TABLE]

The Fourier transform can be obtained using Izubuchi:2018srq

[TABLE]

as

[TABLE]

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.3) in $b^{z}\ll L$ prior to Fourier transforming, one would instead obtain

[TABLE]

so dropping the phase in Eq. (196) is equivalent to performing a small-distance expansion $b^{z}\ll L$ in position space.

Diagram (d)

gives rise to a factor $2$ to lift the symmetry factor. Due to the finite separation $\vec{b}_{T}$, it is UV finite, so we can let $\epsilon\to 0$. In position space, we obtain

[TABLE]

Here, it is quite difficult to take the Fourier transform while keeping the exact $b^{z}$ dependence. It is easier to use the first line of Eq. (C.3) and leave the $k^{z}$ integration until the end,

[TABLE]

While the second line is manifestly finite as $k^{z}\to 0$, the last line has an apparent singularity, which is as usual treated as a plus distribution,

[TABLE]

For comparison, the $b^{z}\ll L$ limit of Eq. (C.3) is given by

[TABLE]

This still has a nontrivial $b^{z}$ dependence, so in order to obtain the $\delta(1-x)$ term in Eq. (C.3) we further need to expand in $b^{z}\ll b_{T}$, giving

[TABLE]

Hence dropping the phase term in Eq. (C.3) corresponds to a small-distance expansion $b^{z}\ll L$ and $b^{z}\ll b_{T}$.

Combined result.

Combining Eqs. (192), (196) and (C.3), we obtain the full exact tadpole diagram as

[TABLE]

where the terms suppressed for $b^{z}\ll b_{T}\ll L$ are

[TABLE]

C.4 Power-suppressed contributions to the matching kernel

The unphysical contributions Eqs. (C.1), (C.2) and (C.3) contain fastly oscillating phases $\sim e^{-\mathrm{i}P^{z}L(1-x)}$, and one may thus expect that these give vanishing contributions as $LP^{z}\to\infty$. However, these phases vanish as $x\to 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}\to\infty$.

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.1), then the unphysical part of the matching is given by

[TABLE]

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.4) to infinity by implicitly assuming that $f^{\mathrm{TMD}}(y)=0$ for $|y|\geq 1$. Changing the integration variable in Eq. (C.4) from $y$ to $x/y$, we obtain

[TABLE]

Next note that $\Delta\tilde{q}(y)$ is entirely given as a plus distribution, so we can rewrite Eq. (206) as

[TABLE]

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 $\Delta\tilde{q}$ is given by an exponential phase factor, a regular function $r(y)$ and potentially a divergence as $y\to 1$,

[TABLE]

For this case, Eq. (C.4) becomes

[TABLE]

The function $h(y)$ is by construction regular as $z\to 1$ due to the subtraction term.

Without giving an exact proof for all functions $h(y)$ appearing at one loop, we can give strong arguments that Eq. (C.4) vanishes as $P^{z}L\to\infty$ for all relevant functions $h(y)$. The strongest criterion to prove this behavior is for functions $h$ that satisfy $\int\mathrm{d}y\,|h(y)|<\infty$, in which case the Riemann-Lebesgue lemma applies. Typically, this is too restrictive, as for example $h(x)=\ln|x|$ is not integrable, but its Fourier transform is known to vanish as $1/(P^{z}L)$. A less restrictive and more intuitive argument is that as long as $h(y)$ is sufficiently smooth such that $P^{z}L\gg|h^{\prime}(y)|$, then the rapid oscillation of the Fourier kernel suppresses the integral. This holds for all $(1-y)P^{z}L\gtrsim 1$, and thus the contribution of the integral not suppressed by the exponential phase must be of order $1/(P^{z}L)$.

All other functions appearing in Eqs. (C.1), (C.2) and (C.3) have an explicit suppression in $b_{T}P^{z}\to\infty$. While they may contain divergences in $1/(1-y)$, they can be regulated the same way as the above, after which the result is clearly suppressed by positive powers of $1/(b_{T}P^{z})$ relative to the physical terms in the matching calculation.

Since all terms in Eqs. (C.1), (C.2) and (C.3) give power-suppressed contributions to the matching, they can be dropped in the final result for the unsubtracted beam function. Upon adding only the physical contributions from Eqs. (C.1), (C.2) and (203), one then obtains Eq. (53).