# Factorization connecting continuum & lattice TMDs

Markus A. Ebert  Stella T. Schindler  Iain W. Stewart  and Yong Zhao

###### Abstract

Transverse-momentum-dependent parton distribution functions (TMDs) can be
studied from first principles by a perturbative matching onto lattice-
calculable quantities: so-called lattice TMDs, which are a class of equal-time
correlators that includes quasi-TMDs and TMDs in the Lorentz-invariant
approach. We introduce a general correlator that includes as special cases
these two Lattice TMDs and continuum TMDs, like the Collins scheme. Then, to
facilitate the derivation of a factorization relation between lattice and
continuum TMDs, we construct a new scheme, the Large Rapidity (LR) scheme,
intermediate between the Collins and quasi-TMDs. The LR and Collins schemes
differ only by an order of limits, and can be matched onto one another by a
multiplicative kernel. We show that this same matching also holds between
quasi and Collins TMDs, which enables us to prove a factorization relation
between these quantities to all orders in $\alpha_{s}$. Our results imply that
there is no mixing between various quark flavors or gluons when matching
Collins and quasi TMDs, making the lattice calculation of individual flavors
and gluon TMDs easier than anticipated. We cross-check these results
explicitly at one loop and discuss implications for other physical-to-lattice
scheme factorizations.

##  1 Introduction

Nuclear and particle physicists have striven to ascertain the full three-
dimensional structure of hadrons for decades, through a combination of
experimental measurements and their theoretical description via parton
distribution functions (PDFs). The study of PDFs in the longitudinal momentum
direction has reached a high level of maturity within many frameworks, such as
global fitting efforts [1, 2, 3, 4, 5, 6, 7, 8], lattice QCD calculations
using moments [9, 10, 11], quasi-distributions in large-momentum effective
theory (LaMET) [12, 13, 14], and various other techniques [15, 16, 17, 18, 19,
20, 21, 22]. A major thrust of current research is to extend this progress to
transverse-momentum-dependent PDFs (TMDs). In principle, TMDs can be extracted
from experiments using global fits [23, 24], but due to limited data these
extractions have not yet yielded the same precision as PDFs. However, ongoing
and planned experiments such as the Electron-Ion Collider [25, 26] will have
enhanced capabilities to probe TMDs, and it is important to complement them
with theoretical insight.

TMDs are intrinsically nonperturbative at small transverse momenta and so
lattice QCD provides the only practical method for their first-principles
calculation. Unfortunately, TMDs are defined along nonlocal lightlike or
close-to-lightlike Wilson line paths, which depend explicitly on a real-valued
time variable and thus lie outside the reach of current lattice techniques due
to a sign problem, the speculated non-deterministic non-polynomial hard issue
of numerically averaging over complex weights in a Euclidean path integral. To
circumvent this obstacle, one can construct equal-time correlators that are
calculable on the lattice, which we collectively dub “Lattice TMDs”. We can
then extract information from lattice calculations by relating a lattice TMD
to a physical continuum TMD (that appear in cross sections) through a so-
called factorization formula, which forms the focus of this paper.

There exist several approaches for defining continuum and lattice TMDs. The
first lattice TMD studies used the Lorentz-invariant approach of refs. [27,
28, 29, 30, 31, 32], which we refer to as the Musch-Hägler-Engelhardt-Negele-
Schäfer (MHENS) scheme. This approach defines TMDs which are calculated with
lattice QCD using spacelike Wilson line paths, and then uses Lorentz
invariance to relate them to the path considered for the continuum TMD
definitions, which is spacelike but close to the light-cone. This method has
so far primarily been used to calculate ratios of physical TMD moments. Later
on, LaMET motivated the study of quasi-TMDs [33, 34, 35, 36, 37, 38, 39, 40,
41, 42], which are Euclidean distributions defined using boosted hadron matrix
elements of somewhat different spacelike Wilson line paths. Since the quasi-
TMD obeys the Collins-Soper evolution [33], one can extract the
nonperturbative Collins-Soper (CS) evolution of TMDs using ratios of quasi-
TMDs at different hadron momenta [35], and first lattice results were obtained
in refs. [43, 44, 45, 46, 47, 48]. Methods to calculate the so-called soft
function have also been proposed [38] and implemented on the lattice [46, 47].

To relate lattice TMDs and physical continuum TMDs one should derive a
factorization formula that demonstrates that these TMDs agree in the infrared,
and perhaps differ in the ultraviolet by perturbative matching coefficients
[33, 34, 36, 39, 41, 49]. In this work, we set up a unified notational
framework for lattice and continuum TMDs and derive the factorization formula
between the physical Collins TMD and quasi-TMD to all orders in $\alpha_{s}$,
for both quark and gluon TMDs. Up until now, no lattice-to-continuum
factorization formlula has been proven, but a matching between the quasi-TMD
and Collins-TMD for the non-singlet quark case has been verified at one-loop
order.111 In ref. [39] arguments for the factorization based on an analysis of
leading regions were suggested. The presentation of such an analysis will be
an important complement to the proof given here. Ref. [49] analyzes
factorization for lattice-friendly correlators. The transverse distribution
that they study there is not the same as the quasi-TMD in refs. [33, 34, 36,
39, 41] and this work.  In our analysis here we also fully account for lattice
renormalization, soft function subtractions, and finite Wilson line lengths;
items that are important to account for in analyses for lattice-QCD-friendly
TMDs.

###  1.1 Statement of factorization

Factorization formulae that relate the non-singlet quark quasi-TMD and
physical TMD have been proposed in refs. [35, 36, 38, 39]. The objects we are
interested in here include both a naive quasi-TMD $\tilde{f}_{{\rm ns}}^{\rm
naive}$ that can be directly calculated with lattice QCD, but whose infrared
structure differs from the physical TMD, and a proper quasi-TMD
$\tilde{f}_{{\rm ns}}$ whose infrared structure agrees. The factorization
theorem for both of these quasi-TMDs can be expressed as

$\displaystyle\tilde{f}^{\text{(naive)}}_{\rm
ns}(x,\vec{b}_{T},\mu,x\tilde{P}^{z})$ $\displaystyle=C_{\rm
ns}(x\tilde{P}^{z},\mu)g_{S}^{q}(b_{T},\mu)\exp\biggl{[}{1\over
2}\gamma_{\zeta}^{q}(b_{T},\mu)\ln{(2x\tilde{P}^{z})^{2}\over\zeta}\biggr{]}f_{{\rm
ns}}(x,\vec{b}_{T},\mu,\zeta)$
$\displaystyle\qquad\times\bigg{\\{}1+\mathcal{O}\biggl{[}{1\over(x\tilde{P}^{z}b_{T})^{2}},{\Lambda_{\rm
QCD}^{2}\over(x\tilde{P}^{z})^{2}}\biggr{]}\bigg{\\}}\,.$ (1.1)

Here ${\rm ns}=u-d$ for non-singlet flavor, and $C_{\rm ns}$ is a perturbative
matching coefficient which does not depend on spin [40, 49]. In eq. (1.1) $x$
is the fraction of the hadron’s longitudinal momentum carried by the parton,
$\vec{b}_{T}$ with $b_{T}=|\vec{b}_{T}|$ is the Fourier conjugate of parton
transverse momentum, and $\mu$ is the ${\overline{\mathrm{MS}}}$
renormalization scale. The quasi-TMD depends on the hadron momentum
$\tilde{P}^{z}$, while the TMD depends on the CS scale $\zeta$; also,
$\gamma_{\zeta}^{q}$ is the anomalous dimension for $\zeta$-evolution, which
is often referred to as the CS kernel [50, 51, 52]. Since we are interested in
nonperturbative $b_{T}\sim\Lambda_{\mathrm{QCD}}^{-1}$, for the proper quasi-
TMD $\tilde{f}_{\rm ns}$ we must have $g_{S}^{q}=1$, while for the naive
quasi-TMD $f_{\rm ns}^{\rm naive}$ the function $g_{S}^{q}(b_{T},\mu)$ is
nonperturbative and arises from using a naive quasi-soft function
$\tilde{S}^{R}_{\rm naive}$, but can be calculated from the reduced soft
function with the methods proposed in ref. [39]. The definition of the proper
quasi-TMD we use here is not the same as earlier literature, see section 2.3.1
for more details. This factorization is valid at large but finite
$\tilde{P}^{z}$, with power corrections suppressed by the parton momentum
$x\tilde{P}^{z}$ as shown in eq. (1.1).

In this work, we prove eq. (1.1) at all orders in $\alpha_{s}$ and generalize
it to the quasi-TMDs of all partons $i$, including light quark flavors
$(u,d,s)$ and gluons. Specifically, for the TMDs for hadron $h$ we find

$\displaystyle\tilde{f}_{i/h}^{[s]}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z})$
$\displaystyle=C_{i}(x\tilde{P}^{z},\mu)\exp\biggl{[}\frac{1}{2}\gamma_{\zeta}^{i}(\mu,b_{T})\ln\frac{\tilde{\zeta}}{\zeta}\biggr{]}f_{i/h}^{[s]}(x,\vec{b}_{T},\mu,\zeta)$
$\displaystyle\quad\times\bigg{\\{}1+\mathcal{O}\biggl{[}{1\over(x\tilde{P}^{z}b_{T})^{2}},{\Lambda_{\rm
QCD}^{2}\over(x\tilde{P}^{z})^{2}}\biggr{]}\bigg{\\}}\,,$ (1.2)

where $[s]$ denotes the choice of spin-polarization for the hadronic state and
operator, and $C_{i}$ and $\gamma_{\zeta}^{i}$ depend only on the choice of
fundamental or adjoint color representations, so $C_{i}=C_{q}$ and
$\gamma_{\zeta}^{i}=\gamma_{\zeta}^{q}$ for all quarks, and $C_{i}=C_{g}$ and
$\gamma^{i}_{\zeta}=\gamma^{g}_{\zeta}$ for gluons. The quasi-TMD in eq. (1.1)
differs from eq. (1.1) by the soft factor,

$\displaystyle\tilde{f}_{i/h}^{[s]}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z})$
$\displaystyle=\tilde{f}_{i/h}^{[s]\rm
naive}(x,\vec{b}_{T},\mu,x\tilde{P}^{z})\sqrt{\frac{\tilde{S}^{R}_{\rm
naive}(b_{T},\mu)}{S_{C}^{R}(b_{T},\mu,2y_{n},2y_{B})}}\,,$ (1.3)

where $S_{C}^{R}(b_{T},\mu,2y_{n},2y_{B})$ is the soft function in the Collins
scheme [53], $R=q,g$ for quarks in the fundamental or gluons in the adjoint
representation, and $y_{n}$ and $y_{B}$ are two different rapidities (see
section 2.2.1). Notably, the quasi-TMD depends on a new variable
$\tilde{\zeta}=x^{2}m_{h}^{2}e^{2(y_{\tilde{P}}+y_{B}-y_{n})}$ with hadron
mass $m_{h}$ and rapidity $y_{\tilde{P}}$, which is equivalent to a CS scale.
Comparing Eqs. (1.1) – (1.3), we also show that $g_{S}^{q}$ can be calculated
as a ratio of soft functions, in agreement with the definition of the reduced
soft function in ref. [39].

Two key results of this work are that eq. (1.1) holds for both quark and gluon
quasi-TMDs, and there is no mixing between the quark and gluon channels or
quarks of different flavors. This differs from the case of the factorizations
of quark and gluon longitudinal quasi-PDFs [54]. The matching coefficient
$C_{i}$ takes different values for the two different cases $i=q,\,g$, but is
independent of quark flavor and spin.

###  1.2 Strategy for proving factorization

We begin by developing a unified notational framework that is applicable to
both lattice and continuum TMDs. This framework brings to the forefront the
common underlying structural features of different TMD schemes, allowing one
to more easily construct factorization formulae relating them to one another.
We can write every lattice and continuum TMD as the product of a proton matrix
element (beam function), vacuum matrix element (soft factor), and appropriate
renormalization factors. In the cases studied in the literature, we show that
the beam function can be expressed as the Fourier transform of a generic,
common correlator $\Omega$. Each scheme manifests as a special case of the
correlator, encoded by a choice of $\Omega$’s arguments and of the order in
which we take parameter limits needed for proper TMD regularization and
renormalization. More details are provided later, and we summarize the choices
needed for various schemes in Tables 1 and 2.

Having a unified notational framework is useful for constructing relationships
between continuum TMD schemes that can be connected to physical observables,
with schemes that can be computed on the lattice. From the structure of the
correlator $\Omega$, we observe that the Collins and quasi-TMDs are closely
related. To relate these schemes we begin by constructing a new TMD scheme
that is intermediate between the quasi and Collins TMDs: the Large-Rapidity
(LR) scheme. The LR scheme uses the same ingredients as the Collins scheme,
but performs UV renormalization at large but finite Wilson line rapidities.
Using Lorentz invariance, we show that the quasi- and LR scheme TMDs are
equivalent, up to terms suppressed at large proton rapidity. We then show that
we can relate the LR and Collins TMDs simply by a perturbative matching
kernel, which is flavor-diagonal and spin-independent for both quarks and
gluons. Combining the quasi-to-LR and LR-to-Collins relations leads to eq.
(1.1). We summarize relationships between a number of common TMD schemes in
figure -119.

This proof is a beautiful application of the fundamental principle underlying
LaMET [12, 13, 14]: partons encode the internal degrees of freedom of a highly
energetic hadron, with the hadron momentum limit $P^{z}\to\infty$ taken prior
to UV renormalization. In contrast, on the lattice one must carry out UV
renormalization prior to the large-momentum limit. This different order of
limits in an asymptotically free theory like QCD induces a nontrivial matching
kernel between a parton observable and its corresponding lattice construction.

Peculiarly, in TMDs this noncommutativity of orders of limits appears
naturally, since at intermediate steps of the calculation one must regulate
so-called rapidity divergences [53, 50, 55, 56, 57, 58, 59, 60, 61, 62, 63,
64], usually by choosing to deviate from lightcone kinematics or by
introducing an additional regulator on the lightcone. These divergences cancel
when combining the hadronic and soft matrix elements, allowing one to take the
lightcone limit (or infinite rapidity-regulator limit), but also forcing one
to decide whether to take this limit before or after UV renormalization. Many
schemes for constructing TMDs take the lightcone limit first, e.g. the Collins
scheme [53]; others, such as the Ji-Ma-Yuan (JMY) scheme, do the opposite
[56]. Once again, exchanging these orders of limits is a UV effect, and thus
induces a nontrivial matching relation between Collins and JMY TMDs.

This paper is structured as follows: section 2 introduces our new notational
framework for defining TMDs. We then provide an overview of the definitions of
common physical and lattice TMD schemes. Section 3 presents the proof of the
factorization statement in eq. (1.1). We then discuss the physical
implications of this proof, in particular the lack of flavor mixing in the
matching and utility for calculating ratios of TMDs. We also confirm the
momentum evolution RGE equation for the hard matching coefficient, and give a
complete solution for it with next-to-leading-logarithmic resummation. From
our proof, we also gain intuition for factorization relations between other
lattice and continuum TMD schemes, in particular for the MHENS and Collins
schemes which differ due to the presence of Wilson line cusp angles in the
former that are not in the latter, and by the need for different soft factors
in these two TMDs. We make concluding remarks and outline future directions in
section 4. In the appendices, we verify our quasi-to-Collins factorization
results analytically at one-loop order, present a one-loop comparison of the
continuum JMY and Collins TMDs, and also discuss Wilson line self-energies.

Figure -119: An overview of schemes and their relationships, including the LR
scheme introduced in this work. See section 3 for details.

##  2 Definition of TMDs

Let us consider the process of hard hadron-hadron scattering:

$h_{1}(P_{1})+h_{2}(P_{2})\to L(q)+X\,,$ (2.1)

where $h_{1,2}$ are colliding hadrons with momenta $P_{1,2}$, $L$ is a
detected color-singlet final state with momentum $q^{\mu}$ (such as
$L=Z/\gamma^{*}$, which decays to leptons), and $X$ denotes additional final-
state particles. Let $L$ have invariant mass $Q$, rapidity $Y$, and transverse
momentum $\vec{q}_{T}$ with $q_{T}=|\vec{q}_{T}|$. For $q_{T}\ll Q$, we can
factorize the cross-section of eq. (2.1) as [50, 51, 52, 53]

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}Q^{2}\mathrm{d}Y\mathrm{d}^{2}\vec{q}_{T}}$
$\displaystyle=\sigma_{0}\sum_{i,j}H_{ij}(Q,\dots)\int\\!\frac{\mathrm{d}^{2}\vec{b}_{T}}{(2\pi)^{2}}e^{\mathrm{i}\vec{q}_{T}\cdot\vec{b}_{T}}f_{i/h_{1}}\bigl{(}x_{1},\vec{b}_{T},\dots\bigr{)}f_{j/h_{2}}\bigl{(}x_{2},\vec{b}_{T},\dots\bigr{)}\,.$
(2.2)

Here, $\sigma_{0}$ is the Born cross-section; the sum runs over all parton
flavors $i,j$ contributing to the Born process $ij\to L$; $H_{ij}$ is the hard
function, which encodes virtual corrections to the Born process; and $f_{i/h}$
are the TMDs, functions which describe the dynamics of partonic quarks and
gluons inside the parent hadron $h$. A struck hadron carries a fraction
$x_{1,2}=Qe^{\pm Y}/E_{\rm cm}$ of its parent hadron’s longitudinal momentum,
with $E_{\rm cm}=\sqrt{(P_{1}+P_{2})^{2}}$ the center-of-mass energy of the
incoming hadrons. The ellipses in eq. (2.1) denote additional parameters
related to UV and rapidity renormalization, whose precise forms are scheme
dependent. Note that we suppress indices related to spin-dependent processes
and contributions.

The literature is rife with schemes for defining TMDs, each of which has
different strengths for different types of calculations. This section reviews
schemes relevant for lattice studies; in particular, we only discuss schemes
based on off-lightcone Wilson lines. Schemes with intrinsically lightlike
Wilson lines [65, 58, 62, 60, 63, 64] are not accessible on a Euclidean
lattice, but many are equivalent to the Collins scheme once limits needed to
obtain TMD PDFs are taken; see refs. [36, 66] for an overview. Because each
scheme in the literature employs its own conventions and notation, in section
2.1 we begin by introducing new unified TMD Lorentz-invariant correlators for
which all schemes follow as special cases. Then, we provide definitions of
physical and lattice schemes in sections 2.2 and 2.3, respectively.

###  2.1 Unified TMD notation

A TMD $f_{i/h}$ generally contains two pieces: a hadronic matrix element
(called the beam function or unsubtracted TMD), which encodes partonic
radiation associated with the initial hadrons; as well as a vacuum matrix
element (the soft function). These matrix elements involve open and closed
staple-shaped Wilson lines, for which we develop a generic notation. Let us
first define a Wilson line along a path $\gamma$ in color representation $R$
as:

$\displaystyle
W^{R}[\gamma]=P\exp\biggl{[}\mathrm{i}g\int_{\gamma}\mathrm{d}x^{\mu}\mathcal{A}_{\mu}^{a}(x)\,T_{R}^{a}\biggr{]}\,,$
(2.3)

where $R=q$ in the fundamental and $R=g$ in the adjoint representation. It is
useful to define a general class of Wilson lines using the three-sided staple
shape shown in figure -118,

$\displaystyle W_{\sqsupset}^{R}(b,\eta
v,\delta)=W^{R}\left[\frac{b}{2}~{}\to~{}\frac{b}{2}+\eta
v-\frac{\delta}{2}~{}\to~{}-\frac{b}{2}+\eta
v+\frac{\delta}{2}~{}\to-\frac{b}{2}\right]\,.$ (2.4)

The length of the staple is relevant for its renormalization properties; here
we have

$\displaystyle L_{\rm staple}$ $\displaystyle=|\eta v-\delta/2|+|\eta
v+\delta/2|+|b-\delta|\,,$ (2.5)

