One-loop renormalization of staple-shaped operators in continuum and lattice regularizations

In this paper we present one-loop results for the renormalization of nonlocal quark bilinear operators, containing a staple-shaped Wilson line, in both continuum and lattice regularizations. The continuum calculations were performed in dimensional regularization, and the lattice calculations for the Wilson/clover fermion action and for a variety of Symanzik-improved gauge actions. We extract the strength of the one-loop linear and logarithmic divergences (including cusp divergences), which appear in such nonlocal operators; we identify the mixing pairs which occur among some of these operators on the lattice, and we calculate the corresponding mixing coefficients. We also provide the appropriate RI'-like scheme, which disentangles this mixing nonperturbatively from lattice simulation data, as well as the one-loop expressions of the conversion factors, which turn the lattice data to the MS-bar scheme. Our results can be immediately used for improving recent nonperturbative investigations of transverse momentum-dependent distribution functions (TMDs) on the lattice. Finally, extending our perturbative study to general Wilson-line lattice operators with n cusps, we present results for their renormalization factors, including identification of mixing and determination of the corresponding mixing coefficients, based on our results for the staple operators.


I. INTRODUCTION
One of the main research directions of nuclear and particle physics is the study of the rich internal structure of hadrons, which are the building blocks of the visible Universe. Quantum chromodynamics (QCD) is the theory governing the strong interactions, which are responsible for binding partons (quark and gluons) together into hadrons. Despite the various theoretical models that have been developed for the investigation of hadron structure (e.g., diquark spectator and chiral quark models), an ab initio calculation is desirable to capture the full QCD dynamics. Due to the complexity of the QCD Lagrangian, an analytic solution is not possible, and numerical simulations (lattice QCD) may be used as a first principle formulation to study the properties of fundamental particles.
Distribution functions consist of a set of nonperturbative quantities that describe hadron structure and have the advantage of being process-independent and accessible both experimentally and theoretically. They are expressed in terms of variables defined in the longitu-Electronic addresses: a marthac@temple.edu, b haris@ucy.ac.cy, c spanoudes.gregoris@ucy.ac.cy dinal and transverse directions with respect to the hadron momentum. Based on this, the distribution functions may be classified into parton distribution functions (PDFs), generalized parton distributions (GPDs) and transverse-momentum dependent parton distribution functions (TMDs). Important information is still missing for all three types of distributions: The most well-studied are PDFs, which are single-variable functions, while the TMDs are only very limitedly studied due to the difficulty in extracting them experimentally and theoretically. However, TMDs are crucial for the complete understanding of hadron structure as they complement, together with GPDs, the three-dimensional picture of a hadron.
Due to their light cone nature, distribution functions cannot be computed directly on a Euclidean lattice and typically are parametrized in terms of local operators that give their moments. The distribution functions can thus be recovered from an operator product expansion (OPE), which is, however, a very difficult task: Signal-to-noise ratio decreases with the addition of covariant derivatives in the operators, and an unavoidable power-law mixing under renormalization appears for higher moments. Nevertheless, information on distribution functions (mainly PDFs and, to a lesser extent GPDs) from lattice QCD was obtained from their first moments, via calculations of matrix elements of local operators. These moments are directly related to measurable quantities, for example, the axial charge and quark momentum fraction.
Novel approaches for an ab initio evaluation of distribution functions on the lattice have been employed in recent years. In these approaches, nonlocal operators, including a Wilson line, are involved. While local operators have been used extensively in perturbative and nonperturbative calculations, nonlocal operators were limitedly studied. In particular, calculations using nonlocal operators with Wilson lines in a variety of shapes appear in the literature within continuum perturbation theory. Starting from the seminal work of Mandelstam [1], Polyakov [2], Makeenko-Migdal [3], there have been investigations of the renormalization of Wilson loops for both smooth [4] and nonsmooth [5] contours. Due to the presence of the Wilson line, power-law divergences arise for cutoff regularized theories, such as lattice QCD. It has been proven that in the case of dimensional regularization and in the absence of cusps and self-intersections, all divergences in Wilson loops can be reabsorbed into a renormalization of the coupling constant [4]. Wilson-line operators have been studied with a number of approaches, including an auxiliary-field formulation [6,7], and the Mandelstam formulation [8]. Particular studies of Wilson-line operators with cusps in one and two loops, can be found in Refs. [5,9,10]. There is also related work, in the context of the heavy quark effective theory (HQET) 1 [12][13][14][15], including investigations in three loops [16].
Computations of matrix elements using nonlocal operators with a straight Wilson line have been revived in lattice QCD and phenomenology mainly due to their connection to PDFs via the quasi-PDFs approach proposed by X. Ji [17] 2 . Several aspects of the properties of nonlocal Wilson-line operators have been addressed, including the feasibility of a calculation from lattice QCD [22][23][24][25], their renormalizability [26][27][28][29][30][31][32] and appropriate renormalization prescriptions [33][34][35]. The renormalization has proven to be a challenging and delicate process in which a number of new features emerge, as compared to the case of local operators: There appears an additional power-law divergence, and the matrix elements are nonlocal and contain an imaginary part.
While information on physical quantities is obtained from hadron matrix elements, calculated nonperturbatively in numerical simulations of lattice QCD, perturbation theory has played a crucial role in the development of a complete renormalization prescription based on Ref. [33]. In the latter work the renormalization was addressed in lattice perturbation theory and a finite mixing was identified among nonlocal operators of twist-2 and twist-3. The complete mixing pattern discussed in Refs. [33] led to the proposal of a nonperturbative RI-type scheme [34,36], also employed in Ref. [37]. This development of the renormalization of nonlocal operators has been a crucial aspect in state-of-the-art numerical simulations, e.g. the work of Refs. [38,39] (for a recent overview on lattice QCD calculations see Ref. [21] and references therein).
In this work we generalize the calculation of Ref. [33] to include nonlocal operators with a staple-shaped Wilson line. We compute their Green's functions to one-loop level in perturbation theory using dimensional (DR) and lattice (LR) regularizations. The functional form of the Green's functions reveals the renormalization pattern and mixing among operators of different Dirac structure, in each regularization. We find that these operators renormalize multiplicatively in DR, but have finite mixing in LR. Results for both regularizations have been combined to extract the renormalization functions in the lattice MS scheme. In addition, the results in DR have been used to obtain the conversion factor between RI-type and MS schemes. We also present an extension to operators containing a Wilson line of arbitrary shape on the lattice, with n cusps. Preliminary results of the current work have been presented in Ref. [40].
Staple-shaped nonlocal operators (see Fig. 1) are crucial in studies of TMDs, which encode important details on the internal structure of hadrons. In particular, they give access to the intrinsic motion of partons with respect to the transverse momentum, through the formalism of QCD factorization, that can be used to link experimental data to the three-dimensional partonic structure of hadrons. An operator with a staple of infinite length, η→∞, (see Fig. 1) enters the analysis of semi-inclusive deep inelastic scattering (SIDIS) processes 3 in a kinematical region where the photon virtuality is large and the measured transverse momentum of the produced hadron is of the order of Λ QCD [41].
FIG. 1: Staple-shaped gauge links as used in analyses of SIDIS and Drell-Yan processes. For notation, see Ref. [42].
To date, only limited studies of TMDs exist in lattice QCD (see, e.g., Refs. [42][43][44][45] and references therein), such as the generalized Sivers and Boer-Mulders transverse momentum shifts for the SIDIS and Drell-Yan cases. These studies include staple links of finite length that is restricted by the spatial extent of the lattice volume. To recover the desired infinite length one checks for convergence as the length increases, and an extrapolation to η→∞ is applied. More recently, the connection between nonlocal operators with staple-shaped Wilson line and orbital angular momentum [46,47] has been discussed. This relies on a comparison between straight and staple-shaped Wilson lines, with the staple-shaped path yielding the Jaffe-Manohar [48,49] definition of quark orbital angular momentum, and the straight path yielding Ji's definition [49][50][51]. The difference between these two can be understood as the torque experienced by the struck quark as a result of final state interactions [49,50].
An important aspect of calculations in lattice QCD is the renormalization that needs to be applied on the operators under study (unless conserved currents are used). As is known from older studies [4,[6][7][8]52], the renormalization of Wilson-line operators in continuum theory (except DR) includes a divergent term e −δmL , where δm is a dimensionful quantity whose magnitude diverges linearly with the regulator, and L is the total length of the contour. For staple-shaped operators, L = (2|y|+|z|), where y ≶ 0 and z ≶ 0 define the extension of the staple in the y−z plane, chosen to be spatial. The existing lattice calculations of staple-shaped operators assume that the lattice operators have the same renormalization properties as the continuum operators, in particular that there is no mixing present. This allows one to focus on ratios between such operators [42][43][44][45] in order to cancel multiplicative renormalization, which is currently unknown 4 . However, as we show in this paper, this is not the case for operators where finite mixing is present and must be taken into account.
One of the main goals of this work is to provide important information that may impact nonperturbative studies of TMDs and potentially lead to the development of a nonperturbative renormalization prescription similar to the case of quasi-PDFs discussed above. The paper is organized in five sections including the following: In Sec. II we provide the set of operators under study, the lattice formulation, the renormalization prescription for nonlocal operators that mix under renormalization and the basics of the conversion to the MS scheme. Section III presents our main results in dimensional and lattice regularization. This includes both the renormalization functions and conversion factors between the RI and MS schemes. An extension of the work to include general nonlocal Wilson-line operators with n cusps is presented in Sec. IV, while in Sec. V we give a summary and future plans. For completeness we include two appendices where we give the expressions for the Green's functions in dimensional regularization (Appendix A) as well as the expressions related to the renormalization of the fermion fields (Appendix B).

