Transverse-momentum-dependent factorization for lattice observables

Using soft collinear effective field theory, we derive the factorization theorem for the quasi-transverse-momentum-dependent (quasi-TMD) operator. We check the factorization theorem at one-loop level and compute the corresponding coefficient function and anomalous dimensions. The factorized expression is built from the physical TMD distribution, and a nonperturbative lattice related factor. We demonstrate that lattice related functions cancel in appropriately constructed ratios. These ratios could be used to explore various properties of TMD distributions, for instance, the nonperturbative evolution kernel. A discussion of such ratios and the related continuum properties of TMDs is presented.


I. INTRODUCTION
Over the last years, continuous progress in theory and phenomenology of a transverse-momentum-dependent (TMD) factorization theorem made it a valuable tool for analysis and prediction of many observables (for a review see [1]). It has been demonstrated that the TMD factorization approach accurately describes the data in a broad range of energies and a wide spectrum of processes [2][3][4][5][6]. Conceptually, TMD factorization [7,8] is different from the collinear factorization and gives rise to a number of specific novel effects. In this article, we apply TMD factorization to a certain class of operators suitable for evaluation by QCD lattice methods.
There are eight leading twist TMD distributions, each of which depends on a transverse variable and longitudinal momentum fraction x. The purely experimental determination of all of these TMD distributions is a highly nontrivial task. Therefore, the prospects for obtaining complementary information from QCD lattice simulations look extremely promising, in particular, due to the possibility of measuring correlators directly in coordinate space. The latter point is advantageous because the TMD factorization theorem we discuss and TMD distributions are naturally formulated in coordinate space, despite the fact that their interpretation is usually given in momentum space. From experimental data, one can extract coordinate space information only via a Fourier transformation, resulting in a significant systematic error and model-bias. A good example for the encountered problems is the TMD evolution kernel DðbÞ (also known as Collins-Soper (CS) kernel [9]). To extract the CS kernel from data, one has to combine data from many experiments performed at varying energies. The current global pool of data gives access to energies from 1 to 150 GeV [5]. However, the precision of most of these data is quite limited, and their interpretation depends nontrivially on D. The later is known up to α 3 s -order in perturbation theory [10], but is poorly constrained beyond perturbative values of b. Even the shape of DðbÞ is questionable (compare for instance the extractions made in [5] and [6], and their discussion in Ref. [11]). This problem can be resolved, or at least reduced, by lattice simulations.
Suggestions for lattice studies of TMD observables were made long ago [12,13]. At that time, however, some crucial assumptions were rather conjectural. Recently, such efforts were promoted to a higher level with the formulation of appropriate factorization theorems [14][15][16]. In all cases, one considers an equal-time analog of a TMD operator, which turns into an ordinary TMD operator after the boost. In this paper, we present a different analysis of the same operator within the TMD factorization approach, based on the q T -dependent soft-collinear effective field theory (SCET II). We demonstrate that the suggested lattice observables are more closely related to TMD hadron tensor rather than to the TMD distributions. Using the TMD hadron tensor, we present a construction that stresses the analogy between lattice observables and physical quantities used in the description of conventional processes like the Drell-Yan process, utilizing the same terminology. Although the route of the derivation of the factorized expression differs from that of [15,16], we arrive at an equivalent result. Checking the factorized expression at the one-loop level, we found that the perturbative parts are in complete agreement with [15].
The expression derived here is only the leading term of the factorization theorem. The subsequent terms are formally suppressed by powers of hadron momentum. This fact should not be over interpreted because a closer analysis reveals the potential breakdown of this expansion. Namely, each next order term has a stronger small-x singularity than the preceding ones. Such a problem is quite common for factorization theorems for lattice observables. For example, the quasiparton distribution functions (PDF) [17], and pseudo-PDFs [18] also suffer from this problem, as is discussed, e.g., in Ref. [19]. In the TMD case, the small-x divergent terms are more troublesome because they are enhanced by TMD evolution. This is unavoidable since the hard scale is the parton momentum ∼xP leading to strong factorization breaking at small x.
The paper is split into three sections. In Sec. II we define matrix elements suitable for lattice simulations and derive the factorized expressions using SCET II. In Sec. III we check factorization at the one-loop level and derive the corresponding coefficient function at the next-to-leading order (NLO). Finally, in Sec. IV we discuss ratios of matrix elements that have simpler properties and thus could serve to determine TMD distributions with less effort. We emphasize, in particular, the advantages of ratios at small-longitudinal separation.