where the length of a four-vector is given by $|X|=\sqrt{|X^{2}|}$. At the red
points in figure -118, the staple has cusp angles $\gamma_{\pm}$, which can be
computed from

$\displaystyle\cosh\gamma_{\pm}$ $\displaystyle=\frac{(\eta
v\pm\delta/2)\cdot(b-\delta)}{|\eta v\pm\delta/2||b-\delta|}\,,$ (2.6)

where for space-like separations
$\gamma_{\pm}\in[-\mathrm{i}\pi,\mathrm{i}\pi]$.222Note that we develop our
generic TMD framework with a _three_ -sided staple. Adding more than three
sides will induce extra Wilson line cusps, that create additional
complications for renormalization.

Figure -118: Generic staple-shaped Wilson line defined in eq. (2.4). Black
double-lines extend along $\eta v$, and the blue segment along $b-\delta$
closes the staple. For certain choices of $\delta$, red points can be cusps.
Edges may extend along the conjugate direction $P$, which is not shown.

Generic quark and gluon beam function correlators take the form

$\displaystyle\Omega_{q_{i}/h}^{[\Gamma]}(b,P,\epsilon,\eta v,\delta)$
$\displaystyle=\Bigl{\langle}h(P)\Big{|}\bar{q}_{i}\Bigl{(}\frac{b}{2}\Bigr{)}\frac{\Gamma}{2}W^{F}_{\sqsupset}(b,\eta
v,\delta)q_{i}\Bigl{(}-\frac{b}{2}\Bigr{)}\Big{|}h(P)\Bigr{\rangle}\,,$
$\displaystyle\Omega_{g/h}^{\mu\nu\rho\sigma}(b,P,\epsilon,\eta v,\delta)$
$\displaystyle=\Bigl{\langle}h(P)\Big{|}G^{\mu\nu}\Bigl{(}\frac{b}{2}\Bigr{)}W^{A}_{\sqsupset}(b,\eta
v,\delta)G^{\rho\sigma}\Bigl{(}-\frac{b}{2}\Bigr{)}\Big{|}h(P)\Bigr{\rangle}\,.$
(2.7)

In eq. (2.1), $q_{i}(x)$ is a quark field of flavor $i$, and $G^{\mu\nu}(x)$
is the gluon field strength tensor. The quark and gluon fields are spatially
separated by $b$, which is Fourier-conjugate to the momentum of the struck
parton. In the quark correlator, $\Gamma$ denotes a generic Dirac structure,
while for the gluon correlator $\mu,\,\nu,\,\rho,$ and $\sigma$ are Lorentz
indices. See refs. [67, 40] for decompositions of different choices of
$\Gamma$ into independent spin structures for quark TMDs, and refs. [68, 69]
for the decomposition for gluon TMDs. In both cases, $h$ denotes the struck
hadron with momentum $P$, $\epsilon$ is the UV regulator, and $\eta v$ and
$\delta$ characterize the longitudinal and transverse segments of the Wilson
line, which we illustrate in figure -118.

We define the generic soft vacuum matrix element as