II. CALCULATION SETUP
In this section we briefly introduce the setup of our calculation, along with the notation used in the paper. We give the definitions of the operators and the lattice actions; we also provide the renormalization prescriptions that we use in the presence of operator mixing.

A. Operator setup
The staple-shaped Wilson-line operators have the following form: where W denotes a staple with side lengths |z| and |y|, which lies in the plane specified by the directionsμ 1 andμ 2 (see Fig. 2); it is defined by The symbol Γ can be one of the following Dirac matrices: 1 1, γ 5 , γ µ , γ 5 γ µ , σ µν (where µ, ν = 1, 2, 3, 4 and σ µν = [γ µ , γ ν ]/2). For convenience, we adopt the following notation for each Dirac matrix: The fermion and antifermion fields appearing in O Γ can have different flavor indices. Operators with different flavor content cannot mix among themselves; further, for massindependent renormalization schemes, flavor-nonsinglet operators which differ only in their flavor content will have the same renormalization factors and mixing coefficients. Results for the flavor-singlet case will be identical to those for the flavor-nonsinglet case at one loop, but they will differ beyond one loop and nonperturbatively; however, the setup described below [Eqs. (7 -12)] will be identical in both cases. x x + zμ 1

B. Lattice actions
In our lattice calculation we make use of the Wilson/clover fermion action [53]. In standard notation it reads where a is the lattice spacing and Following common practice, we henceforth set the Wilson parameter r equal to 1. The clover coefficient c SW will be treated as a free parameter, for wider applicability of the results. The mass term (∼ m f 0 ) will be irrelevant in our one-loop calculations, since we will apply massindependent renormalization schemes. The above formulation, and thus our results, are also applicable to the twisted mass fermions [54] in the massless case. One should, however, keep in mind that, in going from the twisted basis to the physical basis, operator identifications are modified (e.g., the scalar density, under "maximal twist", turns into a pseudoscalar density, etc.).
For gluons, we employ a family of Symanzik improved actions [55], of the form, where U plaq. is the 4-link Wilson loop and U rect. , U chair , U paral. are the three possible independent 6-link Wilson loops (see Fig. 3). The Symanzik coefficients c i satisfy the following normalization condition: c 0 + 8c 1 + 16c 2 + 8c 3 = 1.
plaquette rectangle chair parallelogram For the numerical integration over loop momenta we selected a variety of values for c i ; for the sake of compactness, in what follows we will present only results for some of the most frequently used sets of values, corresponding to c 2 = c 3 = 0, as shown in Table I.