A. Definition of lattice observables
The considered lattice observables read W ½Γ f←h ðb; l; L; v; P; SÞ ¼ where jP; Si is a single-hadron state with momentum P and spin S, and Γ is a Dirac matrix. The hadron species and the flavor of the quark field are specified by the labels h and f. ½x; y is the straight gauge link between points x and y, The same object has been considered in [14][15][16]20,21]. Often, the matrix element (1) is called a quasi-TMD distribution in analogy to the quasiparton distribution functions [17]. However, as we demonstrate in the next section, the structure of (1) does not remind of a TMD distribution but rather of a TMD hadronic tensor, such as the hadronic tensor for Drell-Yan or Semi-Inclusive Deep Inelastic Scattering (SIDIS). For that reason, we avoid the term quasi-TMD distribution, and denote (1) by the letter W. The space-time orientation of different quantities in W is given in Fig. 1. The lattice operator must be an equaltime operator and thus the vectors v μ and b μ do not have time components. Consequently, b 2 < 0 and v 2 < 0.
The vectors P μ and v μ define an analog of the scattering plane. The transverse vector b μ is orthogonal to them, The vector b μ defined by (3) is restricted to a line, due to the absence of time components (except for the special case that the vectors v μ , P μ and the time direction lie in a plane). The situation is different for physical kinematics where the scattering plane is formed by two timelike vectors and thus the vector b μ is restricted to a plane.

B. Factorization limit
The clear separation of collinear and soft-field modes within the hadron is a prerequisite for any TMD factorization theorem. It can be achieved by considering a fast moving hadron, for which anticollinear components of field momenta are suppressed in comparison to collinear ones. To quantify this condition we write the momentum of a hadron as where M is the mass of the hadron, n andn are lightlike vectors n 2 ¼n 2 ¼ 0 (see also Fig. 1), normalized according to ðnnÞ ¼ 1. Here, we use the standard notation of light-cone components of a vector a μ : So, the factorization limit requires with λ being the generic small parameter of SCET. In this regime the hadron momentum is almost lightlike. We also assume that the staple-shaped gauge links contour (1) is much longer than broad Under this assumption, the effects caused by the interaction with the transverse gauge link ½b þ Lv; Lv are suppressed as b=L and l=L, and thus can be neglected. We then introduce an "instant" (formally infinitely heavy) scalar field, HðxÞ, with the Lagrangian and approximate the gauge links ½x; x þ Lv by the H propagator. In Eq. (8), D μ is the QCD covariant derivative. The field H differs from a usual scalar heavy quark field [22] only by the fact that v 2 < 0. (This is why we call it "instant" field.) In the notation (8), the similarity of the matrix element (1) with the ordinary hadron tensor for TMD factorization becomes transparent. We rewrite (1) as where W ½Γ ¼ 1 2 TrðWΓÞ and J i is the heavy-to-light current The structure of the first term in (9) resembles the structure of the hadron tensor for TMD cross sections, with the main difference that there is only a single hadron. The second hadron is replaced by the "instant" field H.

C. Field modes factorization and SCET current
The analogy with the TMD hadron tensor (9) allows us to recapitulate the main points of TMD factorization and apply it to the lattice case. To derive the factorized expression, we use the soft-collinear effective field theory with finite q T (SCET II) approach, similarly to Ref. [8].
In SCET II, one distinguishes collinear, anticollinear and soft fields. In leading approximation, the fast-moving hadron is a composition of collinear fields (ξ for quarks and A c;μ for gluons). Their momentum components are of the structure The separation of transverse and collinear modes requires a clear hierarchy between corresponding momentum components. In the present case, the typical transverse momenta in diagrams are ∼b −1 . And, therefore, we have an additional constraint A priori it is not evident how to count the field H in terms of λ since a Wilson line does not carry a momentum. However, the situation becomes clear if one boosts the system such that P − → 0. Then the Wilson line H is turned towards the light-cone direction n. In the boosted frame, the v þ -component of v μ can be ignored and the field H can be approximately considered as an anticollinear Wilson line. Therefore, the fields H and collinear fields have no direct interaction but only through soft exchanges.
Using these counting rules we write the leading power SCET operator that corresponds to the current J i (10) as In this expression the collinear Wilson line contains all gluons radiated by the collinear quark field. The anticollinear Wilson line Wn contains the gluons radiated by the H-field and is given by a similar expression with n →n. The Wilson lines Y are the result of the decoupling transformation [23]. They have analogous expression to (14) but build with soft gluon fields. The coefficient C H is the matching coefficient between SCET and QCD operators. It depends on the momentum of field ξ (in position spacep is an operator), and is independent of quark flavor. At leading order C H ¼ 1.
Expression (13) applies for ðvPÞ > 0. If ðvPÞ < 0 the gluon fields are summed with the opposite sign [in comparison to (13)] and thus they form Wilson lines pointing to þ∞n. So, for the case ðvPÞ < 0 the SCET operator reads Combining (13) and (15) we express the QCD current (10) by