$\displaystyle S^{R}(b,\epsilon,\eta v,\bar{\eta}\bar{v})$
$\displaystyle=\frac{1}{d_{R}}\Big{\langle}0\Big{|}\mathrm{Tr}\Bigl{[}S^{R}_{\leavevmode\hbox
to4.69pt{\vbox to7.65pt{\pgfpicture\makeatletter\hbox{\hskip
0.2134pt\lower-0.327pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {}{{}}{}{{}}{}{{}}{}
{}{}{{}}{} {}{}{{}}{}{{}}{} {}{}{{}}{}{{}}{}{{}}{} {}{}{{}}{} {}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.42679pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{4.26773pt}{2.84544pt}\pgfsys@lineto{0.0pt}{5.69089pt}\pgfsys@lineto{0.0pt}{7.1136pt}\pgfsys@lineto{4.26773pt}{4.26816pt}\pgfsys@lineto{0.0pt}{1.42271pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.34143pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{-0.15628pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}(b,\eta
v,\bar{\eta}\bar{v})\Bigr{]}\Big{|}0\Bigr{\rangle}\,,$ (2.8)

where the trace is over color. The color averaging factor $d_{R}$ takes values
$d_{q}=N_{c}$ and $d_{g}=N_{c}^{2}-1$. The soft Wilson line is given by

$\displaystyle S^{R}_{\leavevmode\hbox to4.69pt{\vbox
to7.65pt{\pgfpicture\makeatletter\hbox{\hskip 0.2134pt\lower-0.327pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {}{{}}{}{{}}{}{{}}{}
{}{}{{}}{} {}{}{{}}{}{{}}{} {}{}{{}}{}{{}}{}{{}}{} {}{}{{}}{} {}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.42679pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{4.26773pt}{2.84544pt}\pgfsys@lineto{0.0pt}{5.69089pt}\pgfsys@lineto{0.0pt}{7.1136pt}\pgfsys@lineto{4.26773pt}{4.26816pt}\pgfsys@lineto{0.0pt}{1.42271pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.34143pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{-0.15628pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}(b,\eta
v,\bar{\eta}\bar{v})$
$\displaystyle=W^{R}\bigg{[}\frac{b}{2}~{}\to~{}\frac{b}{2}+\bar{\eta}\bar{v}~{}~{}~{}\to~{}-\frac{b}{2}+\bar{\eta}\bar{v}~{}\to~{}-\frac{b}{2}$
$\displaystyle\qquad\qquad~{}\to~{}-\frac{b}{2}+\eta
v~{}\to~{}\frac{b}{2}+\eta v~{}~{}~{}\to~{}\frac{b}{2}\bigg{]}\,,$ (2.9)

as shown in figure -117. $S_{\leavevmode\hbox to4.69pt{\vbox
to7.65pt{\pgfpicture\makeatletter\hbox{\hskip 0.2134pt\lower-0.327pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {}{{}}{}{{}}{}{{}}{}
{}{}{{}}{} {}{}{{}}{}{{}}{} {}{}{{}}{}{{}}{}{{}}{} {}{}{{}}{} {}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.42679pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{4.26773pt}{2.84544pt}\pgfsys@lineto{0.0pt}{5.69089pt}\pgfsys@lineto{0.0pt}{7.1136pt}\pgfsys@lineto{4.26773pt}{4.26816pt}\pgfsys@lineto{0.0pt}{1.42271pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.34143pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{-0.15628pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}$ consists of two
beam function staples glued together at the points $\pm b/2$; the long sides
of the staples run along the $\bar{\eta}\bar{v}$ and $\eta v$ directions. The
dependence on two conjugate directions arises from the appearance of two TMDs
in the physical cross section in eq. (2.2). The length of the soft function
path is $L_{\leavevmode\hbox to4.69pt{\vbox
to7.65pt{\pgfpicture\makeatletter\hbox{\hskip 0.2134pt\lower-0.327pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {}{{}}{}{{}}{}{{}}{}
{}{}{{}}{} {}{}{{}}{}{{}}{} {}{}{{}}{}{{}}{}{{}}{} {}{}{{}}{} {}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.42679pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{4.26773pt}{2.84544pt}\pgfsys@lineto{0.0pt}{5.69089pt}\pgfsys@lineto{0.0pt}{7.1136pt}\pgfsys@lineto{4.26773pt}{4.26816pt}\pgfsys@lineto{0.0pt}{1.42271pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.34143pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{-0.15628pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}=2|\bar{\eta}\bar{v}|+2|\eta
v|+2|b|$.

Figure -117: Wilson line structure of the soft function, eq. (2.8) for
$\eta,\bar{\eta}<0$. Figure adapted from ref. [63].

We define the transverse direction with respect to the plane spanned by $P$
and $v$, taking $P_{\perp}=v_{\perp}=0$. Formally, this can be expressed as
$b_{\perp}^{\mu}=g_{\perp}^{\mu\nu}b_{\nu}$ with

$\displaystyle g_{\perp}^{\mu\nu}$
$\displaystyle=g^{\mu\nu}-\frac{1}{1+\hat{\zeta}^{2}}\biggl{[}\frac{v^{\mu}v^{\nu}}{v^{2}}+\frac{P^{\mu}P^{\nu}}{P^{2}}+\frac{\hat{\zeta}^{2}}{P\cdot
v}\bigl{(}P^{\mu}v^{\nu}+v^{\mu}P^{\nu}\bigr{)}\biggr{]}\,,\quad\hat{\zeta}=\frac{v\cdot
P}{\sqrt{|v^{2}|P^{2}}}\,.$ (2.10)

We always take $v$ and $P$ to span the same plane as $v$ and $\bar{v}$. It
follows that $v_{\perp}=\bar{v}_{\perp}=0$.

Our unified notation facilitates the comparison of different TMD schemes,
particularly when we examine their Lorentz invariants. In the most generic
case, the beam function correlator in eq. (2.1) is specified by four
independent vectors: $b^{\mu},P^{\mu},\eta v^{\mu}$, and $\delta^{\mu}$. From
these vectors we can construct ten independent Lorentz invariants, which we
choose to be

$\displaystyle P^{2}\,,\quad b^{2}\,,\quad\eta^{2}v^{2}\,,\quad P\cdot
b\,,\quad\frac{P\cdot(\eta v)}{\sqrt{P^{2}|(\eta
v)^{2}|}}\,,\quad\frac{b\cdot(\eta v)}{\sqrt{|b^{2}(\eta v)^{2}|}}\,,\quad$
$\displaystyle\frac{\delta^{2}}{b^{2}}\,,\quad\frac{b\cdot\delta}{b^{2}}\,,\quad\frac{P\cdot\delta}{P\cdot
b}\,,\quad\frac{\delta\cdot(\eta v)}{b\cdot(\eta v)}\,.$ (2.11)

None of the TMD schemes we study in this paper contains a vector
$\delta^{\mu}$ that is linearly independent of $b^{\mu},P^{\mu},$ and $\eta
v^{\mu}$; thus, these schemes have six independent Lorentz invariants.
However, the quasi and MHENS TMDs do not follow from the same correlator
defined with six invariants, since they fix $\delta^{\mu}$ in different ways,
as we will see below. Hence, even if the first six invariants in eq. (2.1) are
fixed to be the same, the two approaches have different values for the last
four invariants, and thus the quasi- and MHENS TMDs belong to distinct
schemes.

###  2.2 Continuum TMD schemes

In this section, we provide an overview of physical TMD schemes, which are
defined on a continuous spacetime and have infinitely long Wilson lines, with
$|\eta|=|\bar{\eta}|=\infty$.

##### Lightcone coordinate conventions.

It is convenient to work in a frame where the hadron momenta $P_{1,2}$ in eq.
(2.1) are close to the lightlike unit vectors

$\displaystyle n_{a}^{\mu}=\frac{1}{\sqrt{2}}(1,0,0,1)\,,\qquad
n_{b}^{\mu}=\frac{1}{\sqrt{2}}(1,0,0,-1)\,,$ (2.12)

which obey $n_{a}^{2}=n_{b}^{2}=0$ and $n_{a}\cdot n_{b}=1$. We define the
lightcone decomposition of an arbitrary four-vector $p^{\mu}$ as

$\displaystyle
p^{\mu}=(p^{+},p^{-},p_{\perp})=p^{+}n_{a}^{\mu}+p^{-}n_{b}^{\mu}+p_{\perp}^{\mu}\,,$
(2.13)

where $p^{\pm}=(p^{0}\pm p^{z})/\sqrt{2}$ and
$p_{\perp}^{\mu}=(0,p_{x},p_{y},0)=(0,\vec{p}_{T},0)$. Here $p_{\perp}^{\mu}$
is a Minkowski vector, and $\vec{p}_{T}$ is the corresponding transverse
Euclidean vector with magnitude
$p_{T}\equiv(\vec{p}_{T}^{2})^{1/2}=(-p_{\perp}^{2})^{1/2}$. In lightcone
coordinates, the incoming hadrons in eq. (2.1) have momenta

$\displaystyle
P_{1}^{\mu}=P_{1}^{+}\bigl{(}1,e^{-2y_{1}},0_{\perp}\bigr{)}\,,\qquad
P_{2}^{\mu}=P_{2}^{-}\bigl{(}e^{+2y_{2}},1,0_{\perp}\bigr{)}\,,$ (2.14)

where $y_{1,2}$ are the hadron rapidities.

####  2.2.1 Collins scheme

In the Collins TMD scheme [53], the factorization formula in eq. (2.1) takes
the form

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}Q^{2}\mathrm{d}Y\mathrm{d}^{2}\vec{q}_{T}}$
$\displaystyle=\sigma_{0}\sum_{i,j}H_{ij}(Q,\mu)\\!\int\\!\\!\frac{\mathrm{d}^{2}\vec{b}_{T}}{(2\pi)^{2}}e^{\mathrm{i}\vec{q}_{T}\cdot\vec{b}_{T}}f_{i/h_{1}}\bigl{(}x_{1},\vec{b}_{T},\mu,\zeta_{1}\bigr{)}f_{j/h_{2}}\bigl{(}x_{2},\vec{b}_{T},\mu,\zeta_{2}\bigr{)},$
(2.15)

where $\mu$ is the renormalization scale, and the CS scales [51, 50] are

$\displaystyle\zeta_{1}=2(x_{1}P_{1}^{+})^{2}e^{-2y_{n}}\,,\qquad\zeta_{2}=2(x_{2}P_{2}^{-})^{2}e^{+2y_{n}}\,.$
(2.16)

Here, $y_{n}$ is an arbitrary scheme-dependent parameter that cancels in eq.
(2.15). In particular, we have that
$\zeta_{1}\zeta_{2}=(2x_{1}x_{2}P_{1}^{-}P_{2}^{+})^{2}=Q^{2}$.

The Collins scheme is characterized by spacelike Wilson lines with directions

$\displaystyle n_{A}^{\mu}(y_{A})$ $\displaystyle\equiv
n_{a}^{\mu}-e^{-2y_{A}}n_{b}^{\mu}$
$\displaystyle=(1,-e^{-2y_{A}},0_{\perp})\,,$ $\displaystyle
n_{B}^{\mu}(y_{B})$ $\displaystyle\equiv n_{b}^{\mu}-e^{2y_{B}}n_{a}^{\mu}$
$\displaystyle=(-e^{2y_{B}},1,0_{\perp})\,,$ (2.17)

parametrized by the rapidities $y_{A}$ and $y_{B}$. The Collins TMD for a
hadron $h$ moving along $n_{a}$ with rapidity $y_{P}$ is

$\displaystyle f^{C}_{i/h}(x,\vec{b}_{T},\mu,\zeta)=\lim_{\epsilon\to
0}Z_{\mathrm{uv}}^{R}(\epsilon,\mu,\zeta)\lim_{\begin{subarray}{c}y_{B}\to-\infty\end{subarray}}\frac{B^{C}_{i/h}(x,\vec{b}_{T},\epsilon,y_{P}-y_{B})}{\sqrt{S_{C}^{R}(b_{T},\epsilon,2y_{n},2y_{B})}}\,,$
(2.18)

where $B^{C}_{i/h}$ is the beam function and $S_{C}^{R}$ is the soft function.
$Z_{\mathrm{uv}}^{R}$ absorbs $\epsilon$-poles that result from working in
$d=4-2\epsilon$ dimensions to regulate UV divergences. $Z_{\rm uv}^{R}$ and
$S_{C}^{R}$ depend on the color representation $R$ of the parton $i$
(fundamental $R=q$ for quarks and adjoint $R=g$ for gluons) but are
independent of parton flavor. We emphasize that in eq. (2.18) the lightcone
limit $y_{B}\to-\infty$ is taken _before_ UV renormalization.

The Collins beam and soft functions for quarks and gluons are defined as

$\displaystyle B_{q_{i}/h}^{C}(x,\vec{b}_{T},\epsilon,y_{P}-y_{B})$
$\displaystyle=\int\\!\frac{\mathrm{d}b^{-}}{2\pi}e^{-\mathrm{i}b^{-}(xP^{+})}\Omega_{q_{i}/h}^{[\gamma^{+}]}\bigl{[}b,P,\epsilon,-\infty
n_{B}(y_{B}),b^{-}n_{b}\bigr{]}\,,$ $\displaystyle
B_{g/h}^{C\rho\sigma}(x,\vec{b}_{T},\epsilon,y_{P}-y_{B})$
$\displaystyle=\int\\!\frac{\mathrm{d}b^{-}}{2\pi}\frac{e^{-\mathrm{i}b^{-}(xP^{+})}}{xP^{+}}\Omega_{g/h}^{-\rho-\sigma}\bigl{[}b,P,\epsilon,-\infty
n_{B}(y_{B}),b^{-}n_{b}\bigr{]}\,,$ $\displaystyle
S_{C}^{R}(b_{T},\epsilon,y_{A},y_{B})$
$\displaystyle=S^{R}[b_{\perp},\epsilon,-\infty n_{A}(y_{A}),-\infty
n_{B}(y_{B})]\,,$ (2.19)

where $\Omega_{i/h}$ and $S^{R}$ are the correlators in eqs. (2.1) and (2.8).
The beam function path is

$\displaystyle b=(0,b^{-},b_{\perp})\,,\qquad\delta=(0,b^{-},0)\,,\qquad
v=n_{B}(y_{B})\quad{\rm with}\quad|\eta|\to\infty\,.$ (2.20)

This implies that $b-\delta=b_{\perp}$, and hence the Wilson line’s tranverse
segment is perpendicular to its longitudinal segments. This transverse segment
is important in singular gauges [70]. Note that the transverse segment is
often not specified in the literature: in nonsingular gauges such as Feynman
gauge, a Wilson line at lightcone infinity does not make contributions, and
its self-energy cancels against the corresponding piece in the soft function
[53]. The longitudinal Wilson line segments extend along $\eta
v\pm\tfrac{\delta}{2}=(\eta e^{2y_{B}},\eta\pm\tfrac{\delta}{2},0_{\perp})$
and thus only depend on the rapidity $y_{B}$ in the limit $|\eta|\to\infty$.
The limit $\eta\to-\infty$ taken in eq. (2.2.1) applies to Drell-Yan
kinematics, whereas SIDIS kinematics uses $\eta\to\infty$.

Finally, we remark that due to taking the lightcone limit prior to UV
renormalization, the Collins scheme is equivalent to schemes defined with
rapidity regulators on the lightcone [65, 58, 62, 60, 63, 64] that are often
employed in higher-order perturbative calculations and higher-order resummed
phenomenological analyses, see e.g. ref. [66] for a discussion.

####  2.2.2 Ji-Ma-Yuan (JMY) scheme

The JMY scheme [56] was introduced for the semi-inclusive deep-inelastic
scattering (SIDIS) process, $e^{-}p\to e^{-}hX$. The factorization theorem eq.
(2.2) for Drell-Yan-like processes takes the form

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}Q^{2}\mathrm{d}Y\mathrm{d}^{2}\vec{q}_{T}}$
$\displaystyle=\sigma_{0}\sum_{i,j}H_{ij}(Q,\mu,\rho)$ (2.21)
$\displaystyle\qquad\times\int\\!\frac{\mathrm{d}^{2}\vec{b}_{T}}{(2\pi)^{2}}\,e^{\mathrm{i}\vec{q}_{T}\cdot\vec{b}_{T}}\,f_{i/h_{1}}^{\mathrm{JMY}}\bigl{(}x_{1},\vec{b}_{T},\mu,x_{1}\zeta_{v},\rho\bigr{)}\,f_{j/h_{2}}^{\mathrm{JMY}}\bigl{(}x_{2},\vec{b}_{T},\mu,x_{2}\zeta_{\tilde{v}},\rho\bigr{)}\,.$

This scheme is characterized by timelike Wilson lines with directions

$\displaystyle v^{\mu}$
$\displaystyle=v^{+}n_{a}^{\mu}+v^{-}n_{b}^{\mu}=(v^{+},v^{-},0_{\perp})\,,$
$\displaystyle v^{-}$ $\displaystyle\gg v^{+}>0\,,$
$\displaystyle{\tilde{v}}^{\mu}$
$\displaystyle={\tilde{v}}^{+}n_{a}^{\mu}+{\tilde{v}}^{-}n_{b}^{\mu}=({\tilde{v}}^{+},{\tilde{v}}^{-},0_{\perp})\,,$
$\displaystyle\tilde{v}^{+}$ $\displaystyle\gg\tilde{v}^{-}>0\,.$ (2.22)

The definitions below always use these hierarchies, but not as strict limits;
for simplicity we leave them implicit. The offshellness of Wilson lines is
encoded in the parameters

$\displaystyle\zeta_{v}^{2}=\frac{(2P_{1}\cdot v)^{2}}{v^{2}}$
$\displaystyle=2(P_{1}^{+})^{2}\frac{v^{-}}{v^{+}}\,,\quad\zeta_{\tilde{v}}^{2}=\frac{(2P_{2}\cdot{\tilde{v}})^{2}}{{\tilde{v}}^{2}}=2(P_{2}^{-})^{2}\frac{{\tilde{v}}^{+}}{{\tilde{v}}^{-}}\,,\quad\rho^{2}=\frac{4(v\cdot{\tilde{v}})^{2}}{v^{2}{\tilde{v}}^{2}}=\frac{v^{-}{\tilde{v}}^{+}}{v^{+}{\tilde{v}}^{-}}\,.$
(2.23)

We define the JMY scheme TMD as

$\displaystyle
f^{\mathrm{JMY}}_{i/h}(x,\vec{b}_{T},\mu,x\zeta_{v},\rho)=\frac{B^{\mathrm{JMY}}_{i/h}(x,b_{T},\mu,\zeta_{v})}{\sqrt{S_{\mathrm{JMY}}^{R}(b_{T},\mu,\rho)}}\,.$
(2.24)

Here $B^{\mathrm{JMY}}_{i/h}$ and $S_{\mathrm{JMY}}$ are _renormalized_ beam
and soft functions. This is a crucial distinction from the Collins scheme, in
which we first combine the beam and soft functions, then take the lightlike
limit $\rho\to\infty$, and only thereafter carry out renormalization, cf. eq.
(2.18).

For a hadron moving in the $n_{a}$ direction, the JMY quark beam and soft
functions are

$\displaystyle B_{q_{i}/h}^{\mathrm{JMY}}(x,\vec{b}_{T},\mu,\zeta_{v})$
$\displaystyle=\int\\!\frac{\mathrm{d}b^{-}}{2\pi}e^{-\mathrm{i}b^{-}(xP^{+})}\Omega_{q_{i}/h}^{[\gamma^{+}]}\bigl{[}b,P,\mu,-\infty
v,b^{-}n_{b}\bigr{]}\,,$ $\displaystyle
S_{\mathrm{JMY}}^{R}(b_{T},\mu,y_{A},y_{B})$
$\displaystyle=S^{R}[b_{\perp},\mu,-\infty v,-\infty{\tilde{v}}]\,,$ (2.25)

where $\Omega_{i/h}$ and $S^{R}$ are the renormalized generic correlators in
eqs. (2.1) and (2.8). Just as in the Collins scheme, we take
$b^{\mu}=(0,b^{-},b_{\perp})$ and $\delta=(0,b^{-},0_{\perp})$ (note that this
is usually not specified in the literature). The JMY correlator differs from
the Collins correlator in eq. (2.2.1) by the presence of timelike ($v^{2}>0$)
rather than spacelike ($n_{B}^{2}<0$) Wilson lines. We can perturbatively
match the JMY and Collins scheme TMDs, see appendix B.

JMY gluon TMDs are not explicitly defined in the literature; nonetheless, one
can define them in an analogous manner to the quark TMDs through eq. (2.1).

####  2.2.3 Large Rapidity (LR) scheme

Finally, we introduce a new continuum TMD scheme, the LR scheme. The LR scheme
uses the same beam and soft functions as the Collins scheme, but a different
order of carrying out UV renormalization and approaching the lightcone.
Specifically,

$\displaystyle
f^{\mathrm{LR}}_{i/h}(x,\vec{b}_{T},\mu,\zeta,y_{P}-y_{B})=\lim_{\begin{subarray}{c}-y_{B}\gg
1\end{subarray}}\,\lim_{\epsilon\to
0}\,Z_{\mathrm{uv}}^{\mathrm{LR}}(\epsilon,\mu,y_{n}-y_{B})\frac{B^{C}_{i/h}(x,\vec{b}_{T},\epsilon,y_{P}-y_{B})}{\sqrt{S_{C}^{R}(b_{T},\epsilon,2y_{n},2y_{B})}}\,,$
(2.26)

where $\epsilon\to 0$ implements UV renormalization. The UV counterterm
$Z_{\rm uv}^{\mathrm{LR}}$ does not depend on the rapidities in the beam
function because the beam function does not involve nontrivial cusp angles,
and the renormalization of its staple-shaped Wilson line and quark operators
are rapidity-independent. It does however explicitly depend on the Wilson line
rapidity difference $y_{n}-y_{B}$ through the renormalization of the soft
function. In particular, the UV renormalization constant $Z_{\rm
uv}^{\mathrm{LR}}$ is the product of those of the bare beam function and soft
factor. As will be shown in section 3.1.2, the bare beam function at large,
but finite, $(-y_{B})$ is equal to the quasi-beam function in a hadron state
with momentum given by $y_{P}-y_{B}$. Using the auxiliary field formalism of
the Wilson lines, one can show [71, 72, 37, 73] that renormalization of the
quasi-beam function is multiplicative in coordinate space and independent of
the external hadron momentum. Therefore, $Z_{\rm uv}^{\mathrm{LR}}$ is
independent of $y_{P}-y_{B}$, and only depends on $y_{n}-y_{B}$ due to the
renormalization of the soft factor. This should be contrasted with the Collins
scheme, where $Z_{\rm uv}^{\rm C}$ depends on $\zeta\propto
e^{2(y_{P}-y_{n})}$, cf. eq. (2.18).

We note that the renormalized LR scheme TMD depends on $y_{P}-y_{B}$; hence,
the limit $-y_{B}\gg 1$ is to be understood as taking $(-y_{B})$ large but
finite instead of taking the limit $y_{B}\to-\infty$. The LR and Collins TMDs
use the same $\zeta=(2xm_{h}e^{y_{P}-y_{n}})^{2}$ to encode $y_{n}$
dependence. We can also view the LR scheme as the JMY scheme defined with
spacelike instead of timelike Wilson lines. See section 3.1 for further
elaboration on and derivation of LR scheme properties.

###  2.3 Lattice TMD matrix elements

Next, we provide a brief overview of TMD functions that are amenable to
calculation using lattice QCD. Unlike the continuum schemes in section 2.2,
lattice TMDs are defined using finite-length Wilson lines. We can obtain their
matrix elements from paths involving equal-time spacelike Wilson lines.

####  2.3.1 Quasi-TMDs

Quasi-TMDs are objects that share the same infrared physics as TMDs, but have
finite-length spacelike Wilson lines and are computable on the lattice [34,
36]. The general structure of a quasi-TMD looks quite similar to a TMD:

$\displaystyle\tilde{f}_{i/h}^{[s]}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z},\tilde{\eta})$
$\displaystyle=\lim_{a\to 0}Z^{\prime}_{\rm uv}(\mu,\tilde{\mu})\,Z_{\rm
uv}(a,\tilde{\mu},y_{n}-y_{B})$
$\displaystyle\qquad\quad\times\tilde{B}_{i/h}^{[s]}(x,\vec{b}_{T},a,\tilde{\eta},x\tilde{P}^{z})\,\tilde{\Delta}_{S}^{R}(b_{T},a,\tilde{\eta},y_{n},y_{B})\,,$
(2.27)

where $\tilde{B}_{i}$ is the quasi-beam function; $\tilde{\Delta}_{S}^{R}$ is
the quasi-soft factor; $Z_{\rm uv}(a,\tilde{\mu},\tilde{\zeta})$ implements a
lattice renormalization with the corresponding scale $\tilde{\mu}$;
$Z^{\prime}_{\rm uv}(\mu,\tilde{\mu},\tilde{\zeta})$ implements a conversion
to the ${\overline{\mathrm{MS}}}$ scheme with the ${\overline{\mathrm{MS}}}$
scale $\mu$; $\tilde{\eta}$ is the extent of the Wilson lines in the quasi-
beam and soft functions; and the dependence on $\tilde{\zeta}$ and the
rapidities $y_{n,B}$ of the Wilson lines is explained below. As
$|\tilde{\eta}|\to\infty$, the leading term in $\tilde{f}_{i/h}$ becomes
independent of $\tilde{\eta}$,

$\displaystyle\tilde{f}_{i/h}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z})\equiv\lim_{\tilde{\eta}\to\infty}\tilde{f}_{i/h}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z},\tilde{\eta})\,,$
(2.28)

with corrections of
$\mathcal{O}\bigl{[}{b_{T}}/{\tilde{\eta}},{1}/{(\tilde{P}^{z}\tilde{\eta})}\bigr{]}$.

The quasi-beam functions for quarks and gluons can be expressed using the
generic correlator in eq. (2.1) as

$\displaystyle\tilde{B}_{q_{i}/h}^{[\tilde{\Gamma}]}(x,\vec{b}_{T},a,\tilde{\eta},x\tilde{P}^{z})$
$\displaystyle=N_{\tilde{\Gamma}}\int\frac{\mathrm{d}\tilde{b}^{z}}{2\pi}\,e^{\mathrm{i}\tilde{b}^{z}(x\tilde{P}^{z})}\,\Omega_{q_{i}/h}^{[\tilde{\Gamma}]}(\tilde{b},\tilde{P},a,\tilde{\eta}\hat{z},\tilde{b}^{z}\hat{z})\,,$
$\displaystyle\tilde{B}_{g/h}^{\alpha\rho\beta\sigma}(x,\vec{b}_{T},a,\tilde{\eta},x\tilde{P}^{z})$
$\displaystyle=N^{\alpha\rho\beta\sigma}\int\frac{\mathrm{d}\tilde{b}^{z}}{2\pi}\,\frac{e^{\mathrm{i}\tilde{b}^{z}(x\tilde{P}^{z})}}{x\tilde{P}^{z}}\,\Omega_{g/h}^{\alpha\rho\beta\sigma}(\tilde{b},\tilde{P},a,\tilde{\eta}\hat{z},\tilde{b}^{z}\hat{z})\,.$
(2.29)

Here, $\tilde{\Gamma}$ is a Dirac structure, $\alpha$ and $\beta$ are generic
Lorentz indices, $\rho$ and $\sigma$ are transverse indices, and
$N_{\tilde{\Gamma}}$ and $N^{\alpha\rho\beta\sigma}$ are normalization
factors. To enable calculations on the lattice, we use equal-time paths

$\displaystyle\tilde{b}=(0,b_{T}^{x},b_{T}^{y},\tilde{b}^{z})\,,\qquad\tilde{\eta}v=\tilde{\eta}\hat{z}=(0,0,0,\tilde{\eta})\,,\qquad\delta=\tilde{b}^{z}\hat{z}=(0,0,0,\tilde{b}^{z})\,,$
(2.30)

as illustrated in figure -116. This choice of $\delta$ guarantees that the
transverse Wilson line segment is perpendicular to the longitudinal segments
and that the total length of the staple $\ell=2\tilde{\eta}+b_{T}$ is
independent of $\tilde{b}^{z}$. Therefore, the cusp angles in the quasi-beam
function are always trivially $\pi/2$, so that the UV renormalization factor
is independent of $\tilde{b}^{z}$ and can be pulled out of the Fourier
integral in eq. (2.3.1) [37]. This is a key difference to the MHENS scheme
discussed in section 2.3.2.

Figure -116: Wilson line structure of (left) the quasi-beam function in eq.
(2.3.1), and (right) the MHENS scheme in eq. (2.39). Quasi-TMD staple legs
extend along the $z$ direction and are closed by a perpendicular segment,
whereas MHENS staple legs extend along a generic spacelike direction $v^{\mu}$
and are closed by a segment with nontrivial cusp angle $\gamma$.

The construction of the quasi-beam function in eq. (2.3.1) is guided by the
observation that by boosting the operator and taking
$|\tilde{\eta}|\to\infty$, one recovers eq. (2.2.1). In contrast, the soft
function depends on two almost-lightlike directions and thus cannot be
obtained by boosting an equal-time operator [36]; several potential
(quasi-)soft functions have been proposed [33, 34, 36] despite the fact that
these proposals cannot recover eq. (2.2.1) under any Lorentz boost. Here, we
construct the lattice soft function as a finite-length version of the Collins
soft function in eq. (2.2.1),

$\displaystyle\tilde{S}^{R}(b_{T},a,\tilde{\eta},y_{A},y_{B})=S^{R}\bigg{[}b_{\perp},a,-\tilde{\eta}\frac{n_{A}(y_{A})}{|n_{A}(y_{A})|},-\tilde{\eta}\frac{n_{B}(y_{B})}{|n_{B}(y_{B})|}\bigg{]}\,.$
(2.31)

Here the length of the soft function path is $L_{\leavevmode\hbox
to4.69pt{\vbox to7.65pt{\pgfpicture\makeatletter\hbox{\hskip
0.2134pt\lower-0.327pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {}{{}}{}{{}}{}{{}}{}
{}{}{{}}{} {}{}{{}}{}{{}}{} {}{}{{}}{}{{}}{}{{}}{} {}{}{{}}{} {}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.42679pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{4.26773pt}{2.84544pt}\pgfsys@lineto{0.0pt}{5.69089pt}\pgfsys@lineto{0.0pt}{7.1136pt}\pgfsys@lineto{4.26773pt}{4.26816pt}\pgfsys@lineto{0.0pt}{1.42271pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.34143pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{-0.15628pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}=2(2\tilde{\eta}+b_{T})=2L_{\rm
staple}$. The choice of the minus sign in the last two arguments, $\eta v$ and
$\bar{\eta}\bar{v}$, allows Lorentz-invariant products obtained from these
choices to be more easily related to continuum schemes.

Combining eqs. (2.3.1) and (2.31) as required by eq. (2.3.1) gives

$\displaystyle\tilde{f}_{i/h}^{[\tilde{\Gamma}]}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z})$
$\displaystyle=\lim_{\begin{subarray}{c}\tilde{\eta}\to\infty\\\ a\to
0\end{subarray}}Z^{\prime}_{\rm uv}(\mu,\tilde{\mu})Z_{\rm
uv}(a,\tilde{\mu},y_{n}-y_{B})\,\frac{\tilde{B}_{i/h}^{[\tilde{\Gamma}]}(x,\vec{b}_{T},a,\tilde{\eta},x\tilde{P}^{z})}{\sqrt{\tilde{S}^{R}(b_{T},a,\tilde{\eta},2y_{n},2y_{B})}}$
$\displaystyle=\lim_{\tilde{\eta}\to\infty}\frac{\tilde{B}_{i/h}^{[\tilde{\Gamma}]}(x,\vec{b}_{T},\mu,\tilde{\eta},x\tilde{P}^{z})}{\sqrt{\tilde{S}^{R}(b_{T},\mu,\tilde{\eta},2y_{n},2y_{B})}}\,.$
(2.32)

Here
$\tilde{\zeta}=\big{(}xm_{h}e^{y_{\tilde{P}}+y_{B}-y_{n}}\big{)}^{2}=(2x\tilde{P}^{z}e^{y_{B}-y_{n}})^{2}$,
and the second equality holds for large $\tilde{P}^{z}$.

In practice, calculating (quasi-)TMD soft functions poses a significant
challenge for the lattice. It is possible to construct the quasi-soft function
indirectly through the spacelike meson form factor and quasi-wavefunction
[38]; promising first results using this approach have been reported in refs.
[46, 47].

Prior to this work, the literature has studied different proposals of the
quasi-soft function which are constructed from equal-time Wilson lines [33,
34, 35, 36, 38, 39]. The naive quasi-soft function features a rectangle-shaped
Wilson loop along the $z$ direction,

$\displaystyle\tilde{S}^{R}_{\rm naive}(b_{T},a,\tilde{\eta})$
$\displaystyle\equiv
S^{R}\bigl{[}b_{\perp},a,\tilde{\eta}\hat{z},-\tilde{\eta}\hat{z}\bigr{]}\,,$
(2.33)

whose renormalized continuum version with $\tilde{\eta}=\infty$ in the
${\overline{\mathrm{MS}}}$ scheme is denoted

$\displaystyle\tilde{S}^{R}_{\rm naive}(b_{T},\mu)$ $\displaystyle\equiv
S^{R}\bigl{[}b_{\perp},\mu,\infty\hat{z},-\infty\hat{z}\bigr{]}\,.$ (2.34)

However, it has been shown at one-loop level [36] that $\tilde{S}^{R}_{\rm
naive}(b_{T},\mu)$ does not have the correct IR physics for the quasi-TMD to
be perturbatively matchable to the Collins TMD. Although Refs. [34, 36]
proposed a bent quasi-soft function that works at one-loop order, it was
argued that the factorization utilizing this function will break down at two
loops [38].

Nevertheless, the naive quasi-soft function $\tilde{S}^{R}_{\rm naive}$ can
still serve a useful purpose for lattice calculations, where it can be used to
cancel linear power divergences proportional to $\tilde{\eta}$ and $b_{T}$.
For this reason it is useful to define the naive quasi-TMD as

$\displaystyle\tilde{f}_{i/h}^{[\tilde{\Gamma}]\rm
naive}(x,\vec{b}_{T},\mu,x\tilde{P}^{z})=\lim_{\tilde{\eta}\to\infty}\frac{\tilde{B}_{i/h}^{[\tilde{\Gamma}]}(x,\vec{b}_{T},\mu,\tilde{\eta},x\tilde{P}^{z})}{\sqrt{\tilde{S}^{R}_{\rm
naive}(b_{T},\mu,\tilde{\eta})}}\,.$ (2.35)

This can then be compared to the quasi-TMD $\tilde{f}_{i/h}$ defined with a
quasi-soft function $\tilde{S}^{R}$ that yields the correct infrared
structure. Since the $\eta$-dependence cancels exactly between the quasi beam
and soft functions, we have

$\displaystyle\lim_{\tilde{\eta}\to\infty}\frac{\tilde{S}^{R}_{\rm
naive}(b_{T},\mu,\tilde{\eta})}{\tilde{S}^{R}(b_{T},\mu,\tilde{\eta},2y_{n},2y_{B})}=\frac{\tilde{S}^{R}_{\text{naive}}(b_{T},\mu)}{S_{C}^{R}(b_{T},\mu,2y_{n},2y_{B})}\,,$
(2.36)

which leads to the relationship between $\tilde{f}^{\rm naive}_{i}$ and
$\tilde{f}_{i}$ in eq. (1.3). A further discussion of the equivalence in eq.
(2.36) is provided in section 3.1.2.

Moreover, at large rapidity the Collins soft function behaves as [50, 51, 52]

$\displaystyle\lim_{-y_{B}\gg 1}S^{R}_{\mathrm{C}}(b_{T},\mu,2y_{n},2y_{B})$
$\displaystyle=S_{I}(b_{T},\mu)\>e^{2(y_{n}-y_{B})\gamma_{\zeta}^{q}(b_{T},\mu)}\,,$
(2.37)

where $S_{I}(b_{T},\mu)$ is a rapidity-independent component of the soft
function. Although the entire $S^{R}_{\mathrm{C}}$ exponentiates due to the
non-Abelian exponentiation theorem for Wilson line operators [74, 75], the
important aspect of eq. (2.37) is the particular dependence on rapidity. Using
eq. (2.37), the multiplicative factor relating $\tilde{f}^{\rm naive}_{i}$ and
$\tilde{f}_{i}$ can be simplified as

$\displaystyle\lim_{\tilde{y}_{B}\to-\infty}\sqrt{\frac{\tilde{S}^{R}_{\text{naive}}(b_{T},\mu)}{S_{C}^{R}(b_{T},\mu,2y_{n},2y_{B})}}$
$\displaystyle=\sqrt{\frac{\tilde{S}^{R}_{\text{naive}}(b_{T},\mu)}{S_{I}(b_{T},\mu)}}e^{-2(y_{n}-y_{B})\gamma_{\zeta}^{q}(b_{T},\mu)}$
$\displaystyle\equiv\left[g_{S}^{q}(b_{T},\mu)\right]^{-1}\exp\left[-{1\over
2}\gamma_{\zeta}^{q}(b_{T},\mu)\ln{(2x\tilde{P}^{z})^{2}\over\zeta}\right]\,,$
(2.38)

Plugging eq. (2.3.1) into eq. (1.1) gives the original factorization formula
eq. (1.1) proposed in refs. [35, 36, 38, 39]. We identify
$g_{S}^{q}(b_{T},\mu)$ as the same factor introduced in ref. [36], which is
equivalent to the square root of the _reduced soft function_
$S_{r}(b_{T},\mu)$ in refs. [38, 39].

####  2.3.2 Musch-Hägler-Engelhardt-Negele-Schäfer (MHENS) scheme

TMDs were first studied on the lattice in refs. [27, 28, 29, 30, 31, 32] using
a Lorentz-invariant approach. We discuss this scheme as formulated in ref.
[29], and name it after the authors as the MHENS scheme. The goal of this
scheme is to calculate the beam function

$\displaystyle
B_{q_{i}/h}^{\mathrm{MHENS}~{}[\Gamma]}(x,\vec{b}_{T},P,a,\eta,v)$
$\displaystyle=N_{\Gamma}\int\frac{\mathrm{d}b^{-}}{2\pi}e^{-\mathrm{i}x(P\cdot
b)}\Omega_{q_{i}/h}^{[\Gamma]}(b,P,a,\eta v,0)\bigg{|}_{b^{+}=0}$ (2.39)

using lattice QCD. Here, $\Omega_{q/h}^{[\Gamma]}$ is our usual correlator
defined in eq. (2.1), supplemented with the special choice $\delta=0$ that
reduces the number of Lorentz invariants in eq. (2.1) to six. Note that the
MHENS scheme is based on a Lorentz-invariant formulation of the TMD correlator
in the $\delta=0$ case. Their correlator $\Phi$ is related to our $\Omega$ by

$\displaystyle\tilde{\Phi}^{[\Gamma]}_{\rm unsubtr.}(b,P,a,\eta v)$
$\displaystyle=\Omega_{q_{i}/h}^{[\Gamma]}(b,P,a,\eta v,0)\,.$ (2.40)

To make the connection between Minkowksi and Euclidean spaces for the
correlator in eq. (2.39), ref. [29] decomposes it for generic choices of $b,P$
and $v$ as

$\displaystyle\frac{1}{2}\Omega_{q_{i}/h}^{[\gamma^{\mu}]}(b,P,a,\eta v,0)$
$\displaystyle=P^{\mu}\tilde{A}_{2}\Bigl{(}b^{2},b\cdot P,\frac{v\cdot
b}{v\cdot P},\frac{v^{2}}{(v\cdot P)^{2}},\eta v\cdot P\Bigr{)}$
$\displaystyle\quad+\frac{P^{2}v^{\mu}}{v\cdot
P}\tilde{B}_{1}\Bigl{(}b^{2},b\cdot P,\frac{v\cdot b}{v\cdot
P},\frac{v^{2}}{(v\cdot P)^{2}},\eta v\cdot P\Bigr{)}+\ldots\,,$ (2.41)

where the ellipses stand for other spin-dependent and higher-twist structures.
For the full parametrization and similar decompositions for all $\Gamma$, see
ref. [29].

The scalar amplitudes $\tilde{A}_{2}$ and $\tilde{B}_{1}$ depend on all six
Lorentz invariants, though we leave implicit dependence on the hadron mass
$P^{2}=m_{h}^{2}$. For the beam function in eq. (2.39), it remains to specify
$b^{+}=0$ and $v_{\perp}=P_{\perp}=0$, which in a Lorentz-invariant fashion
reads [29]

$\displaystyle\frac{v\cdot b}{v\cdot P}=\frac{b\cdot
P}{P^{2}}R\bigl{(}\hat{\zeta}^{2}\bigr{)}\,,\quad
R\bigl{(}\hat{\zeta}^{2}\bigr{)}=1-\sqrt{1+\hat{\zeta}^{-2}}\,.$ (2.42)

Thus, the third argument of $\tilde{A}_{2}$ and $\tilde{B}_{1}$ in eq. (2.3.2)
is not independent of the other arguments, and on the lattice one is forced to
choose $b$ and $v$ such that eq. (2.42) is fulfilled. The CS-like parameter
$\hat{\zeta}$ entering eq. (2.42) is defined as

$\displaystyle\hat{\zeta}=\frac{v\cdot P}{\sqrt{|v^{2}|P^{2}}}\,.$ (2.43)

Parameterizing $P$ and $v$ by their rapidities shows that $\hat{\zeta}$
essentially is a rapidity difference,

$\displaystyle
P^{\mu}=\frac{m_{h}}{\sqrt{2}}(e^{y_{P}},e^{-y_{P}},0_{\perp})\,,\quad
v^{\mu}\propto(e^{y_{v}},-e^{-y_{v}},0_{\perp})\quad\Rightarrow\quad\hat{\zeta}=\sinh(y_{v}-y_{P})\,.$
(2.44)

Note the minus sign for $v^{-}$ is required for spacelike $v^{2}<0$; thus,
$y_{v}$ is not an actual rapidity. To connect to the continuum TMDs one
considers the large rapidity limit, equivalent to $v\cdot P\to\infty$ with
$P^{2}$ fixed, so $\hat{\zeta}\to\infty$.

Inserting eq. (2.42) into eq. (2.3.2) and specifying $\Gamma=\gamma^{+}$ as
required for the unpolarized TMD, at leading twist we obtain

$\displaystyle\frac{1}{2}\Omega_{q_{i}/h}^{[\gamma^{+}]}(b,P,a,\eta
v,0)\Big{|}_{b^{+}=0}$ $\displaystyle=P^{+}\tilde{A}_{2B}\Bigl{[}b^{2},b\cdot
P,\frac{v\cdot b}{v\cdot P}=\frac{b\cdot
P}{P^{2}}R\bigl{(}\hat{\zeta}^{2}\bigr{)},\frac{v^{2}}{(v\cdot P)^{2}},\eta
v\cdot P\Bigr{]}\,,$ (2.45)

where $\tilde{A}_{2B}=\tilde{A}_{2}+R(\hat{\zeta}^{2})\tilde{B}_{1}$.
Inserting this result into eq. (2.39), the unpolarized MHENS beam function is
given by

$\displaystyle\frac{1}{2}B_{q_{i}/h}^{\mathrm{MHENS}~{}[\gamma^{+}]}(x,\vec{b}_{T},P,a,\eta
v)$
$\displaystyle=\int\frac{\mathrm{d}(b{\cdot}P)}{2\pi}e^{-\mathrm{i}x(P{\cdot}b)}\tilde{A}_{2B}\Bigl{[}b^{2},b{\cdot}P,\frac{b{\cdot}P}{P^{2}}R\bigl{(}\hat{\zeta}^{2}\bigr{)},\frac{v^{2}}{(v{\cdot}P)^{2}},\eta
v{\cdot}P\Bigr{]}.$ (2.46)

One can obtain other leading-twist and spin-dependent beam functions in a
similar fashion. Due to its Lorentz-invariant formulation, eq. (2.46) can be
evaluated in any frame so long as eq. (2.42) is satisfied. This includes a
frame where $b^{0}=v^{0}=0$, as required on the lattice.

The above method has been applied to calculating the moments of TMDs. Since
integrating over $x$ sets $P\cdot b=0$, the moments can be calculated in a
frame where $b=(0,0,b_{\perp})$. In this case the beam functions
$\Omega_{q/h}^{[\gamma^{+}]}$ agree between the MHENS and Collins schemes
since we have $\delta=0$ in both cases. We elaborate on this further in
section 3.2.4.

To calculate the $x$-dependence, one must evaluate eq. (2.46) for generic
$P\cdot b\neq 0$. In this case, the MHENS and Collins staples are shaped
differently, and thus their beam functions are _not_ equivalent. There are two
key differences. First, the Collins staple is closed along
$b-\delta=b_{\perp}$, and thus its tranverse and longitudinal Wilson line
segments are perpendicular to one another. In contrast, the MHENS staple
closes along $b-\delta=b$, which for $v\cdot b\neq 0$ induces a $(P\cdot
b)$-dependent cusp angle according to eq. (2.42). Second, in the frame where
$b^{0}=0$, the MHENS staple length is not $b^{z}$-independent, leading to
nontrivial Wilson-line self-energies that depend on $b^{z}$ or $P\cdot b$.
Overall, this leads to a nontrivial Wilson line renormalization that depends
on $P\cdot b$ which cannot be factored out of the Fourier integral. We
demonstrate this at one-loop order in appendix C, and discuss the relation
between the MHENS and Collins TMDs in depth in section 3.2.4.

##  3 Factorization between physical and lattice TMDs

This section proves the factorization between the quasi- and Collins TMDs, eq.
(1.1). Our proof of factorization takes two steps, as depicted in figure -119:
first connecting the quasi-TMD to the intermediate LR scheme, then connecting
the LR and Collins schemes.

We begin with a bird’s eye view of various TMD schemes and their relationships
to one another. By using the general correlator $\Omega$ that we introduced in
eq. (2.1), in table 1 we see that TMDs take on a similar form in all schemes.
Key differences manifest in the order of limits taken to form the TMD, as well
as the specific arguments of the beam function correlator $\Omega$. Notably,
we can express $\Omega$ in terms of Lorentz-invariant combinations of its
arguments. Comparing the values that these Lorentz invariants take on in each
scheme, as shown in table 2, provides a useful way of relating different
schemes to one another. These tables are central to our proof, which we
present in section 3.1. We discuss implications of our results in section 3.2.

| TMD | Beam function | Soft function
---|---|---|---
Collins | $\lim\limits_{\epsilon\to 0}Z^{R}_{\text{UV}}\lim\limits_{y_{B}\to-\infty}\cfrac{\Omega_{i/h}}{\sqrt{S^{R}}}$ | $\Omega^{[\gamma^{+}]}_{q/h}\left[b,P,\epsilon,-\infty n_{B}(y_{B}),b^{-}n_{b}\right]$ | $S^{R}\left[b_{\perp},\epsilon,-\infty n_{A}(y_{A}),-\infty n_{B}(y_{B})\right]$
LR | $\lim\limits_{-y_{B}\gg 1}\lim\limits_{\epsilon\to 0}Z^{R}_{\text{UV}}\cfrac{\Omega_{i/h}}{\sqrt{S^{R}}}$ | $\Omega^{[\gamma^{+}]}_{q/h}\left[b,P,\epsilon,-\infty n_{B}(y_{B}),b^{-}n_{b}\right]$ | $S^{R}\left[b_{\perp},\epsilon,-\infty n_{A}(y_{A}),-\infty n_{B}(y_{B})\right]$
JMY | $\lim\limits_{\frac{v^{-}}{v^{+}}\gg 1}\lim\limits_{\epsilon\to 0}Z^{R}_{\text{UV}}\cfrac{\Omega_{i/h}}{\sqrt{S^{R}}}$ | $\Omega^{[\gamma^{+}]}_{q/h}\left[b,P,\mu,-\infty v,b^{-}n_{b}\right]$ | $S^{R}\left[b_{\perp},\mu,-\infty v,-\infty\tilde{v}\right]$
Quasi | $\lim\limits_{a\to 0}Z_{\text{UV}}\cfrac{B_{i/h}}{\sqrt{\tilde{S}^{R}}}$ | $\Omega_{q/h}^{[\gamma^{0,z}]}(\tilde{b},\tilde{P},a,\tilde{\eta}\hat{z},\tilde{b}^{z}\hat{z})$ | $S^{R}\left[b_{\perp},a,-\tilde{\eta}\cfrac{n_{A}(y_{A})}{|n_{A}(y_{A})|},-\tilde{\eta}\cfrac{n_{A}(y_{A})}{|n_{A}(y_{A})|}\right]$
MHENS |  | $\Omega^{[\Gamma]}_{q/h}(b,P,a,\eta v,0)$ |
Table 1: Overview of TMD schemes, as presented in section 2. The correlator
$\Omega$ is a function of Lorentz invariants constructed from its arguments.
See Table. 2 for a comparison of parameter values, Wilson line definitions,
and Lorentz invariants in each scheme.

###  3.1 Proof

We now present a proof of the quasi-to-Collins TMD factorization in detail for
the unpolarized quark TMD case. The proof of factorization for other leading-
twist TMDs follows naturally using the same framework, with only minor,
straightforward modifications for gluon TMDs or other spin structures. We also
remark that our proof employs dimensional regularization and the same UV
regulator for the quasi and LR schemes.

We begin our proof in section 3.1.1 by considering the correlator $\Omega$ as
a function of Lorentz invariants, and examining the values that these Lorentz
invariants take on in various schemes. In section 3.1.2, we see that the quasi
and LR scheme Lorentz invariants are identical at large proton momenta by
evaluating the quasi-TMD in a boosted frame. We thus can move from the quasi
to the LR scheme through a large rapidity expansion. In section 3.1.3 we
demonstrate that reversing the renormalization and lightcone limits to go from
the LR to the Collins scheme gives rise to a perturbative matching
coefficient. The combination of expansion and matching leads to the desired
factorization relation.

####  3.1.1 Beam correlators as a function of Lorentz invariants

Let us begin by examining the structure of the quasi-TMD. In dimensional
regularization, the quark quasi-beam function in eq. (2.3.1) reads

$\displaystyle\tilde{B}_{q_{i}/h}^{[\tilde{\Gamma}]}(x,\vec{b}_{T},\epsilon,\tilde{\eta},x\tilde{P}^{z})$
$\displaystyle=N_{\tilde{\Gamma}}\int\frac{\mathrm{d}\tilde{b}^{z}}{2\pi}\,e^{\mathrm{i}\tilde{b}^{z}(x\tilde{P}^{z})}\Omega_{q_{i}/h}^{[\tilde{\Gamma}]}\bigl{(}\tilde{b},\tilde{P},\epsilon,\tilde{\eta}\hat{z},\tilde{b}^{z}\hat{z}\bigr{)}\,,$
(3.1)

where $\tilde{b}^{\mu}=(0,\vec{b}_{T},\tilde{b}^{z})$. To study an unpolarized
Collins TMD, we must set $\Gamma=\gamma^{+}$ in eq. (2.2.1). To compare this
to the quasi-TMD, we must take $\tilde{\Gamma}=\gamma^{0}$ or $\gamma^{z}$,
which require normalization factors

$\displaystyle N_{\gamma^{z}}=1\,,\qquad\qquad
N_{\gamma^{0}}={\tilde{P}^{z}\over\tilde{P}^{0}}=\tanh(y_{\tilde{P}})\stackrel{{\scriptstyle
y_{\tilde{P}}\gg 1}}{{=}}1\,.$ (3.2)

We can decompose the coordinate-space correlator with arbitrary $b,P,v$ and
$\delta$ into Lorentz-covariant structures as333For the full parameterization
including spin-dependent terms, see e.g. ref. [29]. Note however that they
work with the correlator $\Omega$ where $\delta=0$. The more general analysis
carried out with our $\Omega$ at $\delta\neq 0$ gives rise to additional
terms.

$\displaystyle\Omega_{q_{i}/h}^{[\gamma^{\mu}]}(b,P,\epsilon,\tilde{\eta}v,\delta)$
$\displaystyle=P^{\mu}\Omega_{q_{i}/h}+{b^{\mu}\over-b^{2}}\Omega_{q_{i}/h}^{b}+{v^{\mu}\sqrt{P^{2}}\over\sqrt{|v|^{2}}}\Omega_{q_{i}/h}^{v}+{\delta^{\mu}\over-b^{2}}\Omega_{q_{i}/h}^{\delta}$
$\displaystyle=P^{\mu}\Omega_{q_{i}/h}+\mathrm{higher~{}twist}\,,$ (3.3)

where the dimensionless form factors $\Omega$ on the right-hand side are
functions of the 10 Lorentz invariants in eq. (2.1), which we suppress for
brevity. The prefactors share the same mass dimension and are finite as
$\delta\to 0$ or $b\cdot v/\sqrt{|v^{2}b^{2}|}\to 0$. In the second line, we
neglect terms that are suppressed at large momentum $P$, which do not
contribute at leading power.

| Collins / LR | JMY | Quasi | MHENS
---|---|---|---|---
$b^{\mu}$ | $(0,b^{-},b_{\perp})$ | $(0,b^{-},b_{\perp})$ | $(0,b_{T}^{x},b_{T}^{y},\tilde{b}^{z})$ | $(0,b_{T}^{x},b_{T}^{y},\tilde{b}^{z})$
$v^{\mu}$ | $(-e^{2y_{B}},1,0_{\perp})$ | $(v^{-}e^{2y_{B}^{\prime}},v^{-},0_{\perp})$ | $(0,0,0,-1)$ | $(0,v^{x},v^{y},v^{z})$
$\delta^{\mu}$ | $(0,b^{-},0_{\perp})$ | $(0,b^{-},0_{\perp})$ | $(0,0,0,\tilde{b}^{z})$ | $(0,0,0_{\perp})$
$P^{\mu}$ | $\dfrac{m_{h}}{\sqrt{2}}(e^{y_{P}},e^{-y_{P}},0_{\perp})$ | $\dfrac{m_{h}}{\sqrt{2}}(e^{y_{P}},e^{-y_{P}},0_{\perp})$ | $m_{h}(\cosh y_{\tilde{P}},0,0,\sinh y_{\tilde{P}})$ | $m_{h}\Bigl{(}\cosh y_{P},\dfrac{P^{x}}{m_{h}},\dfrac{P^{y}}{m_{h}},\sinh y_{P}\Bigr{)}$
$b^{2}$ | $-b_{T}^{2}$ | $-b_{T}^{2}$ | $-b_{T}^{2}-(\tilde{b}^{z})^{2}$ | $-b_{T}^{2}-(\tilde{b}^{z})^{2}$
$(\eta v)^{2}$ | $-2\eta^{2}e^{2y_{B}}$ | $2\eta^{2}(v^{-})^{2}e^{2y^{\prime}_{B}}$ | $-\tilde{\eta}^{2}$ | $-\eta^{2}\vec{v}^{\,2}$
$P\cdot b$ | $\dfrac{m_{h}}{\sqrt{2}}b^{-}e^{y_{P}}$ | $\dfrac{m_{h}}{\sqrt{2}}b^{-}e^{y_{P}}$ | $-m_{h}\tilde{b}^{z}\sinh y_{\tilde{P}}$ | $m_{h}\sinh y_{P}\tilde{b}^{z}+P^{x}b^{x}_{T}+P^{y}b^{y}_{T}$
x $\dfrac{b\cdot(\eta v)}{\sqrt{|(\eta v)^{2}b^{2}|}}$ | $-\dfrac{b^{-}e^{y_{B}}}{\sqrt{2}\,b_{T}}\,{\rm sgn}(\eta)$ | $\dfrac{b^{-}e^{y_{B}^{\prime}}}{\sqrt{2}\,b_{T}}\,{\rm sgn}(\eta)$ | $\dfrac{\tilde{b}^{z}}{\sqrt{(\tilde{b}^{z})^{2}+b_{T}^{2}}}\,{\rm sgn}(\eta)$ | $\dfrac{b_{T}^{x}v^{x}+b_{T}^{y}v^{y}+\tilde{b}^{z}v^{z}}{\sqrt{v_{T}^{2}+(v^{z})^{2}}\sqrt{b_{T}^{2}+(\tilde{b}^{z})^{2}}}$
x $\dfrac{P\cdot(\eta v)}{\sqrt{P^{2}|\eta v|^{2}}}$ | $\sinh(y_{P}\\!-\\!y_{B})\,{\rm sgn}(\eta)$ | $\cosh(y_{P}\\!-\\!y_{B}^{\prime})\,{\rm sgn}(\eta)$ | $\sinh y_{\tilde{P}}\,{\rm sgn}(\eta)$ | $\dfrac{P^{x}v^{x}+P^{y}v^{y}+m_{h}v^{z}\sinh y_{P}}{\sqrt{v_{T}^{2}+(v^{z})^{2}}\sqrt{m_{h}^{2}+P_{x}^{2}+P_{y}^{2}}}$
$\dfrac{\delta^{2}}{b^{2}}$ | $0$ | $0$ | $\dfrac{(\tilde{b}^{z})^{2}}{b_{T}^{2}+(\tilde{b}^{z})^{2}}$ | $0$
$\dfrac{b\cdot\delta}{b^{2}}$ | $0$ | $0$ | $\dfrac{(\tilde{b}^{z})^{2}}{b_{T}^{2}+(\tilde{b}^{z})^{2}}$ | $0$
$\dfrac{P\cdot\delta}{P\cdot b}$ | $1$ | $1$ | $1$ | $0$
$\cfrac{\delta\cdot(\eta v)}{b\cdot(\eta v)}$ | $1$ | $1$ | $1$ | $0$
$P^{2}$ | $m_{h}^{2}$ | $m_{h}^{2}$ | $m_{h}^{2}$ | $m_{h}^{2}$
Table 2: Overview of the Lorentz invariants entering the generic TMD
correlator as specified by eq. (2.1). Note that the Collins and LR schemes use
the same four-vectors.

Combining eqs. (3.1) and (3.1.1) and using $\tilde{P}_{\perp}=0$, we have

$\displaystyle\tilde{B}_{q_{i}/h}^{[\tilde{\Gamma}]}(x,\vec{b}_{T},\epsilon,\tilde{\eta},x\tilde{P}^{z})$
$\displaystyle=\int\frac{\mathrm{d}(\tilde{b}\cdot\tilde{P})}{2\pi}\,e^{-\mathrm{i}x(\tilde{b}\cdot\tilde{P})}\,\Omega_{q_{i}/h}\bigl{(}\tilde{b},\tilde{P},\epsilon,\tilde{\eta}\hat{z},\tilde{b}^{z}\hat{z}\bigr{)}\,.$
(3.4)

Note that the integration measure, the Fourier phase, and $\Omega_{q_{i}/h}$
are Lorentz invariants. We can write the LR/Collins beam function similarly:

$\displaystyle B_{q_{i}/h}^{C}(x,\vec{b}_{T},\epsilon,y_{P}-y_{B})$
$\displaystyle=\int\frac{\mathrm{d}(b\cdot P)}{2\pi}\,e^{-\mathrm{i}x(b\cdot
P)}\,\Omega_{q_{i}/h}\bigl{[}b,P,\epsilon,-\infty
n_{B}(y_{B}),b^{-}n_{b}\bigr{]}\,.$ (3.5)

The only differences between eqs. (3.4) and (3.5) lie in the beam function
parametrization, as well as the gauge link’s direction, length, and closure at
infinity, as seen in table 2:

Quasi-TMD:
$\displaystyle\tilde{b}^{\mu}=(0,b_{T}^{x},b_{T}^{y},\tilde{b}^{z})\,,\quad$
$\displaystyle\tilde{P}^{\mu}=\frac{m_{h}}{2}(e^{y_{\tilde{P}}},e^{-y_{\tilde{P}}},0_{\perp})\,,\quad$
$\displaystyle\tilde{\eta}\hat{z}^{\mu}=(0,0,0,\tilde{\eta})\,,$ (3.6) LR
scheme: $\displaystyle b^{\mu}=(0,b^{-},b_{\perp})\,,\quad$ $\displaystyle
P^{\mu}=\frac{m_{h}}{2}(e^{y_{P}},e^{-y_{P}},0_{\perp})\,,\quad$
$\displaystyle\eta n_{B}^{\mu}=\eta(-e^{2y_{B}},1,0_{\perp})\,.$

Note that we distinguish quasi components by tildes. Both schemes have the
same transverse components $b_{\perp}^{\mu}=(0,b_{T}^{x},b_{T}^{y},0)$.

####  3.1.2 Relating LR and quasi-TMDs

As we saw in section 2.1, we can express the quasi-, LR, and Collins TMDs in
terms of the same Lorentz-invariant function $\Omega$ in eq. (2.1), albeit
with different parametrizations for its arguments. We can relate the LR and
quasi-TMDs to one another by considering Lorentz-transforms of their
arguments. Boosting the quasi-TMD four-vectors in eq. (3.6) by a rapidity
$y_{B}$, we have

$\displaystyle\tilde{b}^{\mu}$
$\displaystyle=(0,b_{T}^{x},b_{T}^{y},\tilde{b}^{z})=\Bigl{(}\frac{\tilde{b}^{z}}{\sqrt{2}},-\frac{\tilde{b}^{z}}{\sqrt{2}},b_{\perp}\Bigr{)}$
$\displaystyle\quad\stackrel{{\scriptstyle\rm
boost}}{{\longrightarrow}}\quad\Bigl{(}\frac{\tilde{b}^{z}e^{y_{B}}}{\sqrt{2}},-\frac{\tilde{b}^{z}e^{-y_{B}}}{\sqrt{2}},b_{\perp}\Bigr{)}\,,$
$\displaystyle\tilde{P}^{\mu}$
$\displaystyle=\frac{m_{h}}{\sqrt{2}}(e^{y_{\tilde{P}}},e^{-y_{\tilde{P}}},0_{\perp})$
$\displaystyle\quad\stackrel{{\scriptstyle\rm
boost}}{{\longrightarrow}}\quad\frac{m_{h}}{\sqrt{2}}(e^{y_{\tilde{P}}+y_{B}},e^{-(y_{\tilde{P}}+y_{B})},0_{\perp})\,,$
$\displaystyle\tilde{\eta}\hat{z}$
$\displaystyle=(0,0,0,1)=\frac{\tilde{\eta}}{\sqrt{2}}(1,-1,0_{\perp})$
$\displaystyle\quad\stackrel{{\scriptstyle\rm
boost}}{{\longrightarrow}}\quad-\frac{\tilde{\eta}e^{-y_{B}}}{\sqrt{2}}(-e^{2y_{B}},1,0_{\perp})\,,$
$\displaystyle\delta^{\mu}$
$\displaystyle=\tilde{b}^{z}\hat{z}^{\mu}=(0,0,0,\tilde{b}^{z})=\frac{\tilde{b}^{z}}{\sqrt{2}}(1,-1,0_{\perp})$
$\displaystyle\quad\stackrel{{\scriptstyle\rm
boost}}{{\longrightarrow}}\quad-\frac{\tilde{b}^{z}e^{-y_{B}}}{\sqrt{2}}(-e^{2y_{B}},1,0_{\perp})\,,$
(3.7)

where we recall that $\delta$ encodes the geometry of the transverse Wilson
line. We use lightcone coordinates to make the boost manifest. Comparing eq.
(3.1.2) to eq. (3.6), we can match all but one component of the boosted-quasi
and LR schemes if we take

$\displaystyle b^{-}\equiv-\frac{\tilde{b}^{z}e^{-y_{B}}}{\sqrt{2}}\,,\qquad
y_{P}\equiv
y_{\tilde{P}}+y_{B}\,,\qquad\eta\equiv-\frac{\tilde{\eta}e^{-y_{B}}}{\sqrt{2}}\,;$
(3.8)

that is, except for $\tilde{b}^{+}$ of the boosted $\tilde{b}^{\mu}$.
Fortunately, this does not hold us back: to make the correspondence we need to
fix the parameters $b^{-}$, $y_{P}$, and $\eta$ to their LR scheme values in
eq. (3.8), whereas $\tilde{b}^{+}$ is a derived quantity. Additionally, in the
large rapidity limit we can neglect $\tilde{b}^{+}$:

$\displaystyle\tilde{b}^{+}\equiv\frac{\tilde{b}^{z}e^{y_{B}}}{\sqrt{2}}=-b^{-}e^{2y_{B}}\quad\stackrel{{\scriptstyle
y_{B}\to-\infty}}{{\longrightarrow}}\quad 0\,.$ (3.9)

In this limit, we also have that

$\displaystyle\tilde{b}^{z}$
$\displaystyle=-\sqrt{2}b^{-}e^{y_{B}}\quad\stackrel{{\scriptstyle
y_{B}\to-\infty}}{{\longrightarrow}}\quad 0\,,\qquad
y_{\tilde{P}}=y_{P}-y_{B}\quad\stackrel{{\scriptstyle
y_{B}\to-\infty}}{{\longrightarrow}}\quad\infty\,,$ (3.10)

so at large $\tilde{P}^{z}$ and small $\tilde{b}^{z}$ relative to $b_{T}$, the
quasi-correlator is equivalent to the LR correlator.

Note that we could have alternatively demonstrated the equivalence of the
quasi- and LR/Collins Lorentz invariants by transforming $\tilde{b}^{z}$,
$y_{P}$, and $\tilde{\eta}$ from the values in table 2 to those in eq. (3.8),
and then applying the limit $y_{B}\to-\infty$. For example,

$\displaystyle\tilde{P}\cdot\tilde{b}$ $\displaystyle=-m_{h}\tilde{b}^{z}\sinh
y_{\tilde{P}}$
$\displaystyle=m_{h}\sqrt{2}e^{y_{B}}b^{-}\sinh(y_{P}-\tilde{y}_{B})$
$\displaystyle\quad\stackrel{{\scriptstyle
y_{B}\to-\infty}}{{\longrightarrow}}\quad\frac{m_{h}}{\sqrt{2}}b^{-}e^{y_{P}}\,,$
$\displaystyle\frac{\delta^{2}}{\tilde{b}^{2}}$
$\displaystyle=\frac{1}{1+(b_{T}/\tilde{b}^{z})^{2}}$
$\displaystyle=\frac{1}{1+\left(\frac{b_{T}e^{-y_{B}}}{\sqrt{2}b^{-}}\right)^{2}}$
$\displaystyle\quad\stackrel{{\scriptstyle
y_{B}\to-\infty}}{{\longrightarrow}}\quad 0\,.$ (3.11)

By definition, the Lorentz invariants of a TMD remain unchanged by Lorentz
boosts. If we expand the quasi-TMD invariants in the boosted frame at large
$-y_{B}$ around those of the LR/Collins scheme, we can write a relationship
between correlators:

$\displaystyle\Omega_{q_{i}/h}(\tilde{b},\tilde{P},\epsilon,\tilde{\eta}\hat{z},\tilde{b}^{z}\hat{z})\Big{|}_{\mathrm{quasi}}$
$\displaystyle=\Omega_{q_{i}/h}(\tilde{b},\tilde{P},\epsilon,\tilde{\eta}\hat{z},\tilde{b}^{z}\hat{z})\Big{|}_{\mathrm{boosted~{}quasi}}\,,$
(3.12)
$\displaystyle\lim_{y_{B}\ll-1}\Omega_{q_{i}/h}(\tilde{b},\tilde{P},\epsilon,\tilde{\eta}\hat{z},\tilde{b}^{z}\hat{z})\Big{|}_{\mathrm{boosted~{}quasi}}$
$\displaystyle=\lim_{y_{B}\ll-1}\Omega_{q_{i}/h}(b,P,\epsilon,\eta
n_{B}(y_{B}),b^{-}n_{b})\Big{|}_{\mathrm{Collins/LR}}\,.$

Making the parameterizations of $b$ and $P$ in both schemes explicit and
shifting $y_{\tilde{P}}\to y_{P}-y_{B}$ (this is not a boost, but rather a
change in the parametrization of the proton’s momentum) we obtain

$\displaystyle\lim_{y_{B}\ll-1}\Omega_{q_{i}/h}\left[\tilde{b}\\!=\\!(0,\vec{b}_{T},-\sqrt{2}b^{-}e^{y_{B}}),\,\tilde{P}\\!=\\!\frac{m_{h}}{2}\bigl{(}e^{y_{P}-y_{B}},e^{-(y_{P}-y_{B})},0_{\perp}\bigr{)},\,\epsilon,\,\tilde{\eta}\hat{z},\,\tilde{b}^{z}\hat{z}\right]$
$\displaystyle\>=\lim_{y_{B}\ll-1}\\!\\!\Omega_{q_{i}/h}\left[b\\!=\\!(0,b^{-},b_{\perp}),\,P\\!=\\!\frac{m_{h}}{2}\bigl{(}e^{y_{P}},e^{-y_{P}},0_{\perp}\bigr{)},\,\epsilon,\,-\frac{\tilde{\eta}e^{-y_{B}}}{\sqrt{2}}n_{B}(y_{B}),\,b^{-}n_{b}\right].$
(3.13)

Here, the first correlator yields the quasi-beam function at the shifted
proton momentum, while the second correlator is that of the Collins/LR scheme
at _finite_ length

$\eta=-\frac{\tilde{\eta}e^{-y_{B}}}{\sqrt{2}}.$ (3.14)

Note that $\eta$ and $\tilde{\eta}$ always have opposite signs, and that
$\eta<0$ corresponds to the TMD PDF for Drell-Yan, while $\eta>0$ corresponds
to the TMD PDF for SIDIS.

Next, we supplement eq. (3.1.2) with a soft subtraction and UV
renormalization. On the lattice we cannot take the strict limit
$y_{B}\to-\infty$, so we must keep $y_{B}$ large but finite. The Collins
scheme entails taking the lightcone limit of $B/\sqrt{S}$ prior to UV
renormalization, but here we must renormalize at finite $y_{B}$. Up until this
point, all statements we made hold for both the bare Collins and LR schemes,
but for the remainder of this subsection, we only compare the renormalized
quasi- and LR TMDs. Let us now write the renormalized quasi- and LR TMDs as

$\displaystyle\phantom{=}\tilde{f}_{q_{i}/h}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z},\tilde{\eta})$
$\displaystyle\qquad=\int\frac{\mathrm{d}(\tilde{P}{\cdot}\tilde{b})}{2\pi}e^{-\mathrm{i}x(\tilde{P}{\cdot}\tilde{b})}\lim_{\epsilon\to
0}Z_{\rm
uv}^{q}(\mu,\epsilon,y_{n}-y_{B})\frac{\Omega_{q_{i}/h}(\tilde{b},\tilde{P},\epsilon,\tilde{\eta}\hat{z},\tilde{b}^{z}\hat{z})}{\sqrt{\tilde{S}^{q}(b_{T},\epsilon,\tilde{\eta},2y_{n},2y_{B})}}\,,$
(3.15)

and

$\displaystyle\phantom{=}f^{\mathrm{LR}}_{q_{i}/h}(x,\vec{b}_{T},\mu,\zeta,y_{P}-y_{B},\eta)$
$\displaystyle\qquad=\int\frac{\mathrm{d}(P{\cdot}b)}{2\pi}e^{-\mathrm{i}x(P{\cdot}b)}\lim_{\epsilon\to
0}Z_{\rm
uv}^{q}(\mu,\epsilon,y_{n}-y_{B})\frac{\Omega_{q_{i}/h}[b,P,\epsilon,\eta
n_{B}(y_{B}),b^{-}n_{b}]}{\sqrt{S^{q}(b_{T},\epsilon,\eta,2y_{n},2y_{B})}}\,.$
(3.16)

Here we define the five argument $S^{q}$ by

$\displaystyle S^{q}(b_{T},\epsilon,\eta,2y_{n},2y_{B})$
$\displaystyle=\tilde{S}^{q}(b_{T},\epsilon,\tilde{\eta},2y_{n},2y_{B})\,.$
(3.17)

The parameter $y_{n}$ governs the amount of soft radiation absorbed into the
TMDs and gives rise to the CS scales

$\displaystyle\tilde{\zeta}=2(x\tilde{P}^{+}e^{y_{B}-y_{n}})^{2}=x^{2}m_{h}^{2}e^{2(y_{\tilde{P}}+y_{B}-y_{n})}\,,\quad\zeta=2(xP^{+}e^{-y_{n}})^{2}=x^{2}m_{h}^{2}e^{2(y_{P}-y_{n})}\,.$
(3.18)

In eq. (3.1.2), following standard notation for quasi-TMDs, we encode
dependence on $y_{P}-y_{B}$ in $\tilde{P}^{z}=m_{h}\sinh(y_{P}-y_{B})$,
whereas in eq. (3.1.2) we state this dependence explicitly. The constraint of
$y_{\tilde{P}}=y_{P}-y_{B}$ leads to $\zeta=\tilde{\zeta}$. For both TMDs, we
use the finite-length soft function in eq. (2.31), repeated here for
convenience:

$\displaystyle\tilde{S}^{R}(b_{T},\epsilon,\tilde{\eta},y_{A},y_{B})=S^{R}\bigg{[}b_{\perp},\epsilon,-\tilde{\eta}\frac{n_{A}(y_{A})}{|n_{A}(y_{A})|},-\tilde{\eta}\frac{n_{B}(y_{B})}{|n_{B}(y_{B})|}\bigg{]}\,.$
(3.19)

The geometric length of the soft function Wilson line is twice of that of the
quasi-beam function, so that all linear divergences from Wilson line self-
energies cancel in eq. (3.1.2).444 Recall that in dimensional regularization
considered here, these linear divergences appear as poles in $1/(d-3)$, and
hence are absent in the ${\overline{\mathrm{MS}}}$ scheme where only poles in
$1/(d-4)$ are subtracted. Hence, these linear divergences are set to zero for
perturbative calculations in the ${\overline{\mathrm{MS}}}$ scheme, but it is
important to take them into account for a definition amenable to lattice
calculations.

Since the hadronic matrix elements in eqs. (3.1.2) and (3.1.2) are related by
a boost, we naturally also employ this soft function for the finite-length LR
scheme.

Finally, we discuss the form of the UV counterterm $Z^{q}$, which is simply
the ratio of the individual counterterms $Z_{\rm uv}^{B}$ and $Z_{\rm uv}^{S}$
for the beam and soft functions,

$\displaystyle Z_{\rm uv}^{q}(\mu,\epsilon,y_{n}-y_{B})$
$\displaystyle=\frac{Z_{\rm uv}^{B}(\mu,\epsilon)}{\sqrt{Z_{\rm
uv}^{S}(\mu,\epsilon,2y_{n}-2y_{B})}}\,.$ (3.20)

Here, we use that in the ${\overline{\mathrm{MS}}}$ scheme, the UV divergences
of the quasi-beam and soft functions are multiplicative and
$x\tilde{P}^{z}$-independent, according to the auxiliary field formalism [73,
40].

Using eq. (3.1.2), we can now relate the renormalized finite-length quasi-TMD
and LR TMD defined in eqs. (3.1.2) and (3.1.2),

$\displaystyle\lim_{y_{B}\ll-1}\tilde{f}_{q_{i}/h}(x,\vec{b}_{T},\mu,\tilde{\eta},\tilde{\zeta},x\tilde{P}^{z})=\lim_{y_{B}\ll-1}f^{\mathrm{LR}}_{q_{i}/h}\Bigl{(}x,\vec{b}_{T},\mu,-\frac{\tilde{\eta}}{2}e^{-y_{B}},\tilde{\zeta},y_{P}-y_{B}\Bigr{)}\,.$
(3.21)

Here, we have accounted for the change of Wilson line length in the LR scheme
and used $\tilde{\zeta}$ as the common CS scale.

The final step is to take the limit $\tilde{\eta}\to\infty$ to relate this
result to the continuum TMD. In the continuum TMD, this limit is taken prior
to UV renormalization (or $\epsilon\to 0$), while on the lattice one is forced
to extrapolate to infinite $\tilde{\eta}$ after renormalization. Thus, we must
show that the limits $\tilde{\eta}\to\infty$ and $\epsilon\to 0$ commute.
First, for $\tilde{\eta}\gg b_{T}\gg\tilde{b}^{z}$ the Wilson-line self-energy
contributions cancel exactly in the ratio $\Omega/\sqrt{S}$, which is not
affected by the order of the $\tilde{\eta}\to\infty$ and $\epsilon\to 0$
limits. The reason is that the staple geometry in the correlator $\Omega$ is
one half of that of the quasi-soft function $S$, and the exchange of gluons
between the two halves in $S$ is exponentially suppressed due to the spacelike
separation and large $\tilde{\eta}$. Second, after the subtraction of Wilson-
line self-energy diagrams, the remaining diagrams will include the eikonal
propagators

$\displaystyle\frac{1-e^{\pm i(n_{B}\cdot(k\pm
i\varepsilon)\tilde{\eta}/|n_{B}|}}{n_{B}\cdot k\pm i\varepsilon}\,,$ (3.22)

which have a singularity at $n_{B}\cdot k=0$. Such a singularity can be
regulated either with a finite $\tilde{\eta}$ and without the imaginary part
$i\varepsilon$, or if we keep $i\varepsilon$ and throw away the second term in
the numerator. For both regulators, the results are the same and independent
of $\tilde{\eta}$ if one takes the $\tilde{\eta}\to\infty$ or $\varepsilon\to
0$ limit in the end, which we have also verified explicitly at one-loop order.
Therefore, the $\tilde{\eta}\to\infty$ and $\epsilon\to 0$ limits also commute
for these diagrams.

In summary, this commutativity leads to the equivalence of the quasi- and LR
TMDs with infinite Wilson lines:

$\displaystyle\tilde{f}_{q_{i}/h}$
$\displaystyle(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z})$
$\displaystyle=\lim_{-y_{B}\gg
1}\int\frac{\mathrm{d}(\tilde{P}{\cdot}\tilde{b})}{2\pi}e^{-\mathrm{i}x(\tilde{P}{\cdot}\tilde{b})}\lim_{\epsilon\to
0}Z_{\rm
uv}^{q}(\mu,\epsilon,y_{n}-y_{B})\lim_{\tilde{\eta}\to\infty}\frac{\Omega_{q_{i}/h}(\tilde{b},\tilde{P},\epsilon,\tilde{\eta}\hat{z},\tilde{b}^{z}\hat{z})}{\sqrt{\tilde{S}^{q}(b_{T},\epsilon,\tilde{\eta},2y_{n},2y_{B})}}$
$\displaystyle=\lim_{-y_{B}\gg
1}\int\frac{\mathrm{d}(P{\cdot}b)}{2\pi}e^{-\mathrm{i}x(P{\cdot}b)}\lim_{\epsilon\to
0}Z_{\rm
uv}^{q}(\mu,\epsilon,y_{n}-y_{B})\lim_{\eta\to\infty}\frac{\Omega_{q_{i}/h}[b,P,\epsilon,\eta
n_{B}(y_{B}),b^{-}n_{b}]}{\sqrt{S^{q}(b_{T},\epsilon,\eta,2y_{n},2y_{B})}}$
$\displaystyle=f^{\mathrm{LR}}_{q_{i}/h}(x,\vec{b}_{T},\mu,\tilde{\zeta},y_{P}-y_{B})\,.$
(3.23)

Here we used the commutativity of the $\epsilon\to 0$ and
$\tilde{\eta}\to\infty$ limits in first step, used eq. (3.21) in the second
step, and used that the $\eta\to\infty$ limit naturally gives the continuum LR
TMD defined in eq. (2.26) in the last step.

####  3.1.3 Matching LR and Collins TMDs

We now know that the quasi and LR schemes are equivalent in the large
$\tilde{P}^{z}$, large rapidity, and $\tilde{\eta}\to\infty$ limits. The next
step is to derive the relation between the LR and Collins schemes. According
to eqs. (2.18) and (2.26), the only difference between these schemes is the
order of their $\epsilon\to 0$ and $y_{B}\to-\infty$ limits. In the LR scheme,
large $(-y_{B})$ corresponds to a momentum scale

$\displaystyle\zeta_{\mathrm{LR}}=4x^{2}m_{h}^{2}\sinh^{2}(y_{P}-y_{B})\,,$
(3.24)

so the limit $y_{B}\to-\infty$ corresponds to $\zeta_{\mathrm{LR}}\to\infty$.

Due to asymptotic freedom in QCD, changing the order of limits $\epsilon\to 0$
and $\zeta_{\mathrm{LR}}\to\infty$ should only affect the UV region while
leaving infrared (IR) physics intact. Using the LaMET formalism [12, 13, 14],
we can relate the two different orders of limits with a factorization formula
or perturbative matching,555The impact of exchanging these limits has also
been pointed out in the discussion of the Sudakov form factor in ref. [53], in
particular Eq. (10.97) therein. which takes the form

$\displaystyle f_{q_{i}/h}(x,\vec{b}_{T},\mu,\zeta)$
$\displaystyle=C_{q}^{-1}\big{(}\sqrt{\zeta_{\mathrm{LR}}}/2,\mu\big{)}\>f^{\mathrm{LR}}_{q_{i}/h}(x,\vec{b}_{T},\mu,\zeta,y_{P}-y_{B})+{\cal
O}(y_{B}^{k}e^{y_{B}})\,,$ (3.25)

where we have expanded at large $(-y_{B})$ or in $\Lambda_{\rm
QCD}^{2}/\zeta_{\mathrm{LR}}$. Here, $y_{B}^{k}$ captures logarithms of
$\zeta_{\mathrm{LR}}$ where $k$ is a positive integer. At large $(-y_{B})$,
the matching coefficient $C$ cancels the overall $y_{B}$ dependence in the
renormalized beam and soft functions through its $\zeta_{\mathrm{LR}}$
dependence. By dimensional analysis, $C_{q}$ could depend on the logarithm of
$\zeta_{\mathrm{LR}}/\mu^{2}$ or the rapidity difference $(y_{n}-y_{B})$.
However, according to eq. (2.37), the dependence on $(y_{n}-y_{B})$ implies
dependence on the CS kernel $\gamma_{\zeta}^{q}(b_{T},\mu)$, which at
$b_{T}\sim 1/\Lambda_{\rm QCD}$ is nonperturbative and hence infrared
sensitive. Therefore, $C_{q}$ can only depend on $\zeta_{\mathrm{LR}}/\mu^{2}$
and $\alpha_{s}(\mu)$, and so the factorization formula for the LR-TMD is

$\displaystyle f^{\mathrm{LR}}_{q_{i}/h}(x,\vec{b}_{T},\mu,\zeta,y_{P}-y_{B})$
$\displaystyle=C_{q}\big{(}\sqrt{\zeta_{\mathrm{LR}}}/2,\mu\big{)}\,f_{q_{i}/h}(x,\vec{b}_{T},\mu,\zeta)+{\cal
O}(y_{B}^{k}e^{y_{B}})\,,$ (3.26)

where $\zeta_{\mathrm{LR}}$ is given by eq. (3.24). Combining this with the
relation derived above in eqs. (1.3) and (2.3.1) for the quasi-TMD, with
$\zeta_{\mathrm{LR}}=(2x\tilde{P}^{z})^{2}$, we have

$\displaystyle\tilde{f}_{q_{i}/h}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z})$
$\displaystyle=C_{q}(x\tilde{P}^{z},\mu)f_{q_{i}/h}(x,\vec{b}_{T},\mu,\tilde{\zeta})+{\cal
O}(\tilde{y}_{P}^{k}e^{-\tilde{y}_{P}})$ (3.27)
$\displaystyle=C_{q}(x\tilde{P}^{z},\mu)\exp\left[\frac{1}{2}\gamma_{\zeta}^{q}(\mu,b_{T})\ln\frac{\tilde{\zeta}}{\zeta}\right]f_{q_{i}/h}(x,\vec{b}_{T},\mu,\zeta)+{\cal
O}(\tilde{y}_{P}^{k}e^{-\tilde{y}_{P}})\,.$

where $\zeta$ can be of any value. All the steps in our analysis work equally
well for gluons, simply replacing subscripts $q_{i}\to g$. Thus, this
completes our proof of the factorization formula in eq. (1.1).

###  3.2 Implications

Next, we discuss implications of the factorization relation. In section 3.2.1,
we show how the factorization implies that in the matching quarks and gluons
do not mix, nor do different quark flavors. In section 3.2.2 we derive the
momentum evolution of the quasi-beam function and hard matching coefficient.
In section 3.2.3 we discuss using ratios of quasi-TMDs or quasi-beam functions
to extract ratios of TMDs. In section 3.2.4 we examine implications of our
analysis for the MHENS scheme, outlining the additional steps that need to be
considered to derive a complete relation to the Collins scheme.

####  3.2.1 Absence of mixing

Our derivation of the factorization formula has not specified the quark
flavor, and the result actually implies that mixings between quarks of
different flavors or quarks and gluons do not exist. This lack of mixing is a
generic feature for quasi-TMDs of all parton species.

The LR and Collins schemes differ only in the order of their $\epsilon\to 0$
and $y_{B}\to-\infty$ limits. When matching the schemes, mixing between quark
and gluon channels could only occur if $y_{B}\to-\infty$ leads to a UV-
divergent counterterm that contracts the flavor indices of the quark fields in
the nonlocal bilinear operator. This cannot happen because these quark fields
are always spacelike separated, and thus their exchange of intermediate
particles in a Feynman diagram is exponentially suppressed and does not
generate a new UV divergence. As a result, the $y_{B}\to-\infty$ limit can
only change UV divergences locally at quark-Wilson line vertices and in Wilson
line wavefunction renormalization, which leaves the parton flavor intact.

In contrast to the quasi-TMDs, the quark quasi-PDFs are defined from bilinear
operators with a straight Wilson line along the $z$ direction. In the infinite
boost limit, the spacelike separated quarks will approach the light-cone, thus
inducing a nonlocal UV divergence that contracts the quark flavor indices and
allows mixing with the gluon quasi-PDF [54].

There is another perspective from which we can understand the lack of mixing.
The LR and JMY schemes are related to each other by analytic continuation
between space- and timelike Wilson lines. Thus, the JMY scheme should
factorize similarly to eq. (3.26), as we check at one-loop order in appendix
A.1. If there were quark-gluon or flavor mixing in the Collins-to-JMY
matching, then such mixing would manifest in the TMD factorization formula for
the Drell-Yan or SIDIS cross-section in either scheme; but it does not.
Therefore, no mixing should occur in the Collins-to-JMY or Collins-to-LR
matching.

Note that this factorization relation holds for quark and gluon quasi-TMDs
with all spin-dependent structures [40, 49], so we can use it to compute the
ratios of spin-dependent TMDs from the quasi-TMDs or quasi-beam functions, an
approach that has been proposed and used in refs. [27, 28, 29, 30, 31, 32]. In
summary, all quark and gluon quasi-TMDs should satisfy the factorization
relation in eq. (1.1). We cross-check our all-orders analysis above by
explicit one-loop order calculations in appendix A.2.

####  3.2.2 Resummed result for the matching coefficient

From eq. (2.3.1) we can write $\ln\tilde{f}_{q_{i}/h}=\ln\lim_{\tilde{\eta}\gg
b_{T}}\tilde{B}^{[\tilde{\Gamma}]}_{q_{i}/h}-\frac{1}{2}\ln\lim_{\tilde{\eta}\gg
b_{T}}\tilde{S}^{q}$, and then using eq. (3.27) in the form
$\ln\tilde{f}_{q_{i}/h}=\ln C_{q}+\ln f_{q_{i}/h}$ we can derive the momentum
evolution equation of the quasi-beam function [33]:

$\displaystyle{\mathrm{d}\over\mathrm{d}\ln(2x\tilde{P}^{z})}\>\ln\lim_{\tilde{\eta}\gg
b_{T}}\tilde{B}_{q_{i}/h}^{[\tilde{\Gamma}]}(x,\vec{b}_{T},\mu,\tilde{\eta},x\tilde{P}^{z})$
$\displaystyle=\gamma_{\zeta}^{q}(b_{T},\mu)+\gamma_{C}^{q}(2x\tilde{P}^{z},\mu)\,.$
(3.28)

Here the limit $\tilde{\eta}\gg b_{T}$ should be understood as expanding in
large $\tilde{\eta}$, where the quasi-beam function has a divergent dependence
on $\tilde{\eta}$ which is however independent of $x\tilde{P}^{z}$, and hence
drops out. In taking the $\ln(2x\tilde{P}^{z})$ derivative we hold
$y_{B}-y_{n}$ fixed, so there is no contribution from the $\tilde{S}^{R}$
term. We also used the large momentum formula
$\tilde{\zeta}=(2x\tilde{P}^{z}e^{y_{B}-y_{n}})^{2}$ to convert the
$\ln\tilde{\zeta}$ derivative of $f_{q_{i}/h}$ to give
$\mathrm{d}f_{q_{i}/h}/\mathrm{d}\ln(2x\tilde{P}^{z})=\gamma_{\zeta}^{q}(b_{T},\mu)$.
The other anomalous dimension appearing in eq. (3.28) is

$\displaystyle\gamma_{C}^{q}(2x\tilde{P}^{z},\mu)$
$\displaystyle={\mathrm{d}\over\mathrm{d}\ln(2x\tilde{P}^{z})}\ln
C_{q}(x\tilde{P}^{z},\mu)\,.$ (3.29)

Since the quasi-beam function $\tilde{B}_{q_{i}/h}^{[\tilde{\Gamma}]}$ has a
local UV counterterm $Z_{\rm uv}^{B}(\mu,\epsilon)$ according to the auxiliary
field formalism, the sum of anomalous dimensions in eq. (3.28) must be $\mu$
independent. The known perturbative structure of the CS kernel
$\gamma_{\zeta}^{q}(b_{T},\mu)$ then implies that

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\ln\mu}\gamma_{C}^{q}(2x\tilde{P}^{z},\mu)=-\frac{\mathrm{d}}{\mathrm{d}\ln\mu}\gamma_{\zeta}^{q}(b_{T},\mu)=2\Gamma_{\rm
cusp}^{q}[\alpha_{s}(\mu)]\,.$ (3.30)

It follows that $\gamma_{C}^{q}$ can be written to all orders as

$\displaystyle\gamma_{C}^{q}(2x\tilde{P}^{z},\mu)=2\int_{2x\tilde{P}^{z}}^{\mu}\\!\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\>\Gamma_{\rm
cusp}^{q}\big{[}\alpha_{s}(\mu^{\prime})\big{]}+\gamma_{C}^{q}\big{[}\alpha_{s}(2x\tilde{P}^{z})\big{]}\,.$
(3.31)

Here, $\Gamma_{\rm cusp}^{q}$ and $\gamma_{C}^{q}$ are the cusp and noncusp
anomalous dimensions, whose series expansions are given by

$\displaystyle\Gamma_{\rm cusp}^{q}[\alpha_{s}]$
$\displaystyle=\sum_{n=0}^{\infty}\Big{(}\frac{\alpha_{s}}{4\pi}\Big{)}^{n+1}\,\Gamma_{n}^{q}\,,\qquad\qquad\gamma_{C}^{q}[\alpha_{s}]=\sum_{n=0}^{\infty}\Big{(}\frac{\alpha_{s}}{4\pi}\Big{)}^{n+1}\,\gamma_{C\,n}^{q}\,,$
(3.32) $\displaystyle\text{with}\ \ \
\Gamma_{0}^{q}=4C_{F}\,,\qquad\Gamma_{1}^{q}=4C_{F}\Big{[}\Big{(}\frac{67}{9}-\frac{\pi^{2}}{3}\Big{)}C_{A}-\frac{20}{9}T_{F}n_{f}\Big{]}\,,\>\ldots\,.$

Here $n_{f}$ is the number of light quark flavors and for QCD $C_{F}=4/3$,
$C_{A}=3$, and $T_{F}=1/2$. We also expand the QCD $\beta$-function,
$\mathrm{d}\alpha_{s}(\mu)/\mathrm{d}\ln\mu=\beta[\alpha_{s}(\mu)]$ as

$\displaystyle\beta[\alpha_{s}]$
$\displaystyle=-2\alpha_{s}\sum_{n=0}^{\infty}(\alpha_{s}/4\pi)^{n+1}\beta_{n}\,,$
(3.33) $\displaystyle\text{with}\ \ \
\beta_{0}=\frac{11}{3}C_{A}-\frac{4}{3}T_{F}n_{f}\,,\quad\beta_{1}=\frac{34}{3}C_{A}^{2}-\Big{(}\frac{20}{3}C_{A}+4C_{F}\Big{)}T_{F}n_{f}\,,\>\ldots\,.$

Solving eq. (3.29) gives the general all-orders resummed result

$\displaystyle C_{q}(x\tilde{P}^{z},\mu)$
$\displaystyle=C_{q}\left[\alpha_{s}(\mu)\right]\>\exp\left[\int_{\mu}^{2x\tilde{P}^{z}}\\!\\!\frac{\mathrm{d}\tau}{\tau}\
\gamma_{C}^{q}(\tau,\mu)\right]$ (3.34)
$\displaystyle=C_{q}\left[\alpha_{s}(\mu)\right]\exp\left[\int^{\alpha_{s}(2x\tilde{P}^{z})}_{\alpha_{s}(\mu)}\\!\frac{\mathrm{d}\alpha}{\beta[\alpha]}\bigg{(}\int_{\alpha}^{\alpha_{s}(\mu)}\frac{\mathrm{d}\alpha^{\prime}}{\beta[\alpha^{\prime}]}2\Gamma_{\rm
cusp}^{q}[\alpha^{\prime}]+\gamma_{C}^{q}[\alpha]\bigg{)}\right]\,,$

where the boundary condition is given by
$C_{q}\left[\alpha_{s}(\mu)\right]=C_{q}(\mu/2,\mu)$. Explicit results at a
given order can be obtained by substituting fixed order series for
$\Gamma_{\rm cusp}[\alpha^{\prime}]$, $\gamma_{C}^{q}[\alpha]$, and
$C_{q}[\alpha_{s}]$.

Using the known one-loop results [35, 33] we have

$\displaystyle\gamma_{C}^{q}(2x\tilde{P}^{z},\mu)$
$\displaystyle=\frac{\alpha_{s}(\mu)C_{F}}{\pi}\bigg{(}-\ln{(2x\tilde{P}^{z})^{2}\over\mu^{2}}+1\bigg{)}+{\cal
O}(\alpha_{s}^{2})\,,$ (3.35) $\displaystyle C_{q}[\alpha_{s}]$
$\displaystyle=1+\frac{\alpha_{s}C_{F}}{2\pi}\left(-2+{\pi^{2}\over
12}\right)+{\cal O}(\alpha_{s}^{2})\,,$

which is consistent with $\Gamma_{0}^{q}=4C_{F}$ and allows us to identify
$\gamma_{C\,0}^{q}=4C_{F}$. To obtain results for eq. (3.34) at next-to-
leading-logarithmic (NLL) order for the double logarithmic series present in
$C_{q}$, we can utilize $C_{q}[\alpha_{s}(\mu)]=1$ together with the two-loop
cusp anomalous dimension, and one-loop regular anomalous dimension. Using the
notation of Ref. [76] for evolution kernels, the matching coefficient at NLL
is then

$\displaystyle C_{q}\big{(}x\tilde{P}^{z},\mu\big{)}^{\rm NLL}$
$\displaystyle=\exp\Big{[}-2K_{\Gamma}^{q}\big{(}2x\tilde{P}^{z},\mu\big{)}-K_{\gamma}^{q}\big{(}2x\tilde{P}^{z},\mu\big{)}\Big{]}\,,$
(3.36) $\displaystyle K_{\Gamma}^{q}(\mu_{0},\mu)$
$\displaystyle=-\frac{\Gamma_{0}^{q}}{4\beta_{0}^{2}}\,\biggl{\\{}\frac{4\pi}{\alpha_{s}(\mu_{0})}\,\Bigl{(}1-\frac{1}{r}-\ln
r\Bigr{)}+\biggl{(}\frac{\Gamma_{1}^{q}}{\Gamma_{0}^{q}}-\frac{\beta_{1}}{\beta_{0}}\biggr{)}(1-r+\ln
r)+\frac{\beta_{1}}{2\beta_{0}}\ln^{2}r\biggr{\\}},$ $\displaystyle
K_{\gamma}^{q}(\mu_{0},\mu)$
$\displaystyle=-\frac{\gamma_{C0}^{q}}{2\beta_{0}}\,\ln r\,,$

where $r=\alpha_{s}(\mu)/\alpha_{s}(\mu_{0})$. Expanding we find agreement
with an earlier ${\cal O}(\alpha_{s}^{2})$ analysis for the terms we can
predict at NLL, given by all the ${\cal
O}(\alpha_{s}^{2}\ln^{j}\frac{2x\tilde{P}^{z}}{\mu})$ terms with $j=2,3,4$ in
Eqs.(25,26) of Ref. [39]. Equation (3.36) can be expanded to higher orders in
$\alpha_{s}$, and then predicts the terms in $\ln C_{q}(x\tilde{P}^{z},\mu)$
of the form $\alpha_{s}^{j}\ln^{j+1}\frac{2x\tilde{P}^{z}}{\mu}$ and
$\alpha_{s}^{j}\ln^{j}\frac{2x\tilde{P}^{z}}{\mu}$ for any $j$.

Results for $C_{q}(x\tilde{P}^{z},\mu)$ beyond NLL can be obtained from eq.
(3.34) by substituting in higher order results for the anomalous dimensions
and boundary condition. (Results for $K_{\Gamma}^{q}$ and $K_{\gamma}^{q}$ in
terms of anomalous dimensions can be found in many places in the literature to
order N3LL, see also ref. [77] for an exact solution.) An RGE equation in the
form in eq. (3.28) will also hold for the $x\tilde{P}^{z}$ anomalous dimension
for the gluon TMD, so a resummed formula for its matching coefficient
$C_{g}(x\tilde{P}^{z},\mu)$ is given by the above expressions with $q\to g$
and replacement by the gluon cusp and non-cusp anomalous dimensions.

####  3.2.3 Ratios of quasi-TMDs

The lack of mixing in the factorization formula eq. (1.1) for quasi-TMDs
allows us to calculate ratios of TMDs of all flavor and spin structures more
easily since there are cancellations between the numerator and denominators.
This approach of studying ratios was pioneered in the Lorentz-invariant method
of Refs. [27, 28] using the MHENS scheme. This has been shown to have great
utility for exploring ratios involving an integral over $x$ and different spin
and flavor choices [28, 29, 30, 31, 32]. We return to discuss the prospects
for including renormalization and matching corrections in the MHENS scheme
approach in section 3.2.4.

For quasi-TMDs the ability to more easily calculate ratios of spin dependent
structure functions was observed for quark non-singlet distributions in Refs.
[40, 41], and occur due to the universality of the quasi-TMD to Collins-TMD
matching coefficient. Our result in eq. (1.1) enable us to extend these
observations to all orders in $\alpha_{s}$, and include singlet quark
distributions and gluon distributions. Since the quasi soft factor
$\tilde{\Delta}^{R}$ in eq. (2.3.1) and the matching coefficients $C_{q,g}$ in
eq. (1.1) only depend on the color representation, we can formulate ratios of
quark or gluon quasi-TMDs where these components cancel, and thus immediately
can be related to the analogous ratios for the quark and gluon TMDs in the
Collins scheme. In particular we have

$\displaystyle\frac{\tilde{f}_{q_{i}/h}^{[\tilde{\Gamma}_{1}]}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z})}{\tilde{f}_{q_{j}/h^{\prime}}^{[\tilde{\Gamma}_{2}]}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z})}$
$\displaystyle=\lim_{\tilde{\eta}\to\infty}\frac{\tilde{B}_{q_{i}/h}^{[\tilde{\Gamma}_{1}]}(x,\vec{b}_{T},\mu,\tilde{\eta},x\tilde{P}^{z})}{\tilde{B}_{q_{j}/h^{\prime}}^{[\tilde{\Gamma}_{2}]}(x,\vec{b}_{T},\mu,\tilde{\eta},x\tilde{P}^{z})}=\frac{f_{q_{i}/h}^{[\Gamma_{1}]}(x,\vec{b}_{T},\mu,\zeta)}{f_{q_{j}/h^{\prime}}^{[\Gamma_{2}]}(x,\vec{b}_{T},\mu,\zeta)}\,,$
(3.37)
$\displaystyle\frac{\tilde{f}_{g/h}^{[\tilde{\mu}_{1}\tilde{\nu}_{1}\tilde{\rho}_{1}\tilde{\sigma}_{1}]}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z})}{\tilde{f}_{g/h^{\prime}}^{[\tilde{\mu}_{2}\tilde{\nu}_{2}\tilde{\rho}_{2}\tilde{\sigma}_{2}]}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z})}$
$\displaystyle=\lim_{\tilde{\eta}\to\infty}\frac{\tilde{B}_{g/h}^{[\tilde{\mu}_{1}\tilde{\nu}_{1}\tilde{\rho}_{1}\tilde{\sigma}_{1}]}(x,\vec{b}_{T},\mu,\tilde{\eta},x\tilde{P}^{z})}{\tilde{B}_{g/h^{\prime}}^{[\tilde{\mu}_{2}\tilde{\nu}_{2}\tilde{\rho}_{2}\tilde{\sigma}_{2}]}(x,\vec{b}_{T},\mu,\tilde{\eta},x\tilde{P}^{z})}=\frac{f_{g/h}^{[\mu_{1}\nu_{1}\rho_{1}\sigma_{1}]}(x,\vec{b}_{T},\mu,\zeta)}{f_{g/h^{\prime}}^{[\mu_{2}\nu_{2}\rho_{2}\sigma_{2}]}(x,\vec{b}_{T},\mu,\zeta)}\,.$