C. Renormalization prescription
The renormalization of nonlocal operators is a nontrivial process. As shown in our study of straight-line operators in Ref. [33], a hidden operator mixing is present in chiralitybreaking regularizations, such as the Wilson/clover fermions on the lattice. This mixing does not involve any divergent terms; it stems from finite regularization-dependent terms, which are not present in the MS renormalization scheme, as defined in dimensional regularization (DR). Thus, our first goal is to compute perturbatively all renormalization functions and mixing coefficients which arise in going from the lattice regularization (LR) to the MS scheme. Ultimately, a nonperturbative evaluation of all these quantities is desirable; to this end, and given that the very definition of MS is perturbative, we must devise an appropriate, RI -type renormalization prescription which reflects the operator mixing. We will proceed with the definition of the renormalization factors of operators, as mixing matrices, in textbook fashion. We modify our prescription in Ref. [33] to correspond to the resulting operator-mixing pairs of the present calculation, which are different from those found in the straight-line operators. The reason behind this difference is explained in detail in Sec. IV. The mixing pairs found from our calculation on the lattice are (see Sec. III B 2): , where i can be any of the three orthogonal directions to theμ 2 direction. The remaining operators do not show any mixing, and thus their renormalization factors have the typical 1 × 1 form [see Eq. (9)]. Taking into account all the above, we define the renormalization factors which relate each bare operator O Γ with the corresponding renormalized one via the following equations: where X(Y ) stands for the regularization (renormalization) scheme: X = DR, LR, . . ., Y = MS, RI , . . . As we will see, in dimensional regularization there is no operator mixing and thus the mixing matrices are diagonal.
As is standard practice, the calculation of the renormalization factors of O Γ stems from the evaluation of the corresponding one-particle-irreducible (1-PI) two-point amputated Green's functions Λ Γ ≡ ψ f O Γψf amp . According to the definitions of Eqs. (7 -9), the relations between the bare Green's functions and the renormalized ones are given by 5 where are the renormalization factors of the external quark fields of flavors f and f respectively, defined through the relation, We note that in the case of massless quarks, the flavor content does not affect the renormalization factors of fermion fields or the Green's functions of O Γ , and thus we omit the flavor index in the sequel. We also note that for regularizations which break chiral symmetry (such as Wilson/clover fermions), an additive mass renormalization is also needed, beyond one loop; however, this is irrelevant for our one-loop calculations. The expressions of Λ Γ depend on the coupling constant g 0 , whose renormalization factor is defined through where µ is related to the MS renormalization scaleμ (μ ≡ µ(4π/e γ E ) 1/2 , γ E is Euler's constant) and D is the number of Euclidean spacetime dimensions (in DR: D ≡ 4 − 2ε, in LR: D = 4). For our one-loop calculations, Z X,Y g is set to 1 (tree-level value).
There are four one-loop Feynman diagrams contributing to Λ Γ , shown in Fig   In our computations we make use of two renormalization schemes: the modified minimalsubtraction scheme (MS) and a variant of the modified regularization-independent scheme (RI ). The second one is needed for the nonperturbative evaluations of the renormalized Green's functions Λ Γ on the lattice, which will be converted to MS, through appropriate conversion factors. For our perturbative lattice calculations, the renormalization factors of O Γ in the MS scheme can be derived by calculating Eqs. (10 -12) for both X = LR and X = DR, and demanding that their left-hand sides are X-independent and, thus, identical in the two regularizations.
For the RI scheme, we extend the standard renormalization conditions for the bilinear operators, consistently with the definitions of Eqs. (10 -12), where Λ tree Γ ≡ Γ exp(iq µ 1 z) is the tree-level value of the Green's functions of O Γ ,q is the RI renormalization scale 4-vector, and N c is the number of colors. Note that the traces appearing in Eqs. (15 -17) regard only Dirac and color indices; in particular, Eqs. (15) and (16) retain their 2 × 2 matrix form, and thus they each correspond to four conditions. We mention that an alternative definition of the RI scheme can be adopted so that the renormalization factors depend only on a minimal set of parameters, (q 2 ,q µ 1 ,q µ 2 ), rather than all the individual components ofq; this can be achieved by taking the average over all allowed values of the indices i, j, in conditions (16) and (17), whenever i, j are present. This alternative scheme is not so useful in lattice simulations, where, besides the two special directions of the plane in which the staple lies, the temporal direction stands out from the remaining spatial directions; this leaves us with only one nonspecial direction, and thus this choice of normalization is not particularly advantageous in this case.
The RI renormalization factors of fermion fields can be derived by imposing the massless normalization condition, where S RI ≡ ψ RI ψ RI is the RI -renormalized quark propagator and S tree ≡ (i / q) −1 is its tree-level value.