D. Factorized expression
Substituting the effective currents (17) into the expression (1) we obtain [here for ðvPÞ > 0] W ½Γ f←h ðb; l; L; v; P; SÞ Here, we have performed a Fiertz transformation to recouple color indices, and have dropped color-covariant structures. The collinear, anticollinear and soft fields operate on different Hilbert spaces, such that the total Hilbert space can be written as a direct product of three distinct Hilbert spaces [24,25]. Doing so one has to correct for overlap of the field definitions in the soft region. The overlap contribution can be removed by the so-called zero-bin subtraction factor [26]. Additionally the fields can be Taylor-expanded in the slow (in comparison to other components) directions, which are determined by the counting rules (11). After these operations we obtain the following result: where The functions areΦ where we use the QCD fields since within its own Hilbert space each sector of SCET is equivalent to QCD. The zerobin factor (denoted as Z.b.) removes the contribution from the overlap of soft and collinear modes. It is not known explicitly except for certain regularizations, for instance for δ-regularization [8] where it is equivalent to the TMD soft factor. The functionΦ is (a Fourier transform of) an unsubtracted TMD distribution. The functionΨ has an analogous structure. The only difference is that it measures a TMD distribution of the field H. For that reason we call it an (unsubtracted) instant-jet TMD distribution. The function S is the TMD soft factor. The expression (19) applies for ðvPÞ > 0 and thus,Φ and S correspond to Drell-Yan kinematics. If ðvPÞ < 0 the Wilson line along n points to þ∞n, andΦ and S correspond to SIDIS kinematics (and the coefficient function is replaced byC H ).

E. Recombination of rapidity divergences
Unsubtracted TMD distributions have rapidity divergences that appear due to the presence of infinite lightlike Wilson lines separated in the transverse plane. Rapidity divergences are associated with the directions of Wilson lines. In the current case, there are two light-cone directions, and thus we introduce two regularization parameters δ þ and δ − . These parameters regularize rapidity divergences associated with the directions n andn, correspondingly. In the product of all functions in (19) rapidity divergences cancel, and therefore, the last step of the factorization approach is to recombine rapidity divergences, and to introduce physical (aka finite) TMD distributions. In this procedure we follow Ref. [27]. In what follows, we use δ-regularization for rapidity divergences, but the same procedure can be performed for different kinds of regulators. The final result is independent of the regularization used.
In the expression (19) the rapidity divergences are present according to the following pattern (in this section we omit all arguments of the different functions, except the ones related to rapidity divergences): In Ref. [27] it has been shown that rapidity divergences are structurally equivalent to ultraviolet divergences, and therefore, can be absorbed into a divergent factor R. Introducing rapidity renormalization factors into (24) we obtain where ν AE are the scales of rapidity-divergence renormalization, and The function S 0 depends on ν 2 ¼ 2ν þ ν − due to Lorenz invariance. This expression is independent of ν AE by definition, and each function here is finite. The dependence on ν AE is given by the renormalization group equation where D is the rapidity anomalous dimension [27], or CS-kernel [9]: The equation forΨ is analogous.
Introducing the boost-invariant variables where P þ is the collinear component of hadron momentum and μ is a factorization scale. The parameter ζ is the standard rapidity evolution parameter [7,8,27]. The parameterζ is the analogous parameter for Ψ. Let us emphasize that the expression forζ is unusual, because typically the scale of rapidity divergences is associated with the collinear component of a momentum. However, there is no momentum that is associated with the field H. The only momentum scale presented in Ψ is the factorization scale. In Sec. III, we confirm (28) by a one-loop calculation. The dependence of the functionΦðζ; ν 2 Þ on ζ follows from (26) and reads The function Ψ depends onζ in the same way. Generally, the function S 0 ðν 2 ; bÞ is a process-dependent and (at large-b) nonperturbative function. To get rid of it, we note that the variable ν 2 decouples from the evolution, and thus the function S 0 can simply be absorbed into the definition of a TMD distribution. In fact, the physical definition of TMD distributions already includes such factors (see discussion in [7,27]). They are built from the remnants of TMD soft factors. So, a physical TMD distribution, such as the one used to describe the Drell-Yan or SIDIS processes, is defined together with an appropriate It is independent on ν 2 [27]. To formulate the factorization in terms of physical TMD distributions, we use the definition (30) and compensate the extra factor ffiffiffiffiffiffiffiffiffiffi ffi S TMD 0 p by an appropriate redefinition of the instant-jet TMD distribution: Note that in a suitably defined regularization scheme (for instance δ-regularization [8]), the zero-bin subtraction factor Z:b ¼ S 2 ðbÞ, and thus S TMD 0 ¼ S 0 . However, generally, these factors can be different in the nonperturbative regime.