Here $q_{i}$ and $q_{j}$ can be different quark flavors, $h$ and $h^{\prime}$
can be different hadrons, and the superscripts can be different spin
structures with Dirac matrices $\Gamma^{1},\,\Gamma^{2}$ for quark (quasi-)
TMDs and Lorentz indices $\mu_{k},\nu_{k},\rho_{k},\sigma_{k}$ with $k=1,2$
for gluon (quasi-)TMDs.

To calculate the ratios in eq. (3.37) as a function of $x$, one must first
compute the matrix elements for the quasi-beam functions at all $b^{z}$, then
take the Fourier transform. Because UV divergences in the bare quasi-beam
function matrix elements are $b^{z}$-independent, they factor out of the
Fourier integral. So, in principle we can skip renormalization and matching to
the ${\overline{\mathrm{MS}}}$ scheme when calculating TMD ratios, _if_ there
are no $b^{z}$-dependent finite operator mixings on the discretized lattice.
However, in the presence of such mixings, lattice renormalization is
necessary, as studied in Refs. [43, 73]. Also, in numerical analyses it can be
advantageous to consider the $\tilde{\eta}\to\infty$ limit separately for the
numerator and denominator of eq. (3.37) separately. This can be accomplished
by utilizing the naive quasi-soft function or quasi-beam function at $b^{z}=0$
to cancel the large $\tilde{\eta}$-dependence.