D. Conversion to the MS scheme
The conversion of the nonperturbative RI -renormalized Green's functions Λ RI Γ to the MS scheme can be performed only perturbatively, since the definition of MS is perturbative in nature. The corresponding one-loop conversion factors between the two schemes are extracted from our calculations, and their explicit expressions are presented in Sec. III. As a consequence of the observed operator-pair mixing, some of the conversion factors will be 2 × 2 matrices, just as the renormalization factors of the operators. Following the definitions of Eqs. (7 -8), they are defined as Being regularization independent, they can be evaluated more easily in X = DR; in this regularization there is no operator mixing, and thus the conversion factors of O Γ turn out to be diagonal. We note in passing that the definition of the MS scheme depends on the prescription used for extending γ 5 to D dimensions 6 ; this, in particular, will affect conversion factors for the pseudoscalar and axial-vector operators. However, such a dependence will only appear beyond one loop.
Given that the conversion factors are diagonal, the Green's functions of O Γ in the RI scheme can be directly converted to the MS scheme through the following relation, valid for all Γ: where is the conversion factor for fermion fields.

III. CALCULATION PROCEDURE AND RESULTS
In this section we proceed with the one-loop calculation of the renormalization factors of the staple operators in the RI and MS renormalization schemes, both in dimensional and lattice regularizations. We apply the prescription described above, and we present our final results. We also include the one-loop expressions for the conversion factors between the two schemes.
A. Calculation in dimensional regularization

Methodology
We calculate the bare Green's functions of the staple operators in D Euclidean spacetime dimensions (where D ≡ 4 − 2ε and ε is the regulator), in which momentum-loop integrals are well-defined. The methodology for calculating these integrals is briefly described in our previous work regarding straight Wilson-line operators [35,61], and it is summarized below: We follow the standard procedure of introducing Feynman parameters. The momentum-loop integrals depend on exponential functions of the µ 1 -and/or µ 2 -component of the internal momentum [e.g., exp(ip µ 1 z), exp(ip µ 1 ζ)]. The integration over the components of momentum without an exponential dependence is performed using standard D-dimensional formulae (e.g., [62]), followed by a subsequent nontrivial integration over the remaining components p µ 1 and/or p µ 2 . The resulting expressions contain a number of Feynman parameter integrals and/or integrals over ζ-variables stemming from the definition of O Γ , which depend on modified Bessel functions of the second kind, K n and which do not have a closed analytic form; they are listed in Appendix A. We expand these expressions as Laurent series in ε and we keep only terms up to O(ε 0 ). The full expressions of the bare Green's functions of O Γ are given in Appendix A; these can be used for applying any renormalization scheme.

Renormalization factors
Our one-loop results for the renormalization factors of the staple operators in both MS and RI schemes are presented below.
In the MS scheme, only the pole parts [O(1/ε) terms] contribute to the renormalization factors. Diagram d 1 has no 1/ε terms, as it is finite in D = 4 dimensions. Also, it gives the same expressions with the corresponding straight-line operators, because it involves only the zero-gluon operator vertex. This statement is true in any regularization. As we expected, the divergent terms arise from the remaining diagrams d 2 − d 4 , in which end point [Eqs. (24,25)] and cusp divergences [Eq. (26)] arise. We provide below the pole parts for each subdiagram: where C F = (N 2 c −1)/(2N c ) and β is the gauge fixing parameter , defined such that β = 0 (1) corresponds to the Feynman (Landau) gauge. It is deduced that diagrams d 2 , d 3 give the same pole terms as in the case y = 0, since only end points affect these diagrams (no cusps). Also, the result for the cusp divergences of angle π/2 [Eq. (26)] agree with previous studies of nonsmooth Wilson-line operators for a general cusp angle θ [5,9,10]: it follows from these studies that the one-loop result corresponding to each of the diagrams d 4 (d) and d 4 (f ) is given by −(g 2 C F )/(16π 2 ε) (2θ cot θ + β). By imposing that the MS-renormalized Green's functions of O Γ are equal to the finite parts (exclude pole terms) of the corresponding bare Green's functions, we derive the renormalization factors of O Γ in MS, using Eqs. (10 -12); the result is given below, where we make use of the one-loop expression for the renormalization factor Z DR,MS ψ , given in Appendix B [Eq. (B1)]. Since the pole parts are multiples of the tree-level values Λ tree Γ , the nondiagonal elements of the MS renormalization factors, defined in Eqs. (7,8), are equal to zero. The diagonal elements, shown in Eq. (27), depend neither on the Dirac structure, nor on the lengths of the staple segments; further, they are gauge invariant.
In the RI scheme, there are additional finite terms, which contribute to the renormalization factors of O Γ [according to the conditions of Eqs. (15 -17)]. These terms depend on the external momentum, and they stem from all Feynman diagrams. They are also multiples of the tree-level values of the Green's functions. As a consequence, the RI mixing matrices, defined in Eqs. (7,8), are also diagonal. Therefore, there is no operator mixing in DR. The results for Z DR,RI Our resulting expressions for the conversion factors are given in the following subsection [Eqs. (29 -33)].