F. Final form of the factorized expression
The final form of the factorized expression reads W ½Γ f←h ðb; l; L; v; P; S; μÞ Here, we restored all arguments of functions including the scale μ.
There are two points in equation (32) which need clarification. The first point concerns the dependence of Ψ on the variable l. The variable l appears in Ψ accompanied by the lightlike vectorn, and thus can enter the function only in a scalar product with some other vector of the problem. This can only be the vector v, and thus the dependence on l can appear only as ðv þ v − lΛ QCD Þ or as ðv þ v − l=LÞ in the presence of a regularization parameter L [compare to the function Φ where the vector l enters via ðlv − P þ Þ]. However, in the factorization limit both of these combinations are negligible. Thus, we conclude that the dependence on l is marginal, This statement is also clear from the boosted frame perspective: boosting P − =P þ → ∞ one automatically gets v þ =v − → 0. So, the dependence on l is negligible, unless l is very large. Note that the function Ψ defined in (31) where r depends on Γ and the Lorentz structure of the TMD. For example, Γ ¼ γ þ , r ¼ 0ð1Þ for unpolarized (spin-flip) TMDs.
Using these facts we rewrite (32) as W ½Γ f←h ðb; l; L; v; P; S; μÞ where P v ¼ v − P þ . This is probably the most practical form of the factorization theorem, and we will use it later. The factorization statement is independent on the Dirac structure, which is standard for the TMD factorization approach. Therefore, using this expression one can describe polarized and unpolarized processes equally well.
It is important to emphasize that the size of power corrections in (36) significantly depends on x. In fact, the typical momentum scale entering factorized expressions iŝ p ∼ xP rather then just P. Therefore, a more reliable estimation of the power corrections is OðP − =x 2 P þ Þ. This is a typical size of power corrections to factorization theorems for lattice observables; see e.g., the case of quasi-PDF power correction which are of order OðΛ 2 =x 2 ðpvÞ 2 Þ as is shown in Ref. [19]. Such a large power correction can undermine the applicability of the whole approach as we discuss below.

III. NLO EXPRESSIONS
In this section we present the computation of elements of the factorization theorem at one loop. The calculation supports the correctness of the construction. The calculation presented here is done in the δ-regularization scheme [8,28], which allows us to reuse results of earlier calculations made in [28,29]. Our results coincide with those of [14], where they were reached in a different manner.