####  3.2.4 Matching MHENS and continuum TMDs

We now consider the relation between the MHENS lattice TMD and Collins
continuum TMD, focusing again on the quark case. In the literature, the MHENS
scheme has primarily been used to study matrix elements evaluated at $P\cdot
b=0$ [27, 28, 29, 30, 31, 32]. In this case, the equal-time-restricted Wilson
line path in the MHENS beam function is the same as that of the quasi-beam
function. This is easily seen by comparing the integral over $x$ of the MHENS
beam function in eq. (2.39), with the integral over $x$ of the quasi-beam
function in eq. (2.3.1), and noting that both give the same correlator
$\Omega^{[\Gamma]}_{q/h}(\vec{b}_{T},\tilde{P},a,\tilde{\eta}\hat{z},0)=\tilde{\Phi}_{\rm
unsubtr.}^{[\Gamma]}(\vec{b}_{T},\tilde{P},a,\tilde{\eta}\hat{z})$ times a
factor of $N_{\Gamma}/P^{z}$. For the integral over $x$ we define

$\displaystyle\int\\!\mathrm{d}x\
\tilde{f}_{q_{i}/h}^{[\Gamma]}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z},\tilde{\eta})$
$\displaystyle=\tilde{f}_{q_{i}/h}^{[\Gamma]}(b^{z}=0,\vec{b}_{T},\mu,\tilde{P}^{z},y_{n}-y_{B},\tilde{\eta})$
(3.38) $\displaystyle=f_{q_{i}/h}^{[\Gamma]{\rm
MHENS}}(b^{z}=0,\vec{b}_{T},\mu,\tilde{P}^{z},y_{n}-y_{B},\tilde{\eta})\,.$