Conversion factors
We present below our results for the conversion factors of staple operators between the RI and MS schemes. Since the renormalization factors of O Γ are diagonal in both MS and RI schemes, the conversion factors will also be diagonal. Our expressions depend on integrals of modified Bessel functions of the second kind K n , over one Feynman parameter and possibly over one of the variables ζ appearing in Eq. (2). These integrals are denoted by P i ≡ P i (q 2 ,q µ 1 , z), Q i ≡ Q i (q 2 ,q µ 1 ,q µ 2 , z, y) and R i ≡ R i (q 2 ,q µ 1 ,q µ 2 , z, y); they are defined in Eqs. (A11 -A28) of Appendix A.
C RI ,MS C RI ,MS Tµ 1 ν = 1 + g 2 C F 16π 2 15 + 2(β + 6)γ E + 7 log µ 2 q 2 + (β + 2) log q 2 z 2 4 + 4 log q 2 y 2 4 We note that the real parts of the above expressions, as well as the bare Green's functions, are not analytic functions of z (y) near z → 0 (y → 0); in particular, the limit z → 0 leads to quadratic divergences, while the limit y → 0 leads to logarithmic divergences. The singular limits were expected, due to the appearance of contact terms beyond tree level. In the case y = 0, the staple operators are replaced by straight-line operators of length |z|, the renormalization of which is addressed in our work of Ref. [33]. In the case z = 0, the nonlocal operators are replaced by local bilinear operators, the renormalization of which is studied, e.g., in Refs. [59,60,63,64].
Since our results for the conversion factors will be combined with nonperturbative data, it is useful to employ certain values of the free parameters mostly used in simulations. To this end, we set:μ = 2 GeV and β = 1 (Landau gauge). For the RI scale we employ values which are relevant for simulations by ETMC [25], as follows: aq = ( 2π L n 1 , 2π L n 2 , 2π L n 3 , 2π T (n 4 + 1 2 )), where a is the lattice spacing, (L 3 × T ) is the lattice size and (n 1 , n 2 , n 3 , n 4 ) is a 4-vector defined on the lattice. A standard choice of values for n i is the case n 1 = n 2 = n 3 = n 4 , in which the temporal component n 4 stands out from the remaining equal spatial components. As an example we apply (n 1 , n 2 , n 3 , n 4 ) = (4,4,4,9), L = 32, T = 64 and a = 0.09 fm. For a better assessment of our results, we plot in Fig. 6 the real and imaginary parts of the quantities C Γ , defined through C RI ,MS Γ = 1 + g 2 C F 16π 2 C Γ + O(g 4 ), as functions of the dimensionless variables z/a and y/a, using the above parameter values. In the case y = 0, we use the expressions of the conversion factors for straight-line operators, calculated in Ref. [33], while in the case z = 0, we use the one-loop expressions of the conversion factors for local bilinear operators, written in Refs. [60,64]. For definiteness, we choose µ 1 = 1 and µ 2 = 2. Graphs for C V 4 = C A 4 , C T 12 = C T 13 = C T 42 = C T 34 and C T 14 = C T 32 are not included in Fig. 6, as their resulting values are very close to those of C V 2 (fractional differences: 10 −3 ).
The real parts of C Γ are even functions of both z/a and y/a. In Fig. 6, one observes that, for large values of z/a, they tend to stabilize, while for large values of y/a they tend to increase; thus, a two-loop calculation of the conversion factors is essential for more sufficiently convergent results. Further, the dependence on the choice of Γ becomes milder for increasing values of z/a and y/a. Regarding the imaginary parts of C Γ , they are odd functions of z/a and even functions of y/a. For large values of z/a or y/a, they tend to FIG. 6: Real (left panels) and imaginary (right panels) parts of the quantities C S = C P , C V 1 = C A 1 and C V 2 = C V 3 = C A 2 = C A 3 , involved in the one-loop expressions of the conversion factors: C RI ,MS Γ = 1 + g 2 C F 16π 2 C Γ + O(g 4 ), as functions of z/a and y/a [for β = 1,μ = 2 GeV , a = 0.09 fm, aq = ( 2π L n 1 , 2π L n 2 , 2π L n 3 , 2π T (n 4 + 1 2 )), L = 32, T = 64, (n 1 , n 2 , n 3 , n 4 ) = (4,4,4,9)]. Here, we choose µ 1 = 1 and µ 2 = 2.
converge to a positive value. In particular, when both z/a and y/a take large values, the imaginary parts tend to zero. For large values of y/a and, simultaneously, small values of z/a, the imaginary parts of C Γ demonstrate a small fluctuation around zero, which differs for each Γ, either in form (e.g., C V 2 and C S have opposite signs for given values of z/a, y/a) or in magnitude (e.g., the fluctuation of C S is bigger and sharper than the fluctuation of C V 1 ). As regards theq dependence, we have not included further graphs for the sake of conciseness; however, testing a variety of values for the components of aq, used in simulations, we find no significant difference, especially for large values of z/a and y/a.