A. Hard matching coefficient
To evaluate the matching coefficient between the QCD current (10) and the SCET current (13) one needs to compute and compare both sides of equations (17). At the same time, one should demonstrate cancellation of collinear and soft divergences. We use δ-regularization for collinear and soft divergences and dimensional regulariza- In the δ-regularization scheme [8,28] the zero-bin subtraction factor coincides with the soft factor squared Z:bj δ−reg ¼ S 2 ðbÞ: Therefore, the matching relation at NLO turns into where we omit arguments for simplicity, and use the following shorthand notation for coefficients of perturbative We have also used that There are four diagrams that contribute to (38), presented in Fig. 3. There exist other diagrams contributing to each term of (38), (these are various selfenergy diagrams), but their contributions exactly cancel in the sum.
In δ-regularization the diagram J QCD reads where Δ and Δ v are parameters of δ-regularization (Δ > Δ v > 0). Evaluating this diagram in the limit Δ, Δ v → 0 and p þ ≫ p − we obtain where the dots stand for power suppressed contributions ∼Δ. The i0-terms are important for proper analytic continuation between the cases ðpvÞ > 0 and ðpvÞ < 0. It reads ð−ðpvÞ − i0Þ ¼ jpvje iðargðpvÞ−πÞ . The first term in brackets represents the soft divergence, whereas the second and the third term are collinear and anticollinear divergences. Note that ϵ > 0 throughout and thus factors Δ −ϵ are divergent at Δ → 0.
Evaluating analogously the rest of the diagrams (note that the results for Φ ½1 and S ½1 in δ-regularization can be found in Refs. [8,28] and [8,29], correspondingly) we obtain where the upper sign corresponds to the geometrical configuration with ðvPÞ > 0 and the lower sign corresponds to configuration with ðvPÞ < 0. The regularized propagators reproduce soft propagators in the soft regime; therefore, the parameters Δ and Δ v are related to δ AE according to The sign of v − is the same as the sign of ðpvÞ, and thus δ þ > 0. Substituting these expressions into (38) we observe that each divergent sector cancels exactly (i.e., at all orders of ϵ-expansion). Clearly, it is important to keep the proper direction of Wilson lines in mind, which leads to different signs of ðvPÞ resulting in these cancellations. Altogether, this confirms the derived factorization theorem at NLO. The hard-matching coefficient is Performing renormalization in the MS-scheme we obtain C H vp μ where L ¼ lnðð2jpvjÞ 2 =μ 2 Þ. Importantly, the coefficient function is the same for ðpvÞ > 0 and ðpvÞ < 0 at this perturbative order. Nonetheless, the continuation between these regions is nontrivial, and for higher orders the coefficient functions could be different. The expression (45) coincides with the one derived for the hard coefficient function in Ref. [14,15], where the calculation has been done differently. We have also performed the same computation with a finite-length H-Wilson line, as in (1) also takes place, although the matching relation between δ þ and L depends on ϵ and does not hold at higher orders of perturbation theory.

B. Anomalous dimensions
The where γ F and γ Ψ are ultraviolet anomalous dimensions, and D is the rapidity anomalous dimension. The integrability condition for these equations gives the Collins-Soper equation where Γ cusp is the cusp anomalous dimension for lightlike Wilson lines. The solution for ultraviolet anomalous dimensions is The anomalous dimension γ V is known up to α 3 s -order, and at LO γ V ¼ −6C F α s =ð4πÞ. The anomalous dimension γ Ψ is the finite part of the heavy-to-light anomalous dimension. Since the vertex diagram (42) contributes only to the double logarithm structure, the finite part of the anomalous dimension is twice the anomalous dimension of a heavy quark field [30] All components of the current J i are renormalized by a single renormalization factor, J ren: i ¼ Z J J i , and thus the matrix element W is renormalized by The corresponding anomalous dimension, is evaluated in [31] at NNLO (for v 2 > 0) and reads Note that the anomalous dimensions γ J and γ Ψ are known from the literature and are related to heavy quarks physics (i.e., with v 2 > 0). They agree with those presented here at LO. However, they could disagree at higher perturbative orders due to v 2 < 0 kinematics. The renormalization group requires Substituting LO anomalous dimensions and the NLO coefficient function (45), we check that this relation is satisfied if the ζ-parameters are related by This fixes the relative freedom in the definition of boostinvariant variables ζ andζ. It also confirms our observation that in the absence of momentum the rapidity divergences inΨ are weighted by the factorization scale μ (28).