The first quasi-TMD here has $x$-dependence in three of its arguments (two
written explicitly and the other in $\tilde{\zeta}$), so it is convenient to
write the $x$-independent result as a new function, whose distinction is
tagged by the first $b^{z}=0$ argument. We adopt the same notation for the
MHENS TMD, as shown. Given this correspondence, we can simply adopt the same
terms used in defining the quasi-TMD in eq. (2.3.1) to define a renormalized
and soft subtracted MHENS TMD for $\tilde{P}\cdot b=0$, $b^{0}=0$, and
$v=\hat{z}$, giving

$\displaystyle f_{q_{i}/h}^{[\Gamma]{\rm
MHENS}}(b^{z}=0,\vec{b}_{T},\mu,\tilde{P},y_{n}-y_{B},\tilde{\eta})$ (3.39)
$\displaystyle\quad\equiv\lim_{a\to 0}Z^{\prime}_{\rm
uv}(\mu,\tilde{\mu})\,Z_{\rm
uv}(a,\tilde{\mu},y_{n}-y_{B})\,\frac{N_{\Gamma}}{P^{z}}\,\tilde{\Phi}_{\rm
unsubtr.}^{[\Gamma]q_{i}/h}(\vec{b}_{T},\tilde{P},a,\tilde{\eta}\hat{z})\,\tilde{\Delta}_{S}^{q}(b_{T},a,\tilde{\eta},y_{n},y_{B})\,.$