Methodology
At first, let us give the lattice version of the staple operators, where upper (lower) signs of the first and third parenthesis correspond to m > 0 (m < 0) and upper (lower) signs of the second parenthesis correspond to n > 0 (n < 0). The calculation of the bare Green's functions of such nonlocal operators on the lattice is more complicated than the corresponding calculation of local operators; the products of gluon links lead to expressions whose summands, taken individually, contain possible additional IR singularities along a whole hyperplane, instead of a single point [terms ∼ 1/sin(p µ /2) or 1/sin 2 (p µ /2)]. Also, the UV-regulator limit, a → 0, is more delicate in this case, as the Green's functions depend on a through the additional combinations z/a, y/a, besides the combination aq (where q is the external quark momentum). Thus, we have to modify the standard methods of evaluating Feynman diagrams on the lattice [65], in order to apply them in the case of nonlocal operators.
The procedure that we used for the calculation of the bare Green's functions of O latt.
Γ is briefly described in our previous work regarding straight Wilson-line operators [33], and it is summarized below: The main task is to write the lattice expressions, in terms of continuum integrals, which are easier to calculate, plus lattice integrals independent of aq; however, the latter will still have a nontrivial dependence on z/a and y/a. To this end, we perform a series of additions and subtractions to the original integrands: we extend the standard procedure of Kawai et al. [65], in order to isolate the possible IR divergences stemming from the integration over the p µ component, which appears on the integrals' denominators [∼ 1/sin(p µ /2) or 1/sin 2 (p µ /2)]. To accomplish this, we add and subtract to the original integrands the lowest order of their Laurent expansion in p µ . Also, in order to end up with continuum integrals, we add and subtract the continuum counterparts of the integrands; then, the integration region can be split up into two parts: the whole domain of the real numbers minus the region outside the Brillouin zone. The above operations allow us to separate the original expressions into a sum of two parts: one part contains integrals which can be evaluated explicitly for nonzero values of a, leading to linear or logarithmic divergences, and a second part for which a naive a → 0 limit can be taken, e.g., The numerical integrations entail a very small systematic error, which is smaller than the last digit presented in all results shown in the sequel.

Green's functions and operator mixing
The results for the bare lattice Green's functions of the staple operators are presented below in terms of the MS-renormalized Green's functions, derived by the corresponding calculation in DR, where α i are numerical constants which depend on the gluon action in use; their values are given in Table II for the Wilson, Tree-level Symanzik and Iwasaki gluon actions 7 . We note that α 2 , α 3 , and α 4 have the same values (up to a sign) as the corresponding coefficients in the straight-line operators [33].  In Eqs. (37,38), we observe that there is a linear divergence [O(1/a)], which depends on the length of the staple line (|z|+2 |y|); this was expected according to the studies of closed Wilson-loop operators in regularizations other than DR [4]. This divergence arises from the tadpolelike diagram d 4 and in particular from the subdiagrams d 4 (a), d 4 (b), d 4 (c). We note that the coefficient α 2 entering the strength of the linear divergence, is given by where D(p) µν is the gluon propagator,ν is the direction parallel to each straight-line segment of the Wilson line andp equals the four-vector momentum p with p ν → 0. Moreover, additional contributions of different Dirac structures than the original operators appear ([Γ, γ µ 2 ] terms); these contributions arise from the "sail" diagrams d 2 , d 3 and in particular from the subdiagrams d 2 (c), d 3 (a). In order to obtain on the lattice the same results for the MS-renormalized Green's functions as those obtained in DR, we have to subtract such regularization dependent terms in the renormalization process. A simple multiplicative renormalization cannot eliminate these terms; the introduction of mixing matrices is therefore necessary. However, for the operators with Γ = S, V µ 2 , A i , T ij , where i = j = µ 2 = i, the contribution [Γ, γ µ 2 ] is zero, and, thus, there is no mixing for these operators. In conclusion, there is mixing between the operators , where i = µ 2 , as we have mentioned previously. This feature must be taken into account in the nonperturbative renormalization of TMDs.