IV. RATIOS OF LATTICE OBSERVABLES
The factorized expression (36) has the generic form of a TMD factorization theorem, and, therefore, it incorporates three nonperturbative functions. These are the TMD distribution Φ, the instant-jet TMD distribution Ψ, and the rapidity anomalous dimension D. The latter is not explicitly presented in the formula, but enters via the scaling properties of the distributions, see (46)-(48). To determine these functions one needs to measure W in a large range of parameters P and l. However, even in this case the function Ψ, which depends on l only weakly (33), would be totally correlated with D.
There are two principal ways to bypass this problem. The first approach is to obtain the values of Ψ in an independent calculation. This could be done in perturbation theory (at small values of b) [15], or performing a separate lattice calculation, such as the one suggested recently in Ref. [32]. The second approach is to consider ratios of lattice observables, such that undesired factors (in particular the function Ψ) cancel. This approach looks more promising because the measurements of ratios are simpler on the lattice. In addition to the cancellation of Ψ one would also profit from the cancellation of various other multiplicative factors such as lattice renormalization constants, and a corresponding reduction of systematic uncertainties for the lattice results.
In this section, we consider ratios of the form In such ratios the contribution of Ψ cancels, as well as a common μ-dependence. Various properties of these ratios have been considered in [13,14,21,33,34]. In the subsequent sections, we discuss particularly interesting combinations of parameters in R, that were not yet mentioned in the literature. These combinations allow to check the validity of the factorization theorem, and estimate the joined systematic uncertainties of the approach and the lattice computation. Additionally, we consider the case l ¼ 0. The l ¼ 0 case is particularly simple to simulate on the lattice, as has been done in Refs. [33,34] for TMD double moments. We show that it grants access to the nonperturbative anomalous rapidity dimension. For brevity of the formulas in this section we denote only those arguments of W that are different, assuming that all the other arguments in the ratio (58) are the same. Also, we universally denote all power corrections by A. Sign-flip The most elementary test of the factorization theorem (36) is the measurement of the famous sign-flip of P-odd TMD distributions between the Drell-Yan and SIDIS process [35]. SIDIS and Drell-Yan kinematics are distinguished by the sign of ðvPÞ. Therefore, the sign-flip can be tested by replacing v → −v.
For example, considering Γ ¼ γ þ . We have two Lorentz structures where ϵ μν T ¼ ϵ þ−μν . These structures can be independently extracted from a lattice simulation [33], and are proportional to unpolarized f 1 and Sivers f ⊥ 1T TMD distributions. The Sivers distribution is P-odd and its sign depends on the direction of the Wilson lines, in contrast to unpolarized distributions. Therefore, the following ratios should hold: These relations are trivial at NLO due to the independence of coefficient jC H j 2 of the sign of ðvPÞ. However, they could be violated by higher perturbative terms, if there are nontrivial effects of analytical continuation in ðvPÞ. Similar measurements have been performed in [33,34], where the ratios W ⊥ 1T =W 1 (and similar for Γ ¼ iσ þμ γ 5 ) has been studied at different values of b and P and for different signs of v. Perfect agreement with (61) has been demonstrated.

B. Power suppressed terms
A great feature of lattice QCD is the possibility to measure objects unaccessible in an experiment directly. In particular, one can compare measurements of different Lorentz structures and check the TMD factorization theorem in a completely controlled environment. The Dirac structures of higher TMD twist must be suppressed due to dominance of the collinear components in the hadron. We have where Γ 0 (20). Despite the apparent triviality of this statement the numerical evaluation of this ratio on the lattice is very important because it allows to estimate systematic uncertainties. In a sense, it directly measures the size of power corrections to the factorization theorem (36). This is very valuable information, because the accessible hadron momenta in state-of-the-art lattice simulations are at most a few GeV, such that power corrections can be large.