The limit $\tilde{\eta}\to\infty$ of eq. (3.39) gives a finite result
independent of $\tilde{\eta}$, since the Wilson line self-energy power law
divergences cancel between $\tilde{\Phi}^{[\Gamma]}_{\rm unsubtr.}$ and
$\tilde{\Delta}_{S}^{q}$. With this definition for the MHENS TMD, our result
in eq. (1.1) relating the quasi- and Collins TMDs also yields a relationship
between the MHENS and Collins TMDs:

$\displaystyle\lim_{\tilde{\eta}\to\infty}\tilde{f}_{q_{i}/h}^{[\Gamma]{\rm
MHENS}}(b^{z}=0,\vec{b}_{T},\mu,\tilde{P},y_{n}-y_{B},\tilde{\eta})$
$\displaystyle=\int\\!\mathrm{d}x\
C_{q}(x\tilde{P}^{z},\mu)\>f_{q_{i}/h}^{[\Gamma]}(x,\vec{b}_{T},\mu,\zeta)$
$\displaystyle\qquad\times\bigg{\\{}1+\mathcal{O}\biggl{[}{1\over(\tilde{P}^{z}b_{T})^{2}},{\Lambda_{\rm
QCD}^{2}\over(\tilde{P}^{z})^{2}}\biggr{]}\bigg{\\}}.$ (3.40)

Here, the MHENS TMD (or quasi-TMD) on the LHS is calculated with states
involving proton momentum $\tilde{P}$, while the Collins TMD on the RHS
utilizes states with proton momentum $P$. We have chosen to relate these two
momenta by the rapidity relation $y_{P}=y_{\tilde{P}}+y_{B}$ that appeared in
our proof of factorization. In the large $\tilde{P}^{z}$ limit, we have
$\zeta=(2xP^{z}e^{-y_{n}})^{2}=(2x\tilde{P}^{z}e^{y_{B}-y_{n}})^{2}=\tilde{\zeta}$,
which eliminates the $\ln(\tilde{\zeta}/\zeta)$-dependent term in eq. (1.1),
which would otherwise appear in the integrand in eq. (3.2.4). Note that for
large $\tilde{P}^{z}$ we also have $\hat{\zeta}\to\infty$, where $\hat{\zeta}$
was defined in eq. (2.43).

We now consider the implications of the factorization in eq. (3.2.4) for
computing Collins TMD ratios. Taking ratios of MHENS TMDs with different
choices of spin structures $\Gamma$, we see that the UV renormalization
factors $Z^{\prime}_{\rm uv}$ and $Z_{\rm uv}$ and soft factor
$\tilde{\Delta}_{S}^{q}$ all drop out. Thus, ratios of MHENS beam functions
give us information about the ratios of Collins-TMDs,

$\displaystyle\lim_{\stackrel{{\scriptstyle\text{\scriptsize$a\to
0$}}}{{\tilde{\eta}\to\infty}}}\ \frac{N_{\Gamma_{1}}\tilde{\Phi}_{\rm
unsubtr.}^{[\Gamma_{1}]q_{i}/h}(\vec{b}_{T},\tilde{P},a,\tilde{\eta}\hat{z})}{N_{\Gamma_{2}}\tilde{\Phi}_{\rm
unsubtr.}^{[\Gamma_{2}]q_{i}/h}(\vec{b}_{T},\tilde{P},a,\tilde{\eta}\hat{z})}=\frac{\int\mathrm{d}x\>C_{q}(x\tilde{P}^{z},\mu)f_{q_{i}/h}^{[\Gamma_{1}]}(x,\vec{b}_{T},\mu,\zeta)}{\int\mathrm{d}x\>C_{q}(x\tilde{P}^{z},\mu)f_{q_{i}/h}^{[\Gamma_{2}]}(x,\vec{b}_{T},\mu,\zeta)}\,,$
(3.41)

albeit with an integration over the matching coefficient
$C_{q}(x\tilde{P}^{z},\mu)$. Here there is power law sensitivity to
$\tilde{\eta}$ in the numerator and denominator of the LHS, but this
sensitivity cancels in the ratio, as does the dependence on $a$ (assuming that
there is no mixing amongst spin structures for the lattice fermion
discretization chosen [78, 43, 73, 79]). As explained in section 3.1.3,
$C_{q}$ arises precisely because of the different orders in which
renormalization and the large rapidity limit are taken to be performed between
the MHENS and Collins TMDs. Equation (3.41) reduces to the relation between
MHENS beam functions and the moment of the Collins TMDs discussed in earlier
literature [28, 29, 30, 31, 32] if one sets $C_{q}=1+{\cal O}(\alpha_{s})$,
i.e. works to tree level in the matching coefficient. Beyond tree-level, the
convolution becomes nontrivial. Nevertheless, one can plug the TMDs from
global analysis into eq. (3.41) to compare with the lattice ratio of the MHENS
beam function. It is worth noting that the CS evolution for the $\zeta$
dependence is multiplicative and independent of $x$, and hence the ratio of
Collins TMDs on the RHS is independent of $\zeta$ as long as the same value of
$\zeta$ is used in the numerator and denominator.

By constructing lattice TMD ratios with the same spin structures but different
momenta $\tilde{P}^{z}$ in the numerator and denominator, one can extract the
CS kernel $\gamma_{\zeta}(\mu,b_{T})$ [35]. Once again one cannot avoid the
need to include the matching coefficient whether working in longitudinal
momentum or position space [37]. Thus a formula for obtaining $\gamma_{\zeta}$
that does not explicitly rely on a truncation of the $\alpha_{s}$ expansion
always requires calculation of the full $b^{z}$ dependence of lattice beam
functions. When the series expansion in $\alpha_{s}$ is utilized, a broader
range of extraction techniques are possible, see for example Ref. [45].

Next we consider the use of the MHENS correlator for obtaining $x$-dependent
information about Collins TMDs. There are two complications in this case
relative to the use of quasi-TMDs. The first is that the MHENS staple-shaped
Wilson line path for the beam functions has non-trivial cusp angles
$\gamma(v,b)$, which from eq. (2.6) is given by

$\displaystyle\cosh[\gamma(v,b)]=\pm\frac{v\cdot b}{|v||b|}\,,$ (3.42)

where the sign is determined by $\pm={\rm sign}(\tilde{\eta})$. For $b^{z}\neq
0$ we have $v\cdot b\neq 0$, and the UV renormalization factor $Z_{\rm
uv}^{\rm MHENS}$ will depend on $\gamma(v,b)$ irrespective of whether or not
we extrapolate towards infinitely long staples, $|\tilde{\eta}|\to\infty$.
This complicates the analysis because the UV renormalization is now $b^{z}$
dependent and hence will not cancel in ratios formed from correlators with the
same $x$.

The second complication for $x$-dependent MHENS calculations is that the
length of the Wilson line path becomes $b^{z}$-dependent, since using eq.
(2.5) we have

$\displaystyle L_{\rm staple}^{\rm MHENS}$
$\displaystyle=2|\tilde{\eta}v|+|b|\,.$ (3.43)

This implies that if we do a Fourier transform in $b\cdot P$ to obtain
$x$-dependent correlators, then the power law dependence on the staple length
does not cancel for ratios taken at finite $\tilde{\eta}$, where the $|b|$
term can not be neglected.

To make these issues more transparent, we consider the generalization of the
MHENS TMD definition in eq. (3.39) that is needed to include
$b^{z}$-dependence. Working in an equal-time configuration with $b^{0}=0$ and
$v^{0}=0$ that is suitable for lattice calculations, we expect that the
definition would take the form

$\displaystyle f_{q_{i}/h}^{[\Gamma]{\rm MHENS}}(b,\mu,\tilde{P},v,\tilde{v})$
$\displaystyle\equiv\lim_{\stackrel{{\scriptstyle\text{\scriptsize$a\to
0$}}}{{|\tilde{\eta}|\to\infty}}}Z^{\prime}_{\rm uv}(\mu,\tilde{\mu})\,Z_{\rm
uv}^{\rm MHENS}\big{[}a,\tilde{\mu},v,\tilde{v},\gamma(v,b),\ldots\big{]}\,$
(3.44) $\displaystyle\qquad\times\frac{N_{\Gamma}}{P^{z}}\,\tilde{\Phi}_{\rm
unsubtr.}^{[\Gamma]q_{i}/h}(b,\tilde{P},a,\tilde{\eta}v)\,\tilde{\Delta}_{S}^{q\,{\rm
MHENS}}(b,a,\tilde{\eta}v,\tilde{\eta}\tilde{v})\,.$

Although $Z_{\rm uv}^{\rm MHENS}$ depends on the cusp angle $\gamma(v,b)$ this
is unlikely to be a fundamental road block, since this renormalization can be
carried out perturbatively (for example, four-loop results for the related
cusp-anomalous dimension are now available in ${\overline{\mathrm{MS}}}$
[80]). It also seems likely that a non-perturbative method of carrying out the
calculation of $Z_{\rm uv}^{\rm MHENS}$ on the lattice could also be
formulated. A greater difficulty will be determining a suitable soft factor
$\tilde{\Delta}_{S}^{q\,{\rm MHENS}}=(S^{q}_{\rm MHENS})^{-1/2}$, which itself
must satisfy three non-trivial constraints. In particular, it should be
constructed from a soft function $S^{q}_{\rm MHENS}$ which is a vacuum matrix
element of a closed Wilson loop that has a total length $2L_{\rm staple}^{\rm
MHENS}$ with eq. (3.43). This ensures that the $|\tilde{\eta}|\to\infty$ limit
for $f_{q_{i}/h}^{[\Gamma]{\rm MHENS}}$ will exist. An additional constraint
is that it should include dependence which compensates for the mismatch in the
Lorentz invariants in the 12th and 13th rows of table 2. Finally, it must have
the proper infrared dependence on $b_{T}$ such that
$\tilde{f}_{q_{i}/h}^{[\Gamma]{\rm MHENS}}$ correctly reproduces the infrared
structure of the Collins TMD. This last constraint is necessary for a
factorization formula relating the MHENS TMD and Collins TMD to exist. The
construction of a suitable $\tilde{\Delta}_{S}^{q\,{\rm MHENS}}$ involves two
steps: finding an operator definition for this quantity satisfying the above
constraints, and then developing a method by which this factor can be computed
with lattice QCD. The argument $\tilde{v}$ that we have written for
$\tilde{\Delta}_{S}^{q\,{\rm MHENS}}$ in eq. (3.44) should be a suitably
chosen space-like vector. The ellipses in $Z_{\rm uv}^{\rm MHENS}$ in eq.
(3.44) denote any further arguments needed due to UV renormalization for the
soft factor (like additional cusp angles).