Renormalization factors
The MS renormalization factors can be derived by the requirement that the terms in Eq. (38) vanish in the renormalized Green's functions. Thus, through Eqs. (10 -12), one obtains the following results for the diagonal and nondiagonal elements of the renormalization factors: (41) where the coefficients e ψ i stem from the renormalization factor of the fermion field Z LR,MS ψ , given in Appendix B [Eq. (B5)].
A number of observations are in order, regarding the above one-loop results: both diagonal and nondiagonal elements of the renormalization factors are operator independent, just as the corresponding renormalization factors in DR. Also, the dependence of the diagonal elements on the clover coefficient c SW is entirely due to the renormalization factor of fermion fields; on the contrary, the dependence of the nondiagonal elements on c SW is derived from the Green's functions of the operators, and in particular it is different for each choice of gluon action. Consequently, tuning the clover coefficient we can set the nondiagonal elements of the renormalization factors to zero and, thus, suppress the operator mixing. At one-loop level, this can be done by choosing c SW = −α 3 /α 4 . For the gluon actions given in this paper, the values of the coefficient c SW , which lead to no mixing at one loop, are 1.7442 for Wilson action, 1.6623 for tree-level Symanzik action and 1.5222 for Iwasaki action; these values are the same as those, which eliminate the mixing in the case of straight-line operators [33].
In the RI scheme, the renormalization factors can be read off our expressions for the conversion factors, given in Eqs. (29 -33), in a rather straightforward way, Since the conversion factors are diagonal, the one-loop nondiagonal elements of the RI renormalization factors are equal to the corresponding MS expressions.

IV. EXTENSION TO GENERAL WILSON-LINE LATTICE OPERATORS WITH n CUSPS
The current study of staple operators, along with our previous work on straight-line operators [33], lead us to some interesting conclusions about nonlocal operators. From these two cases, we can completely deduce the renormalization coefficients of a general Wilson-line operator with n cusps, defined on the lattice; in particular, we determine both the divergent (linear and logarithmic) and the finite parts of multiplicative renormalizations, as well as all mixing coefficients. We can also justify the nature of the mixing in each case.
The coefficients α i are numerical constants, which depend on the Symanzik coefficients of the gluon action in use; their values for Wilson, tree-level Symanzik and Iwasaki gluons are given in Tables II and III. Comparing the above results for the two types of operators, we come to the following conclusions which can be generalized to Wilson-line lattice operators of arbitrary shape: • The linear divergence [O(1/a)] depends on the Wilson line's length.
• Diagram d 1 gives a finite, regulator-independent result in all cases.
• The only contribution of sail diagrams (d 2 and d 3 ) to ∆Λ Γ comes from their end points. This is because any parts of a segment which do not include the end points will give finite contributions to Λ LR Γ , in which the naïve continuum limit a → 0 can be taken, leading to the same result as in DR and thus to a vanishing contribution in ∆Λ Γ . Consequently the shape of the Wilson line is largely irrelevant and, indeed, all numerical coefficients in Eq. (50) coincide with those in Eq. (46). The only dependence on the shape regards the Dirac structure of the operator which mixes with O Γ . The mixing terms depend on the direction of the Wilson line in the end points. For the straight Wilson line, the direction in both end points is sgn(z)μ 1 , which leads to the appearance of the additional Dirac structure sgn(z)(Γγ µ 1 +γ µ 1 Γ) upon adding together sail diagrams d 2 and d 3 . For the staple Wilson line, the direction in the left end point is sgn(y)μ 2 and in the right end point is −sgn(y)μ 2 ; thus, the additional Dirac structure which appears upon adding the two sail diagrams is sgn(y)(Γγ µ 2 − γ µ 2 Γ).
The mixing pairs for each type of nonlocal operator can be also explained (partially) by symmetry arguments. For straight-line operators, there is a residual rotational (or hypercubic, on the lattice) symmetry (including reflections) with respect to the three transverse directions to theμ direction parallel to the Wilson line. As a consequence, operators which transform in the same way under this residual symmetry can mix among themselves, under renormalization; i.e., mixing can occur only among the pairs of operators (O Γ , O Γγµ ). This argument can now be applied to a general Wilson line: given that only end points contribute, mixing can occur only with O Γγµ , whereμ refers to the directions of the two end points of the line. Clearly, the subsets of operators which finally mix depend on the commutation properties between Γ and γ µ . We note that, if the fermion action in use preserves chiral symmetry, then none of the operators will mix with each other.
• The tadpole diagram (d 4 ) for the staple operators gives, aside from the linearly divergent terms, two types of contributions: one corresponds to each of the three straightline segments [first square bracket in Eq. (51)], which is identical to the corresponding contribution from the straight Wilson line, multiplied by a factor of 3, and another contribution for each of the two cusps [second square bracket in Eq. (51)], multiplied by a factor of 2, which cannot be obtained from the study of straight-line operators.
As a consequence of the above, it follows that the difference Λ LR Γ − Λ MS Γ for a general Wilsonline operator with n cusps (and, hence, n + 1 segments), defined on the lattice, can be fully extracted from the combination of our results for the straight-line and the staple operators: the contributions of each straight-line segment and each cusp, appearing in the general operators, are obtained from Eqs. (44 -51). Therefore, without performing any new calculations, the result for the Green's functions of general Wilson-line lattice operators with n cusps, is determined below, 1, 2, 3, 4), (52) where where e Γ = e ψ 1 + 1 − 2α 5 − (n + 1)α 6 − nα 7 and e ψ i are given in Appendix B. It is worth noting that the results in Eqs. (56,57) are both gauge invariant, as was expected.