C. Nonperturbative rapidity anomalous dimension
The most exciting property of the ratios R is their exclusive sensitivity to the rapidity anomalous dimension, which was also pointed out in Refs. [14,16,21]. The properly constructed ratio is almost independent of nonperturbative functions except the rapidity anomalous dimension. To extract D one needs the ratio of W 's at different momenta Using (36) where Note that the scale μ is taken to be the same in the numerator and the denominator in order to cancel the unknown functions Ψ. Evolving the remaining two functions in ζ to the same point ζ 0 , and partially canceling evolution factors we get A similar ratio is considered in detail in Ref. [21], where it is suggested to calculate the Fourier images of denominator and numerator separately. Such a method has bright prospects, but is noticeably more demanding on the lattice side than the one discussed below. The main difficulty comes from the noncancellation of lattice renormalization factors that, therefore, have to be computed separately. An additional, but not smaller, problem comes from the Fourier transformation in l, which has to be deduced from the few measurable points with l ≪ L; Λ −1 QCD . Altogether, these problems could result in a large systematic uncertainty.
To avoid these difficulties we suggest to consider the case of l ¼ 0. Roughly speaking, the plain l ¼ 0 case corresponds to the ratio of the first Mellin moments of TMD distributions. The higher moments can be accessed by taking derivatives with respect to l. Let us denote where the prefactor is chosen such that at b → 0 the ratios become unity R ðnÞ → 1. These ratios give direct access to the rapidity anomalous dimension, as we demonstrate below. The lattice computation of the n ¼ 1 case is relatively straightforward, and some exploratory computations were already made in [13,33,34]. The consideration of higher n is more complicated due to the mixture of operators with different Dirac structures (especially if b is not much smaller than L; see also [36]), and possible problems with restoration of rotational symmetries. Nonetheless, we expect that these problems could be overcome and, at least, the case n ¼ 2 is feasible. The ratio (66) at l ¼ 0 is convenient to present in the form where r ðnÞ ¼ 1 þ Oðα s Þ. The NLO contribution to r ðnÞ is obtained using (45) It is straightforward to check that this expression is independent on μ and ζ 0 . Therefore, it can be further simplified by using the optimal definition of a TMD distribution [37]. Setting ζ 0 ¼ ζ μ ðbÞ such that Φ where Φ ½Γ ðx; bÞ (without scale) is the optimal TMD distribution. The optimal definition is convenient and often used for phenomenological extractions [3][4][5]37,38]. Let us mention, that the NNLO expression for r ðnÞ can be derived using only the NLO anomalous dimensions and the finite part of jC H j 2 at NLO. Thus, it could be possibly reconstructed from already existing calculations. In Fig. 4 we plot M ðnÞ;unpol ln jxj as function of b for different values of n starting from n ¼ 2, using values from [5]. In Fig. 5 we show the function R ð2Þ for typical lattice momenta P þ 1 ¼ 2.5 GeV and P þ 2 ¼ 2 GeV (and jv − j ¼ 1= ffiffi ffi 2 p ). The scale μ is set to bē such that the logarithm in (68) is zero. As input we take two of the most recent extractions of unpolarized TMDPDFs and rapidity anomalous dimensions [5,6]. The computation is done with the ARTEMIDE package [3]. The uncertainty band is due to the uncertainty of the phenomenological parameters. The considered models have essentially different behavior at large values of b, which should be clearly distinguishable on the lattice. Let us also note that at small values of b, the limit R ðnÞ → 1 is definitely violated. This is due to the ðjbjP þ Þ −1 correction that is part of OðλÞ.
To extract the rapidity anomalous dimension D from the lattice data, we have to deal with the function M ln jxj . The value of M is difficult to estimate by lattice methods. However, we argue that the effects caused by nontrivial M can be neglected with reasonable accuracy. First of all, we mention that the corrections to 1 in r ðnÞ have an extra suppressing factor lnðP þ 1 =P þ 2 Þ. Together with α s =ð4πÞ it guarantees that these corrections are numerically small. For example, for the parameter values used in Fig. 5 the contribution of the term with M ð2Þ is ∼0.06. Second, the function M has minor dependence on b. In Fig. 5, the contribution of maximum and minimum values of M ð2Þ (see Fig. 4) differ by ∼0.008, which is a tiny number in comparison to the expected accuracy of lattice computations. Based on this observation, we conclude that the function M ðnÞ could be replaced by a constant, adding a ∼1% systematic uncertainty. This constant can be estimated using the existing phenomenological extractions (with ∼1-2% of systematic uncertainty), or it can simply be neglected (implying a ∼10% systematic uncertainty). Therefore, the proposed method allows the determination of R ðnÞ within a few percent of systematic uncertainty, depending on the selected strategy. Alternatively, one can estimate (the constant) M ðnÞ from the lattice data at a single point.
Such simplified schemes should be applied with caution, because M has different behavior for different types of TMD distributions and different n. In the majority of cases, M ln jxj can be approximated by a constant (in b), however, in some cases not. It is clear that M ln jxj is closer to a constant if the integrand has better convergence properties at x → 0. However, for some cases the integrals in M ln jxj are divergent, such that these cases cannot be used for analysis. We should also keep in mind that higher perturbative terms contain M ln n jxj , and have worse convergence properties.
There are two sources of small-x divergence in M: The first one is the factor jxj −2D . The rapidity anomalous dimension D is greater then zero for b ≳ 2e −γ E =μ. Its asymptotic behavior is unknown, although typically it is expected to be a monotonously growing function. Additionally, the value of D also increases with the increase of μ. The uncertainty of the large-b behavior of modern extractions of D are quite drastic [5,6,11,38]. Nonetheless, all recent extractions agree that D > 1=2 for b ≳ 3-4 GeV −1 (here μ ∼ 2 GeV). Therefore, in this range the factor jxj −2D is singular.
The second one is the TMD distribution itself. Generally, at small-x TMD PDFs behave as x α . The value of α depends on the kind of TMD PDF, as has been studied in [39,40]. It has been shown (in the large-N c approximation) that α < −1 for the unpolarized structure Γ ¼ γ þ , −1 < α < 0 for the helicity structure Γ ¼ γ þ γ 5 , and α > 0 for the transversity structure Γ ¼ σ þμ . In each case only the leading distribution has been considered (i.e., f 1 ,g 1L FIG. 5. Comparison of the functions R ð2Þ evaluated with two different phenomenological models for rapidity anomalous dimension: "SV19" and "Pavia19" that are considered in [5] and [6], correspondingly. The uncertainty band is obtained by variation of model parameters within their uncertainty range. and h 1 ). One can expect weaker singularities with a similar general hierarchy for other distributions (i.e., f ⊥ 1T , g 1T , etc). The power of the small-x singularity also depends on the flavor combination. In particular, the nonsinglet combinations have weaker small-x behavior (see e.g. [39]).
In this way, there is a certain hierarchy of R for different Γ and n, such that some give simpler access to D that other, see Table I. Note that the convergence properties improve for nonsinglet flavor combination, and for larger n. In particular, the unpolarized case is convergent for n ¼ 2 in a large range of b, which is also seen in Fig. 4.
Finally, let us stress again that the rapidity anomalous dimension is universal. Therefore, the ratios R ðnÞ should be almost independent of quark flavor, Dirac structure Γ, hadron type, and the momentum parameter n (for convergent cases). The difference between all these cases is only due to functions M ðnÞ .