Let us assume that a suitable $\tilde{\Delta}_{S}^{q\,{\rm MHENS}}$ has been
determined. In this case the factorization formula for the MHENS TMD to
Collins TMD should take the form

$\displaystyle f_{q_{i}/h}^{[\Gamma]{\rm
MHENS}}\big{(}x,\vec{b}_{T},\mu,\tilde{P},v,\tilde{v}\big{)}$
$\displaystyle\equiv\int\\!\frac{\mathrm{d}(b\cdot\tilde{P})}{2\pi}\>e^{-\mathrm{i}x(\tilde{P}\cdot
b)}\>f_{q_{i}/h}^{[\Gamma]{\rm MHENS}}(b,\tilde{P},\mu,v,\tilde{v})$
$\displaystyle=C_{q}^{\rm
MHENS}(xv\cdot\\!\tilde{P},\ldots,\mu)\,\exp\biggl{[}\frac{1}{2}\gamma_{\zeta}^{q}(\mu,b_{T})\ln\frac{\tilde{\zeta}}{\zeta}\,\bigg{]}\,f_{q_{i}}^{[\Gamma]}(x,\vec{b}_{T},\mu,\zeta)$
$\displaystyle\qquad\times\bigg{\\{}1+\mathcal{O}\biggl{[}{1\over\big{(}xv\cdot\\!\tilde{P}b_{T}\big{)}^{2}},{\Lambda_{\rm
QCD}^{2}\over\big{(}xv\cdot\\!\tilde{P}\big{)}^{2}}\biggr{]}\bigg{\\}}\,,$
(3.45)

where based on results from the JMY scheme we expect
$\tilde{\zeta}=2(x\tilde{P}\cdot
v)^{2}\sqrt{|v^{2}||\tilde{v}^{2}|}/(v^{2}\,v\cdot\tilde{v})$. Our notation
anticipates the fact that the matching coefficient $C_{q}^{\rm MHENS}$ is
likely to differ from the $C_{q}$ in the quasi-TMD factorization. This is
expected due to the fact that the UV behavior of $\tilde{\Delta}_{S}^{q\,{\rm
MHENS}}$ can differ from $\tilde{\Delta}_{S}^{q}$, and the fact that $Z_{\rm
uv}^{\rm MHENS}$ differs from $Z_{\rm uv}$.

Instead of directly trying to match the MHENS-TMD onto the Collins-TMD as in
eq. (3.2.4), one can consider determining the $x$-dependence of ratios. For
example, taking ratios with potentially different Dirac structures, flavors
$q_{i}$ and $q_{j}$, and hadrons $h$ and $h^{\prime}$, but the same $x$ and
$\tilde{P}$ we have

$\displaystyle\frac{f_{q_{i}/h}^{[\Gamma_{1}]}(x,\vec{b}_{T},\mu,\zeta)}{f_{q_{j}/h^{\prime}}^{[\Gamma_{2}]}(x,\vec{b}_{T},\mu,\zeta)}=\frac{f_{q_{i}/h}^{[\Gamma_{1}]{\rm
MHENS}}(x,\vec{b}_{T},\mu,\tilde{P},v,\tilde{v})}{f_{q_{j}/h^{\prime}}^{[\Gamma_{2}]{\rm
MHENS}}(x,\vec{b}_{T},\mu,\tilde{P},v,\tilde{v})}$ (3.46)
$\displaystyle=\lim_{\stackrel{{\scriptstyle\text{\scriptsize$a\to
0$}}}{{|\tilde{\eta}|\to\infty}}}\frac{\int\mathrm{d}(b\cdot\tilde{P})\,e^{-\mathrm{i}x(\tilde{P}\cdot
b)}\,Z_{\rm uv}^{\rm MHENS}(a,\ldots)N_{\Gamma_{1}}\tilde{\Phi}_{\rm
unsubtr.}^{[\Gamma_{1}]q_{i}/h}(b,\tilde{P},a,\tilde{\eta}\hat{z})\,\tilde{\Delta}_{S}^{q\,{\rm
MHENS}}(b,a,\tilde{\eta}v,\tilde{\eta}\tilde{v})}{\int\mathrm{d}(b\cdot\tilde{P})\,e^{-\mathrm{i}x(\tilde{P}\cdot
b)}\,Z_{\rm uv}^{\rm MHENS}(a,\ldots)N_{\Gamma_{2}}\tilde{\Phi}_{\rm
unsubtr.}^{[\Gamma_{2}]q_{j}/h^{\prime}}(b,\tilde{P},a,\tilde{\eta}\hat{z})\,\tilde{\Delta}_{S}^{q\,{\rm
MHENS}}(b,a,\tilde{\eta}v,\tilde{\eta}\tilde{v})}.$

In taking ratios using eq. (3.2.4) the $C_{q}^{\rm MHENS}$ and CS kernel terms
drop out, giving the first equality in eq. (3.46). Using eq. (3.44) then gives
the second equality. For finite $\tilde{\eta}$ the soft factors
$\tilde{\Delta}_{S}^{q\,{\rm MHENS}}(b,a,\tilde{\eta}v,\tilde{\eta}\tilde{v})$
do not cancel out from the numerator and denominator, due to their dependence
on the integration variable, namely the component of $b$ that is parallel to
$\tilde{P}$.666It is possible that one may be able to construct
$\tilde{\Delta}_{S}^{q\,{\rm MHENS}}$ such that the dependence on the
component $b\cdot\tilde{P}/m_{h}$ is subleading as $|\tilde{\eta}|\to\infty$,
just like it is in $L_{\rm staple}^{\rm MHENS}$, in which case these soft
factors can be canceled out in the ratio in eq. (3.46) as
$|\tilde{\eta}|\to\infty$. In addition the UV counterterms $Z_{\rm uv}^{\rm
MHENS}(a,\ldots,\gamma(b,v),\ldots)$ depend on the integration variable
through their dependence on the cusp angle $\gamma(b,v)$, and the dependence
on these variables remains regardless of how large $\tilde{\eta}$ is.777It is
possible that one may be able to define the soft factor $\Delta_{S}^{q\,{\rm
MHENS}}$ so as to cancel the UV dependence on the cusp angles, and thus make
$Z_{\rm uv}^{\rm MHENS}$ independent of $\gamma(b,v)$. However, this makes it
more probable that $\Delta_{S}^{q\,{\rm MHENS}}$ will not cancel between the
numerator and denominator of eq. (3.46) as $|\tilde{\eta}|\to\infty$. Equation
(3.46) can be contrasted with the ratios involving quasi-TMDs in the first
line of eq. (3.37). The situation is simpler for the quasi-TMDs because of the
simpler dependence of the UV counterterm and soft factor, which cancel out in
the ratios at finite $\tilde{\eta}$.

##  4 Conclusion

The central focus of this paper was to derive the factorization formula eq.
(1.1) that relates the lattice-calculable quasi-TMD to physical TMD schemes
through a simple perturbative matching coefficient. This formula is valid at
all orders in $\alpha_{s}$ up to power corrections for any light quark flavor
and for gluons.

We began our derivation by developing a generalized TMD notational framework
applicable to both continuum and lattice TMDs, enabling us to unify various
choices used in the literature and thus more easily unpack the relationships
between various physical and lattice TMD schemes. Comparing the operator
structures and Lorentz invariants appearing in each scheme, we observed a
close relation between the Collins and quasi-TMDs. We then constructed a new
continuum TMD scheme intermediate between the Collins and quasi-TMDs, which we
called the large rapidity (LR) scheme. The LR and Collins schemes differ by
the order of UV renormalization and the lightcone limits, so they are related
by a perturbative matching in the spirit of LaMET. Meanwhile, using Lorentz
invariance we showed that the quasi- and LR TMDs are equivalent. This enabled
us to prove the full factorization relation between the Collins and quasi-
TMDs. For any quark flavor $q_{i}$ and for gluons $g$, the relations are

$\displaystyle\tilde{f}_{q_{i}/h}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z})$
$\displaystyle=C_{q}(x\tilde{P}^{z},\mu)\exp\biggl{[}\frac{1}{2}\gamma_{\zeta}^{q}(\mu,b_{T})\ln\frac{\tilde{\zeta}}{\zeta}\biggr{]}f_{q_{i}/h}(x,\vec{b}_{T},\mu,\zeta)+\ldots\,,$
(4.1)
$\displaystyle\tilde{f}_{g/h}(x,\vec{b}_{T},\mu,\tilde{\zeta},x\tilde{P}^{z})$
$\displaystyle=C_{g}(x\tilde{P}^{z},\mu)\exp\biggl{[}\frac{1}{2}\gamma_{\zeta}^{g}(\mu,b_{T})\ln\frac{\tilde{\zeta}}{\zeta}\biggr{]}f_{g/h}(x,\vec{b}_{T},\mu,\zeta)+\ldots\,,$
(4.2)

where the ellipses indicate power corrections. The quark matching coefficient
$C_{q}$ and CS kernel $\gamma_{\zeta}^{q}$ are both independent of the quark
flavor.

The factorization formula has many implications. First, when matching quasi-
and continuum TMDs, there is no mixing between quarks and gluons, nor is there
flavor mixing. This means that lattice calculations of TMDs for various
flavors and for gluons should be easier than anticipated. We confirmed the
momentum renormalization-group evolution for the matching coefficient, and
solved it to obtain a explicit result at NLL, confirming it agreed with
earlier fixed-order results in the literature. Finally, our proof has
implications for factorization formulas matching the Lorentz-invariant
approach (MHENS scheme) to lattice TMDs. It implies that ratios of MHENS TMDs
with $b=\vec{b}_{T}$ give direct access to information about the ratio of
Collins TMDs integrated against the matching coefficient
$C_{q}(x\tilde{P}^{z},\mu)$ in the numerator and denominator. The treatment of
factorization for $x$-dependent ratios of MHENS correlators is a bit more
complicated, and we identified the additional ingredients as being from cusp
anomalous dimensions in the lattice renormalization and the need for a
different soft factor subtraction in the relation to the Collins TMD.

Our results are an important step in pushing forward the knowledge we can
extract from the lattice about TMD observables, yet much remains to be done.
We have derived a factorization formula at leading power to all orders in
$\alpha_{s}$. However, the perturbative quasi-to-Collins matching coefficient
$C_{q}$ is only known at one-loop, except for certain logarithmic terms
constrained by its RGE, see section 3.2.2. Recently there has been renewed
interest in TMDs at subleading power, with for example a derivation of the
necessary form of the factorization formula for polarized SIDIS at subleading
power [81, 82]. Finding continuum-to-lattice factorization formulas for
subleading power TMDs would be interesting.

As experiments come online that promise to push our knowledge of hadronic
structure to new depths, the need for corresponding first-principles
predictions becomes ever more clear. The challenges faced in calculating TMDs
on the lattice give us a roadmap for efficiently deriving other key hadronic
properties. Constructing a generic operator encompassing all possible physical
and lattice scheme choices for a specific distribution is useful for
understanding the space of possibilities. Comparing the quasi-to-Collins and
MHENS-to-Collins factorizations, we see that we must walk a fine line between
perturbative and numerical challenges when choosing a lattice observable. From
the lattice standpoint, unless headway is made on the sign problem, matrix
element correlators from which distributions are constructed must employ
Wilson lines on equal-time paths. Wilson lines should have as few sides and
cusps as possible to minimize difficulties with renormalization. From the
start, it is also important to account for lattice renormalization, soft
function subtractions, and finite Wilson line lengths. One must also be
careful with different orders of limits and renormalization, which can often
lead to additional perturbative matching kernels.

In the case of TMDs specifically, it is clear from the phase space of possible
lattice correlators that there are additional freedoms in the definitions that
could still be exploited. Quasi-TMDs and MHENS TMDs provide two examples of
how things can differ due to these choices, and we have advocated for some of
the benefits of the quasi-TMD approach. The calculation of MHENS TMDs has an
excellent track record, and it will be important to continue to carry out
calculations with various lattice TMDs, and confirm consistency amongst the
results. Ultimately, this program has the potential to lead to precise
determinations of the full functional dependence of the eight leading-power
spin-dependent TMDs.

## Acknowledgments

We thank Michael Engelhardt, Xiangdong Ji and Yizhuang Liu for useful
discussions. This work was supported by the U.S. Department of Energy, Office
of Science, Office of Nuclear Physics, from DE-SC0011090, DE-AC02-06CH11357
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. S.T.S. was partially supported by
the U.S. National Science Foundation through a Graduate Research Fellowship
under Grant No. 1745302. Y.Z. is partially supported by an LDRD initiative at
Argonne National Laboratory under Project No. 2020-0020.

##  Appendix A Perturbative cross-checks

###  A.1 Matching unpolarized quark TMDs in different schemes at NLO

In this section, we perturbatively test our results from section 3 by first
constructing the matching between the Collins and JMY TMDs at next-to-leading
order (NLO), and then deriving the same matching between the Collins and LR
TMDs. The latter immediately yields the matching between the quasi- and
Collins TMDs.

####  A.1.1 Matching Collins and JMY TMDs

We first relate the Collins and JMY TMDs. Since the physical cross-section is
scheme-independent, we have

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}Q^{2}\mathrm{d}Y\mathrm{d}^{2}\vec{q}_{T}}$
(A.1)
$\displaystyle=\sigma_{0}\sum_{i,j}H_{ij}(Q,\mu)\int\frac{\mathrm{d}^{2}\vec{b}_{T}}{(2\pi)^{2}}e^{\mathrm{i}\vec{q}_{T}\cdot\vec{b}_{T}}f_{i/h_{1}}\bigl{(}x_{1},\vec{b}_{T},\mu,\zeta_{1}\bigr{)}f_{j/h_{2}}\bigl{(}x_{2},\vec{b}_{T},\mu,\zeta_{2}\bigr{)}$
$\displaystyle=\sigma_{0}\sum_{i,j}H_{ij}(Q,\mu,\rho)\int\frac{\mathrm{d}^{2}\vec{b}_{T}}{(2\pi)^{2}}e^{\mathrm{i}\vec{q}_{T}\cdot\vec{b}_{T}}f^{\mathrm{JMY}}_{i/h_{1}}\bigl{(}x_{1},\vec{b}_{T},\mu,x_{1}\zeta_{v},\rho\bigr{)}f^{\mathrm{JMY}}_{j/h_{2}}\bigl{(}x_{2},\vec{b}_{T},\mu,x_{2}\zeta_{\tilde{v}},\rho\bigr{)}\,,$

where the first and second lines employ the Collins and JMY schemes as given
in eqs. (2.15) and (2.21), respectively. We want to choose a frame where the
TMDs appear symmetrically, such that there is no CS evolution from an
asymmetric split of soft radiation. This can be ensured by choosing
$\zeta_{1}=\zeta_{2}$ and $x_{1}\zeta_{v}=x_{2}\zeta_{\tilde{v}}$, which using
their definitions implies888 The one-loop expressions for eq. (A.1) only agree
if one fixes $x_{1}\zeta_{v}=x_{2}\zeta_{\tilde{v}}$. This may simply be a
result of using this particular frame when calculating the NLO ingredients in
ref. [56]; the factorization theorem in ref. [56] does not make a definite
statement on this. For more details, see appendix B.

$\displaystyle\zeta_{1}=\zeta_{2}=Q^{2}\,,\qquad(x_{1}\zeta_{v})^{2}=(x_{2}\zeta_{\tilde{v}})^{2}=\rho
Q^{2}\,.$ (A.2)

The last equality fixing $\rho$ results from the definition in eq. (2.23). Eq.
(A.1) then implies

$\displaystyle
f^{\mathrm{JMY}}_{i/h}\bigl{(}x_{1},b_{T},\mu,\sqrt{\rho}Q,\rho\bigr{)}$
$\displaystyle=C_{i}^{\mathrm{JMY}}(Q,\mu,\rho)f_{i/h}\bigl{(}x_{1},b_{T},\mu,Q^{2}\bigr{)}\,,$
(A.3)

where $C_{i}^{\mathrm{JMY}}$ is a perturbative kernel which is defined as

$\displaystyle
C_{i}^{\mathrm{JMY}}(Q,\mu,\rho)=\sqrt{\frac{H_{i\bar{i}}(Q,\mu)}{H_{i\bar{i}}(Q,\mu,\rho)}}\,.$
(A.4)

Here, we set $j=\bar{i}$ because this relation must hold for the simplest
cases of flavor-diagonal Drell-Yan scattering. The hard function is
independent of quark flavor, so the matching is also flavor-
independent.999This is violated starting at three loops due to closed quark
loops that couple to the vector current. In order for eq. (A.1) to be true,
such contributions must precisely cancel in the ratio in eq. (A.4). We can
write an asymmetric form of eq. (A.3) using CS evolution,

$\displaystyle
f^{\mathrm{JMY}}_{i/h}\bigl{(}x,b_{T},\mu,\sqrt{\rho}Q,\rho\bigr{)}$
$\displaystyle=C_{i}^{\mathrm{JMY}}(Q,\mu,\rho)\exp\left[\frac{1}{2}\gamma_{\zeta}^{i}(b_{T},\mu)\ln\frac{Q^{2}}{\zeta}\right]f_{i/h}\bigl{(}x,b_{T},\mu,\zeta\bigr{)}\,.$
(A.5)

Here, $C_{i}^{\mathrm{JMY}}$ and $\gamma_{\zeta}^{i}$ differ for quarks and
gluons, but do not depend on quark flavor. Appendix B verifies eq. (A.5) at
one loop; here, we only show one-loop results for the hard function,

$\displaystyle H_{q\bar{q}}(Q,\mu)$
$\displaystyle=1+\frac{\alpha_{s}C_{F}}{2\pi}\left[-\ln^{2}\frac{Q^{2}}{\mu^{2}}+3\ln\frac{Q^{2}}{\mu^{2}}+\frac{7}{6}\pi^{2}-8\right]+\mathcal{O}(\alpha_{s}^{2})\,,$
(A.6) $\displaystyle H_{q\bar{q}}(Q,\mu,\rho)$
$\displaystyle=1+\frac{\alpha_{s}C_{F}}{2\pi}\left[-\ln^{2}\frac{Q^{2}}{\mu^{2}}+3\ln\frac{Q^{2}}{\mu^{2}}+\ln^{2}\frac{\rho
Q^{2}}{\mu^{2}}-2\ln\frac{\rho
Q^{2}}{\mu^{2}}+2\pi^{2}-4\right]+\mathcal{O}(\alpha_{s}^{2})\,,$

see eqs. (B.1) and (B.12). We stress that the second line only holds in the
frame defined by eq. (A.2). Inserting eq. (A.6) into eq. (A.4), we obtain the
one-loop quark matching kernel

$\displaystyle C_{q}^{\mathrm{JMY}}(Q,\mu,\rho)$
$\displaystyle=1+\frac{\alpha_{s}C_{F}}{2\pi}\left(-\frac{1}{2}L_{\rho}^{2}+L_{\rho}-\frac{5}{12}\pi^{2}-2\right)+\mathcal{O}(\alpha_{s}^{2})\,,$
(A.7)

where we abbreviate the logarithm as

$\displaystyle L_{\rho}=\ln\frac{\rho
Q^{2}}{\mu^{2}}=\ln\frac{(2x_{1}P_{1}\cdot v)^{2}}{v^{2}\mu^{2}}\,.$ (A.8)

####  A.1.2 Matching Collins and LR/Quasi TMDs