V. CONCLUSIONS AND FUTURE PLANS
In this paper, we have studied the one-loop renormalization of the nonlocal staple-shaped Wilson-line quark operators, both in dimensional regularization (DR) and on the lattice (Wilson/clover massless fermions and Symanzik-improved gluons). This is a follow-up calculation of Ref. [33], in which straight-line nonlocal operators are studied. These perturbative studies are parts of a wider community effort for investigating the renormalization of nonlocal operators employed in lattice computations of parton distributions (PDFs, GPDs, TMDs) of hadronic physics. A novel aspect of this calculation is the presence of cusps in the Wilson line included in the definition of the nonlocal operators under study, which results in the appearance of additional logarithmic divergences. Perturbative studies of such nonsmooth operators had not been carried out previously on the lattice. As in the case of the straight-line operators, certain pairs of these nonlocal operators mix under renormalization, for chirality-breaking lattice actions, such as the Wilson/clover fermion action. The path structure of each type of nonlocal operator (straight-line, staple, . . . ) leads to different mixing pairs. The results of the present study provide additional information on the renormalization of general nonlocal operators on the lattice.

Particular novel outcomes of our calculation are
• The one-loop results for the amputated two-point one-particle-irreducible (1-PI) Green's functions of the staple operators both in DR [Eqs. (A1 -A10)] and on the lattice [Eqs. (37,38)].
• The mixing pairs of the staple operators: , i =µ 2 (for notation, see Sec. II A). We propose a minimal RI -like condition [Eqs. (15 -17)], which disentangles this mixing and which is appropriate for nonperturbative calculations of parton-distribution functions on the lattice.
• The one-loop expressions for the renormalization factors of the staple operators in both dimensional and lattice regularizations, in the MS scheme and the proposed RI scheme [Eqs. (27, 28, 40 -43)].
• An extension of our calculations to general Wilson-line lattice operators with n cusps; we have provided results for their renormalization factors [Eqs. (56,57)].
Our results are useful for improving the nonperturbative investigations of transverse momentum-dependent distribution functions (TMDs) on the lattice. Such an example is the calculation of the generalized g 1T worm-gear shift in the TMD limit (|η|→∞); this quantity involves a ratio between the axial and vector operators. A recent study of TMDs on the lattice [45] reveals tension between results for g 1T in the clover and domain-wall formulations. This is not observed in other structures and is an indication of nonmultiplicative renormalization. Our proposed RI -type scheme can be applied to the nonperturbative evaluation of renormalization factors and mixing coefficients of the unpolarized, helicity and transversity quasi-TMDs; this is expected to fix the inconsistency between the two calculations of g 1T . Also, our one-loop conversion factors can be used to convert the RI nonperturbative results to the MS scheme. Our results for general Wilson-line lattice operators with n cusps can be used in the nonperturbative renormalization of more general continuum nonlocal operators.
Comparing our results for the staple operators with the corresponding ones for the straight-line operators, we deduce that the strength of the linear divergences is the same for both types of operators; the presence of cusps lead to additional logarithmic divergences in the staple operators. Also, the observed mixing pairs among operators with different Dirac structures depend on the direction of Wilson line in the end points, and thus, they are different between the two types of operators: the straight-line operator O Γ mixes with O {Γ,γµ 1 } , while the staple operator O Γ mixes with O [Γ,γµ 2 ] (for notation, see Sec. II A). However, the values of the mixing coefficients are the same in the two cases.
Further perturbative investigations of the staple operators can lead to improved and more robust results. Our future plans include three extensions of the present calculation: • The first one is the one-loop evaluation of lattice artifacts to all orders in the lattice spacing, for a range of numerical values of the external quark momentum, of the momentum renormalization scales, and of the action parameters, which are mostly used in simulations. Such a procedure has been successfully employed to local operators [67][68][69]. The subtraction of the unwanted contributions of the finite lattice spacing from the nonperturbative estimates is essential in order to reduce large cutoff effects in the renormalized Green's functions of the operators and to guarantee a rapid convergence to the continuum limit.
• Secondly, we intend to add stout smearing on gluon links appearing in the definition of the staple operators and to investigate its impact to the elimination of ultraviolet (UV) divergences and of operator mixing; modern simulations employ such smearing techniques for more convergent results.
• Thirdly, a natural continuation of the present work is the two-loop calculation of the conversion factors between the RI and MS schemes; higher-loop corrections will eliminate large truncation effects from the nonperturbative results. Based on our extensive studies for systematic uncertainties on the renormalization functions for the straight Wilson line [34,70], we find empirically that the one-loop conversion factor is sufficient for lattice spacing satisfying |z|/a ≤ 7−8 and (aµ) 2 within the interval [2 − 4]. Outside these regions, a two-loop conversion factor would be called for; clearly, however, other systematic uncertainties will also become more relevant (lattice artifacts, volume effects, etc).
Finally, our perturbative analysis can be also applied to the study of further composite Wilson-line operators, relevant to different quasidistribution functions, e.g., gluon quasi-PDFs, etc.

List of integrals:
In what follows, we use the notation: s ≡ q 2 (1 − x)x.