V. CONCLUSION
In the present article, we have considered quasi-TMD operators, that can be investigated on the lattice. We pointed out the similarity of the lattice observable to the hadronic tensor of TMD processes, such as Drell-Yan or SIDIS. Using the method of soft-collinear effective field theory (SCET II) we derived the factorized expression for the lattice hadronic tensor in terms of physical TMD distributions, and the new instant-jet TMD distribution Ψ defined in (21), (31). The factorized expression generally coincides with expressions derived in [15,16], although the route of derivation is different. We have checked factorization at one-loop level and derived the hard matching coefficient at this order, which coincides with the one derived in [15]. The LO anomalous dimension can be extracted from the literature related to heavy-quark physics, and that value coincides with the results of our calculation. The present derivation is done for arbitrary Dirac structure, and can be easily extended to other interesting cases, such as gluon operators.
Since the factorization formula contains an unknown nonperturbative function Ψ, it is advantageous to consider the ratios of lattice observables with the same geometrical parameters of the operators (i.e., l, b, L and v). In this case, many troublesome factors, such as Ψ and lattice renormalization factors, cancel. The remaining parameters, namely the Dirac structure Γ, hadron momentum, spin and flavor, are enough to extract valuable information on TMD distributions and to estimate the uncertainties of the method. In particular, we pointed out that the ratio of the first derivatives at l ¼ 0 with different hadron momenta can be used to accurately determine the rapidity anomalous dimension (Collins-Soper kernel). In this case, one does not need to evaluate Fourier transformations with respect to l, as suggested in [14]. Evaluating the ratios (62) for suppressed and unsuppressed Dirac structures allows to estimate the systematic uncertainty of the method by lattice simulations.
The hard scale of the derived factorization theorem is the hadron momentum P. Thus, one could expect that the corrections to the factorized term are P −1 -suppressed. However, this is only a crude estimate because the parton fields carry only a fraction of the total hadron momentum. Therefore, the true factorization scale is the parton momentum xP, which is generally much smaller. In contrast to the scattering processes, where the parton momentum is detected, lattice simulations involve all possible parton momenta. This leads to problems caused by low-x divergences. In particular, the power corrections to lattice factorizations are 1=x 2 -enhanced [19]. This observation limits the application range of such factorization approaches. In particular, in the l ¼ 0 case (that was considered in [33,34]), the size of corrections is very strongly dependent on the operator. In certain cases (for instance unpolarized operators) already at NLO level one can encounter divergences. However, one is free to avoid such cases when determining D.
Quite generally, many complications can be avoided if one considers the ratios of lattice observables. The information that could be extracted from ratios is limited, but still valuable. For instance, one can extract the nonperturbative rapidity anomalous dimension. The most simple observable in this case is the ratio of the first (and possibly the second) Mellin moments of quasi-TMDs at different hadron momenta. This ratio is almost exclusively dependent on D. In addition, it depends on ln jxj ðb; μÞ evaluated on TMD PDFs of different kinds. The estimation of ranges for b is made using values from [5,39,40], and gives only the general scale.

Γ
Kind Name Comment
Convergent in some large range of b.