Renormalization of asymmetric staple-shaped Wilson-line operators in lattice and continuum perturbation theory

In this work, we study the renormalization of nonlocal quark bilinear operators containing an asymmetric staple-shaped Wilson line at the one-loop level in both lattice and continuum perturbation theory. These operators enter the first-principle calculation of transverse momentum-dependent parton distribution functions (TMDPDFs) in lattice QCD using the formulation of Large Momentum Effective Theory. We provide appropriate RI$'$-type conditions that address the power and logarithmic divergences, as well as the mixing among staple operators of different Dirac structures, using a number of different possible projectors. A variant of RI$'$, including calculations of rectangular Wilson loops, which cancel the pinch-pole singularities of the staple operators at infinite length and reduce residual power divergences, is also employed. We calculate at one-loop order the conversion matrix, which relates the quasi-TMDPDFs in the RI$'$-type schemes to the reference scheme $\overline{\rm MS}$ for arbitrary values of the renormalization momentum scale and of the dimensions of the staple.


I. INTRODUCTION
One of the directions of research in lattice QCD, which has shown rapid progress in the last decade, is the firstprinciple study of a family of nonperturbative distribution functions (DFs) that describe the internal structure of hadrons: parton distribution functions (PDFs), generalized parton distribution functions (GPDs) and transversemomentum dependent parton distribution functions (TMDPDFs).All three types of DFs are crucial for a comprehensive understanding of the three-dimensional hadron picture.Calculating DFs from first principles has long been a challenge in Hadron Physics due to their nonperturbative and light-cone nature.The latter does not allow for a direct nonperturbative computation of DFs on a Euclidean lattice.In the last decade, a groundbreaking approach by X. Ji [1] has overcome this issue.The approach connects Euclidean equal-time correlation functions (referred to as quasi-DFs), which are accessible by lattice simulations, to the physical light-cone DFs using the framework of Large Momentum Effective Theory (LaMET) [2].This breakthrough has paved the way for extracting, for the first time, DFs from lattice simulations.
The computation of DFs using Ji's approach (also called the "quasi-PDFs" approach) involves a three-step process.Firstly, nonperturbative calculations are performed to determine hadron matrix elements of gauge-invariant quark or gluon nonlocal operators.These operators contain path-ordered Wilson lines with specific shapes, such as straight lines for PDFs and staple-shaped lines for TMDPDFs.Secondly, the nonlocal operators are renormalized to establish a connection with physically measurable quantities.This task is challenging compared to the case of local operators, as explained below.Lastly, the renormalized lattice DFs are perturbatively matched to the corresponding physical light-cone DFs using LaMET.One-loop matching formulae have been extracted in the literature for several quasi DFs [91,92]; in some cases a two-loop formula is also available [93,94].
In this study, we focus on the implementation of the second step in the quasi-PDFs approach regarding the renormalization of nonlocal operators using one-loop perturbation theory in both continuum and lattice regularizations.While a nonperturbative calculation of the renormalization functions in numerical simulations of lattice QCD is desirable, perturbation theory can give important feedback for the development of an appropriate nonperturbative renormalization prescription, which addresses all kinds of divergences, as well as possible mixing with operators of equal or lower dimension allowed by global symmetries.Most importantly, perturbation theory gives us the matching functions (at a given order) between nonperturbative renormalization schemes and continuum perturbative schemes -primarily MS -employed in phenomenology.
Studies of nonlocal operators with Wilson lines in continuum perturbation theory go back decades, including seminal work [95][96][97][98][99][100][101][102][103][104] for the renormalization of open and closed (loops) Wilson lines, with and without singular points (cusps, self-intersections), and having quark or gluon fields at the endpoints.Lattice studies of nonlocal operators have emerged only in the last decade after the development of the quasi-PDFs approach.The first perturbative lattice calculation of Wilson-line operators was made by our group in Ref. [105], to one loop for massless quarks, using the Wilson/clover fermion action and a variety of Symanzik-improved gluon actions.A straight Wilson line with quark fields at the endpoints was employed in order to investigate the renormalization of quark quasi-PDFs.This study showed that the lattice formulation introduces several new complications, such as mixing among operators of equal dimension and different twists, and power-law divergences (even in the absence of mixing with lower-dimensional operators, in contrast to the case of local operators).A number of extensions of this study have been followed by our group regarding the presence of finite quark mass [106] and the calculation of one-loop artifacts to all orders in the lattice spacing a [107].In a different extension of our study, we have considered nonlocal operators with a symmetric staple-shaped Wilson line [108] in order to investigate the renormalization of quark quasi-TMDPDFs.In our current work, we extend further the latter calculation by employing the more general case of an asymmetric staple-shaped Wilson line, where all three segments of the staple can have different lengths.Studies by other groups along these lines have appeared in Refs.[109][110][111][112].
Staple-shaped nonlocal operators have a wide range of applicability.They enter the analysis of semi-inclusive deep inelastic scattering (SIDIS) processes, as well as the Drell-Yan (DY) processes, in a kinematic region where the photon virtuality is large and the measured transverse momentum of the produced hadron is of the order of Λ QCD [113].In these analyses, the segments of the staple that are parallel to the direction in which the hadron is boosted have an infinite length.Thus, while our study focuses on finite lengths of the staple segments, an extrapolation to infinite limit must be taken in the renormalized matrix elements of these operators.The presence or not of an asymmetry in the shape of the staple operators affects their renormalization.Besides the presence of an additional scale in the renormalization conditions, the mixing pattern is different between symmetric and asymmetric staple operators (see Sec. II.2).However, this difference is not visible in one-loop lattice perturbation theory, as concluded by the present calculation.
As stated before, there are a number of challenges to address in order to renormalize the nonlocal Wilson-line operators of arbitrary shape: 1. Power divergences arise for cutoff regularized theories, such as lattice QCD [98].The divergences depend on the total length (L) of the Wilson line and can be absorbed in an exponential factor of the form e −δmL , where δm is a dimensionful regularization-dependent quantity whose magnitude diverges linearly with the regulator.
2. Logarithmic divergences arise not only from contact terms but also from the singular points of the Wilson line.
In the case of a straight line, singular points are just the endpoints.For operators of different shapes, there can also be cusp divergences [100], which depend on the angle and number of cusps present.In the case of the staple line, there are two cusps of angle π/2.These divergences can be addressed by using typical renormalization schemes, such as RI ′ or ratio schemes.
3. Finite operator mixing arises between Wilson-line operators O Γ with different (products of) Dirac matrices Γ depending on the regularization.In the case of straight-line operators, the one-loop computation [105] shows mixing in pairs between nonlocal fermion bilinears of the form (O Γ , O {Γ,γν 1 }/2 ), where ν1 is the direction of the straight line, for chirality-breaking actions.Symmetry arguments (reflections, charge conjugation) [114] confirm that the one-loop mixing pattern is also valid nonperturbatively.In the case of (symmetric and asymmetric) staple operators, the one-loop computation [108] shows mixing in pairs with a different pattern: ), where ν2 is the direction of the sides of the staple line.However, symmetry arguments [68] show a wider mixing First, we define the operators under study along with our conventions2 .The asymmetric staple-shaped Wilson-line fermion bilinear operators are defined in Euclidean space as follows: where W(x, z, y, y ′ ) denotes the staple-shaped Wilson line as given schematically in Fig. 1 and defined as: U(r, ℓ μ) denotes the straight-line path-ordered (P) exponential (Wilson line), expressed in terms of the gluon field A µ , which connects the points r and r + ℓμ.A lattice discretization of U(r, ℓ μ) in terms of gluon links U µ (x) which connect points x and x + aμ is given below: where U −µ (x) ≡ U † µ (x−aμ) and upper (lower) signs correspond to ℓ > 0 (ℓ < 0).Other discretizations involve smeared gluon links, e.g., stout, HYP, Wilson flow.We plan to investigate the impact of stout smearing at the one-loop level in future work.There is a total of 16 different operators which can be extracted from Eq. (II.1.1)depending on the choice of the Dirac matrix Γ inserted between the fermion and antifermion fields: Γ = 1 1 (scalar S), γ 5 ≡ γ 1 γ 2 γ 3 γ 4 (pseudoscalar P ), γ µ (vector V µ ), γ 5 γ µ (axial vector A µ ), σ µν ≡ [γ µ , γ ν ]/2 (tensor T µν ).Of particular interest is the study of vector, axial-vector, and tensor operators, which correspond to the unpolarized, helicity, and transversity types of TMDPDFs, respectively.
FIG. 1: The shape of the asymmetric staple line, as defined in the operator O Γ (x, z, y, y ′ ).Here, the staple is placed at x = 0.
Since we are working in one-loop perturbation theory with massless fermions, the flavor content of the fermion fields ψ and ψ is irrelevant for our calculations up to this order.Thus, our one-loop results are valid for both flavor singlet and nonsinglet operators, given that there is no additional mixing with gluon nonlocal operators [70,119].

II.2. Symmetry properties of staple-shaped operators
We examine below the properties of staple operators under symmetry transformations.Similar investigations can be found in [65,68].Operators with the same behavior under all symmetries can mix among themselves.In this way, one can identify the mixing pattern of the staple operators based on nonperturbative arguments.Since the staple operators depend on two special directions (ν 1 and ν2 ) in which the staple is defined, we consider appropriate versions of the symmetry transformations of the QCD action (in Euclidean space) with respect to the special directions.2. Two-dimensional (2D) rotational (or square, on the lattice) symmetry: The operators are covariant under rotations over the two-dimensional (2D) plane transverse to the plane in which the staple is defined.The 16 staple operators are classified into two representations of the 2D rotational group: (1) scalar: where (ν 1 , ν 2 , ν 3 , ν 4 ) correspond to different orthogonal directions in the 4D Euclidean space.The first representation is one-dimensional while the second is two-dimensional reducible (e.g., (V ν3 , V ν4 ) splits into two different one-dimensional representations: V ν3 + iV ν4 and V ν3 − iV ν4 ).Operator mixing can only occur among operators that support the same irreducible representation.The fact that an operator such as V ν3 + iV ν4 can mix with A ν3 + iA ν4 but not with A ν3 − iA ν4 implies certain relations among the corresponding renormalization and mixing coefficients.Thus, in contrast to the local operators, which are covariant under 4D rotations, the residual rotational symmetry cannot completely prevent mixing between scalar, pseudoscalar, vector, axial-vector, and tensor operators.Also, operators which support the same representation of the 4D rotational group but different representations of the 2D rotational group, e.g.V ν1 (or V ν2 ) and V ρ̸ =(ν1,ν2) , will not share the same renormalization factor (in contrast to the case of local operators); at least, the 2D rotational symmetry prevents the mixing among these operators.

Parity (P):
In Euclidean space, temporal and spatial directions are not distinguished.Thus, parity can be generalized in any direction.The generalized parity transformations P µ for the fermion and gluon fields with respect to the direction µ are defined below [124].Here, ⃗ x is the 3-vector, which is perpendicular to the µ direction.
The transformation of the staple-shaped operators under generalized parity is: In the above relation it is understood that, after the parity transformation P µ , there follows a translation T 2 ⃗ x by an amount 2 ⃗ x ; such a translation is clearly allowed by translational invariance, and it also does not affect the "relative coordinates" y, y ′ , z. Thus:

Time reversal (T ):
As in the case of parity, time reversal in Euclidean space is generalized in any direction.
The generalized time-reversal transformations T µ for the fermion and gluon fields with respect to the direction µ are defined below.We use the same notation as in parity. (II.2.9) The transformation of the staple-shaped operators under generalized time reversal is (a translation T 2 xµ has also been applied): Time reversal does not provide any additional information to the mixing compared to the residual rotational symmetry and generalized parity.
We note that all-order perturbative studies in the continuum show that a multiplicative renormalization can address all the divergences of the nonlocal Wilson-line operators [98,100,103].Thus, the mixing which is allowed by symmetries is absent in minimal subtraction schemes.However, it can occur as finite mixing in nonminimal schemes.On the lattice, the mixing is present when employing MS or any nonperturbative intermediate scheme; while MS is a minimal continuum scheme, on the lattice, additional finite regularization-dependent terms contribute to the renormalization functions in order to be able to match the continuum MS-renormalized Green's functions.In our work, we focus on the wider mixing pattern arising in the case of asymmetric staple operators when a chirality-breaking fermion action is employed.Our goal is to construct a common prescription for renormalizing the staple operators in any regularization (regularization independent).To this end, we consider the following quadruplets of the sixteen independent bilinear operators O Γ (with Γ ∈ S i ) 4 : In the present study, we do not consider possible mixing (on the lattice) with higher dimensional operators multiplied by the appropriate power of the lattice spacing.In the case of local operators, such mixing is present only for finite values of the lattice spacing, and it vanishes when taking the continuum limit.However, in the case of nonlocal operators, where power divergences a −n , n ∈ Z + are present, O(a) effects in the bare Green's functions can contribute to the renormalized Green's functions at two loops 5 .Alternatively, one can suppress these unwanted effects in two ways: (1) by removing power divergences from the Green's function through an appropriate ratio with another Green's function which has the same power divergences, e.g., a closed Wilson loop, (2) by subtracting artifacts from the bare Green's functions calculated in lattice perturbation theory.Our group has successfully applied this method to the renormalization of local quark bilinear operators [125][126][127], and more recently to the renormalization of nonlocal straight Wilson-line operators for quasi-PDFs [107].We plan to study one-loop discretization effects for the staple operators to all orders in the lattice spacing in future work.

II.3. Green's functions of staple-shaped operators with external fermions
As is standard practice, one-particle-irreducible (1-PI) two-point amputated Green's functions of the operators under study with external elementary fields, e.g., fermion fields, can be used for the extraction of renormalization functions6 : A summation over the position of the staple-shaped operator is taken, which is allowed by translational symmetry, in order to simplify the calculations.Such Green's functions using local operators are easily calculated in continuum perturbation theory to very high order.However, due to the nonlocal nature of the staple-shaped operators, additional scales (staple lengths) appear in Green's functions, which make the computation more complex even at the one-loop level.The corresponding calculation on the lattice is even more demanding since the procedure for isolating divergences from the Feynman integrals, as well as the procedure for taking the continuum limit a → 0 (where a is the lattice spacing) are more complicated (see, e.g., Refs.[105,108]).
There are four one-loop Feynman diagrams contributing to Λ Γ (q, z, y, y ′ ), shown in Fig. (2).These diagrams will appear in both continuum and lattice regularizations since all vertices are present in both regularizations.To this perturbative order, zero (d 1 ), one (d 2 − d 3 ) or two (d 4 ) gluons stem from the staple-shaped Wilson line.Due to the shape of the staple, the diagrams are further divided into thirteen subdiagrams, shown in Fig. (3), depending on the side of the staple from which gluons emanate.The transformation properties of Λ Γ (q, z, y, y ′ ) under C, P, T are given below (⃗ q is the momentum 3-vector which is perpendicular to the µ direction appearing in each transformation): By combining P µ , T µ and C, Λ Γ (q, z, y, y ′ ) transforms to: (II.3.5)

II.4. Renormalization conditions
We formulate below different versions of appropriate regularization-independent (RI ′ ) prescriptions which address all possible divergences and mixing of the staple-shaped operators.In contrast to our previous study regarding the FIG. 3: One-loop subdiagrams contributing to the Green's functions of the asymmetric staple-shaped operator with external fermions.The straight (wavy) lines represent fermions (gluons).The operator insertion is denoted by an asymmetric staple-shaped line.
renormalization of symmetric staple operators, we construct renormalization matrices that address the mixing as observed by studying symmetries and not by studying one-loop lattice perturbation theory; thus, the renormalization matrices will be 4 × 4. Several of the nondiagonal elements of these matrices will vanish at one-loop.In this way, the renormalization prescription will be more appropriate for nonperturbative calculations addressing possible mixing that can be seen in higher loops.We first give our conventions regarding the renormalization of the operators under study, as well as of the relevant elementary fields and parameters that enter our perturbative calculations: where ψ X (ψ R ) is the bare (renormalized) fermion field, g X (g R ) is the bare (renormalized) coupling constant.X denotes dimensional (DR) or lattice (LR) regularization, R denotes renormalization schemes of RI ′ or MS, µ 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) 7 .A sum over Γ ′ matrices, which belong to the mixing set of Γ (cf.Eqs.(II.2.31) -(II.2.34)), is implicit in Eq. (II.4.1)8 .Given the above conventions, the renormalized Green's function Λ R Γ (q, z, y, y ′ ) under study is defined through: • RI ′ conditions: We first employ the typical RI ′ scheme as defined for the renormalization of local fermion bilinear operators [128] by extending the renormalization conditions consistently with the mixing and the definition in Eq. (II.4.2): where ⟨ψ X (q) ψX (q)⟩ is the fermion propagator, q is the RI ′ renormalization 4-vector scale and N c is the number of colors.Note that the traces appearing in Eqs.(II.4.3 -II.4.4) regard both Dirac and color indices.Also, Eq. (II.4.3) corresponds to 16 × 4 = 64 conditions, which determine all the elements of the 4 × 4 renormalization matrices for the 4 mixing sets S i .We use two different choices of projectors in Eq. (II.4.3): where ⃗ r ≡ z ν1 + (y − y ′ ) ν2 , ⃗ q L ≡ q ν1 ν1 + q ν2 ν2 and ⃗ q T ≡ ⃗ q − ⃗ q L = q ν3 ν3 + q ν4 ν4 .Use of these two choices of projectors amounts to the implementation of two different RI ′ prescriptions, which we will denote as RI ′ 1 , RI ′ 2 from now on.Compared to the first choice of projectors, the second one can further remove finite contributions of some Dirac structures, allowed by Lorentz symmetry, from the elements of the renormalization matrices.Similar projectors have also been studied in the renormalization of local operators, leading to reduced contributions from hadronic contamination in the nonperturbative data, especially for small values of the renormalization scale q.
The RI ′ scheme can address the power9 and logarithmic divergences, as well as the mixing between different Dirac structures, in the same way as MS does.However, in contrast to the MS scheme, RI ′ can also treat the pinch-pole singularity when y → ∞.This means that this infinite limit can be taken only in the RI ′ -renormalized Green's functions and not in the MS-renormalized Green's functions of the staple operators.This is true because both schemes are defined for finite values of y, where pinch-pole singularities are not present.Then, terms that diverge with y (in the y → ∞ limit) do not contribute to the MS renormalization function, but they do so in RI ′ .Thus, when multiplying the bare Green's functions with their renormalization functions, only in RI ′ these terms are eliminated.However, since RI ′ is just an intermediate scheme entering the procedure of renormalizing the operators in the reference scheme MS, nonperturbatively, we cannot benefit from this additional feature of RI ′ ; after conversion to the MS scheme, the pinch-pole singularity comes back.Thus, the standard MS prescription is not appropriate for renormalizing staple operators in the infinite limit y → ∞.
• Alternative prescriptions: Alternative prescriptions which are also applicable in the limit y → ∞ are described below.The main feature in these prescriptions stems from the fact that the same pinch-pole singularity arises in a closed Wilson loop.Also, this object shows power and cusp divergences as the staple-shaped operators.As proposed in Ref. [65], we can redefine the standard staple-shaped operator O Γ by dividing it with the square root of the vacuum expectation value of a rectangular z × (y + y ′ ) Wilson loop L(z, y + y ′ ): where The dimensions of the Wilson loop are chosen in such a way as to cancel the power10 , cusp and pinch-pole divergences of O Γ .Both the standard MS prescription and the RI ′ prescription (Eq.II.4.3) can now be applied to the operator O Γ .In particular, the condition for RI ′ (we will call it RI ′ -bar; the related renormalization functions are denoted by Z R,X ΓΓ ′ ) now reads: where ΓΓ ′ addresses the remaining end-point divergences, as well as the mixing.
In order to extract Z R,X ΓΓ ′ , we need to calculate the one-loop Feynman diagrams of Fig. 4 contributing to ⟨L X (z, y+y ′ )⟩.Older studies of Wilson loops can be found, e.g., in Refs.[123,129] in both continuum and lattice using Wilson gluons.Here, we extend these calculations to the case of Symanzik-improved gluons.For completeness, we have also repeated the continuum calculation.
Since the end-point divergences do not depend on the dimensions of the staple11 , a nonperturbative determination of Z R,X ΓΓ ′ is expected to exhibit a much milder dependence on the staple lengths z, y, y ′ which lie in the renormalization window: In this way it becomes more acceptable to renormalize the modified operators O Γ defined at large values of the lengths z, y, y ′ , using renormalization functions Z RI ′ ,X ΓΓ ′ defined at smaller values of z, y, y ′ within the perturbative region.To distinguish the lengths appearing in the bare operators from the reference lengths appearing in the renormalization functions, we will call the latter as z, ȳ, ȳ′ .If one is interested in taking the limit y → ∞, then there are only two relevant lengths: z and ȳ − ȳ′ .
Another option that we do not consider in this work is the short-distance ratio (SDR) scheme, described in Ref. [66].In this scheme, ratios of hadron matrix elements of the modified operators at different external momenta are considered.All the good features of RI ′ -bar are also valid in the SDR scheme.However, we stress that SDR is valid when operator mixing is absent or negligible.

III.1. Green's functions:
We provide below our one-loop results for the bare Green's functions Λ DR Γ (q, z, y, y ′ ) in dimensional regularization (DR) up to O(ε 0 ).For the calculation of the one-loop momentum integrals, we make use of the formulae given in Appendix A. The results depend on several Dirac structures allowed by the residual rotational symmetry, and they are expressed in terms of integrals over Feynman parameters and/or over ζ-variables running over the sides of the staples: F i ≡ F i (q, r), G i ≡ G i (q, y, z), Ḡi ≡ Ḡi (q, y ′ , z), H i ≡ H i (q, y, z), Hi ≡ Hi (q, y ′ , z), I i ≡ I i (q, y − y ′ , z).These integrals are listed explicitly in Appendix B. The following notation is employed: Λ DR σµν 1 (q, z, y, y ′ ) = Σ 6,1 σ µν1 + Σ 6,2 σ µν2 + Σ 6,3 q µ 1 1 + Σ 6,4 q µ σ ν1ν2 + Σ 6,5 γ µ/ q + Σ 6,6 ε ν1ν2µρ γ 5 γ ρ/ q + Σ 6,7 q µ γ ν1 / q + Σ 6,8 q µ γ ν2 / q, (III.1.6) + Σ 7,5 γ µ/ q + Σ 7,6 ε ν1ν2µρ γ 5 γ ρ/ q + Σ 7,7 q µ γ ν1 / q + Σ 7,8 q µ γ ν2 / q, (III.1.7) Λ DR γ5γν 2 (q, z, y, y ′ ) = γ 5 Λ DR γν 2 (q, z, y, y ′ ), Λ DR γ5γµ (q, z, y, y ′ ) = γ 5 Λ DR γµ (q, z, y, y ′ ), (III.1.9) Λ DR σν 4 ν 3 (q, z, y, y ′ ) = γ 5 Λ DR σν 1 ν 2 (q, z, y, y ′ ), (III.1.10) where Σ i,j ≡ Σ i,j (μ 2 , q ν1 , q ν2 , q 2 , z, y, y ′ ): and s 0 ≡ s 0 (μ 2 , q 2 , z, y, y ′ ), s i,j ≡ s i,j (q ν1 , q ν2 , q 2 , z, y, y ′ ): .Linear divergences are not present in DR.Also, the term [y/z tan −1 (y/z) + y ′ /z tan −1 (y ′ /z)] appearing in s 0 (Eq.(III.1.12))gives rise to a pinch-pole singularity (linear divergence) in the limit y → ∞ for fixed values of (y − y ′ ).This term stems from the sub-diagram d 4e .Our results agree with previous studies in Refs.[108,123] which consider specific cases of the staple operators.More specifically, when y ′ → y, in which the asymmetric staple gauge link becomes symmetric, the above expressions reproduce the results of Ref. [108], where a slightly different basis of Feynman-parameter integrals has been employed, which is related to our basis (when y ′ = y) through linear combinations and/or integration by parts.Other vanishing limits of the lengths for one or more staple segments, which classically result to simpler shapes, e.g., straight-line gauge link of length y − y ′ when (z → 0), or, straight-line gauge link of length z when (y = y ′ and y → 0), or, single point when (y = y ′ and y → 0 and z → 0), are singular because cusp and endpoints do not vanish smoothly in these limits giving rise to linear or logarithmic divergences (see Eqs. (III.3.34 -III.3.35)).Thus, in these cases a complete consistency check is not applicable.However, by excluding the divergent parts of our results (s 0 ), the aforementioned limits can be taken in order to reproduce the finite parts of the corresponding Green's functions with simpler shapes of Wilson line [105].In particular, some of the form factors Σ i,j are zero resulting to fewer Lorentz structures for each operator, as expected by symmetries.The basis of Feynman-parameter integrals is now reduced due to vanishing terms, or identical integrals.The different Dirac structures appearing in Eqs.(III.1.1 -III.1.10)can be classified into two categories: (1) the structures corresponding to the tree-level Green's functions of the staple operators that belong to the same mixing set according to symmetries (structures multiplied by Σ i,1 or Σ i,2 ), and (2) all the remaining structures.The former structures are multiplied by regularization-dependent form factors (Σ i,j ), while the form factors of the latter structures do not depend on regularization.This is confirmed by comparing our one-loop results in both continuum and lattice regularizations (Eq.(IV.1.2)).In particular, the bare Green's functions in the continuum consist of two instead of four structures of the first category because chiral symmetry is preserved (see discussion in Sec.II.2).On the other hand, the lattice bare Green's functions at one loop also contain a pair of the four first-category structures, which is different from that obtained in DR (see Sec. IV).Since the lattice regularization that we employ breaks chiral symmetry, we expect that the missing two structures of the first category will appear in higher loops.
By taking into account the classification of the Dirac structures in the Green's functions, we conclude that a renormalization prescription that considers a wider mixing pattern, i.e., mixing among all 16 independent operators of different Γ matrices, is not an optimal choice for renormalizing the staple operators.In this case, several form factors Σ i,j will contribute to the 16 × 16 mixing matrix, including pure dependence on the regularization-independent form factors Σ i,j>2 in some nondiagonal elements.Given that we are interested in matching bare lattice Green's functions to the MS scheme (through an intermediate nonperturbative scheme), in which the regularization-independent form factors cannot contribute, it is more economic to consider an intermediate renormalization prescription that includes only the minimum number of operators necessary for disentangling the mixing occurring in the MS scheme on the lattice.In our study, we consider such a minimal set of operators as dictated by symmetries.
The one-loop renormalized Green's functions in the MS scheme can be obtained by removing the O(1/ε) term in Eq. (III.1.11).Also, the renormalized Green's functions in any variant of the RI ′ scheme can be obtained by imposing the corresponding renormalization conditions to the above results.The extraction of the conversion functions among these schemes (RI ′ variants and MS) is then straightforward; they are given in Eqs.(III.3.1,III.3.4 -III.3.32) for the schemes under study.
Specific combinations of the form factors Σ i,j contribute in the renormalization conditions of Eq. (II.4.3) for each choice of projectors P [1] Γ and P [2] Γ .The latter choice has the advantage that in most cases, only Σ i,1 and Σ i,2 survive, which multiply the relevant structures for the study of mixing.However, in the case of projectors P [2] γµ and P [2] γ5γµ , for µ = ν 3 , ν 4 , some additional form factors can survive when taking the trace with the amputated Green's functions Λ γ5γρ and Λ γρ , respectively (where ρ ̸ = (ν 1 , ν 2 , µ)).For instance, Eq. (II.4.3) for the case Γ = Γ ′′ = γ ν3 takes the following form (in terms of Σ i,j ), when employing P γν 3 gives: We see that the second choice of the projector can give a simpler combination of Σ i,j , at least for the term multiplying Z DR,RI ′ γν 3 ,γν 3 .But here, the main advantage of using the second choice of projectors is that the renormalization matrices are independent of the individual components of the momentum scale, which are orthogonal to the plane in which the staple lies.Eq. (III.1.42)depends explicitly on qν3 , while Eq.(III.1.43)only depends on qν3 through q2 T ≡ q2 ν3 + q2 ν4 .Indeed, the same dependence on transverse components of q regards all renormalization functions resulting from the choice of projectors P [2] Γ .Consequently, the nonperturbative determination of renormalization functions in this prescription will only depend on q2 T , rather than on its individual components.This gives us the possibility of increasing statistics in the nonperturbative calculations of renormalization functions on the lattice by averaging over different qν3 and qν4 components which have the same q2 T , according to the residual 2D symmetry.Furthermore, we expect that by eliminating contributions coming from the form factors Σ i,j (j > 2), hadronic contaminations, which may be present in the nonperturbative calculation of the renormalization functions, will be reduced, as it happens in the case of local quark bilinear operators.

III.2. Renormalization functions
The renormalization functions of O Γ and O Γ in the MS scheme have been determined by imposing that the MS-renormalized Green's functions of the two operators are equal to the finite parts (exclude pole terms) of the corresponding bare Green's functions: As expected, the renormalization functions in Eqs.(III.2.1 -III.2.2) are diagonal (there is no mixing) and independent of the Dirac matrices Γ, Γ ′ , and the lengths of the staple.In addition, the result coincides with that obtained for the symmetric staple operators in our previous work [108]; thus, the asymmetry in the shape of the staple does not affect the renormalization in MS.The divergent term 7/ε of Eq. (III.2.1) comes from the sum of the pole terms in the bare Green's function (III.1.11) and in the renormalization factor of the external quark fields (see Appendix C).The divergent term 3/ε of Eq. (III.2.2) includes the additional O(1/ε) terms coming from the vacuum expectation value of the Wilson loop (Eq.III.1.44).In the latter case the resulting 3/ε term coincides with that found in the renormalization of a straight Wilson-line operator.
In the RI ′ -type schemes (RI ′ 1 , RI ′ 2 , RI ′ 1 -bar, RI ′ 2 -bar), the renormalization functions are (in general) nondiagonal matrices due to the operator mixing.Their expressions are given below in terms of the conversion matrices that connect the RI ′ -type schemes with the MS scheme, as defined in the next subsection (Eq.(III.3.1)): 4 ). (III.2.4)

III.3. Conversion matrices
By using our results for the bare Green's functions, we extract the conversion matrices that match the renormalization matrices of the staple operators from RI ′ -type (RI ′ 1 , RI ′ 2 , RI ′ 1 -bar, RI ′ 2 -bar) schemes to MS: The conversion matrices are regularization (X) independent, and thus, our results in DR are also applicable to the lattice.They can be used in lattice simulations in order to translate the renormalized hadron matrix elements of the staple operators, such as the quasi-beam function [68] and the soft function [62], from the intermediate RI ′ -type schemes to MS.According to the mixing pattern, we need to calculate a total of 4 (one for each mixing set S i ) conversion matrices of dimensions 4 × 4 for each RI ′ -type scheme.We found that the matrices are block diagonal in 2 sub-matrices of dimensions 2 × 2 at one-loop level 13 ; these blocks have the same generic structure for all Γ: Introducing a γ 5 matrix in both operator and projector leaves the conversion matrix unaffected to one loop 14 : Thus, the conversion matrix for the mixing set S 2 is identical to the conversion matrix for the mixing set S 1 .Also the conversion matrices for the mixing sets S 3 and S 4 are related through the interchange ν 3 ⇆ ν 4 15 .We provide below all nonzero elements of the conversion matrices for each mixing set for the RI ′ 1 and RI ′ 2 schemes; they are expressed in terms of the coefficients s 0 and s i,j given in Eqs.(III.1.12-III.1.41)by setting q = q.Here we use the notation ε LC ≡ ε ν1ν2ν3ν4 , and µ = ν 3 , ν 4 .
(III.3.33) As observed, the conversion matrices in the RI ′ i and RI ′ i -bar schemes exhibit differences in their diagonal elements while remaining unchanged in the nondiagonal elements up to one loop.
There is a nontrivial dependence of the conversion matrices (as well as of the bare Green's functions) on the staple lengths (z, y, y ′ ) and momentum scale q leading to singular limits at vanishing or infinite values of these parameters due to the appearance of contact terms beyond tree level.In particular, the limit z → 0, despite the fact that classically it results in a straight Wilson line of length (y − y ′ ), gives rise to a linear ∼ 1/|z| and a logarithmic ∼ ln(z 2 ) divergence in the conversion matrices of O Γ coming from the cusp and pinch-pole divergent terms of the original operator: ). (III.3.34) Setting y ′ = y, and considering the limit y → 0, where the staple becomes a straight line of length z, cusp points do not vanish smoothly giving logarithmic divergences ∼ ln(y 2 ): As discussed in Sec.III.1, the limit y → ∞, when (y − y ′ ) is fixed, results in a pinch-pole linear singularity (∼ |y|): . In contrast to the aforementioned limits, the limit (y − y ′ ) → 0, in which the asymmetric staple Wilson line becomes symmetric, is not singular, and thus, our present results can reproduce our previous results in Ref. [108] for the case of symmetric staples.Furthermore, the limit q2 → 0 is nonsmooth for some of the operators O Γ giving a logarithmic divergence (∼ ln(q 2 )), as follows: ), (III.3.37) where d Γ = −3 for (Γ = 1 1, γ 5 ), 1 for (Γ = σ νiνj ), and 0 for (Γ = γ νi , γ 5 γ νi ).Note that the limit is taken by rescaling simultaneously all components of the 4-vector scale q.Moreover, the limit q2 → ∞ also gives a logarithmic divergence, which is independent of the operator O Γ : Another property of the conversion matrices comes from the combination of P µ , T µ , C symmetries: one can prove (see Eq. (II.3.5)) that the real (imaginary) parts of the conversion matrices are even (odd) under (z → −z, y → −y, y ′ → −y ′ ).This is confirmed by our one-loop computation.
We illustrate our results for the conversion matrices in the plots of Figs.
(5 -7), we examine the dependence of some representative conversion matrix elements on the Wilson-line scales y, z, and y − y ′ (in lattice units), respectively.In particular, we plot the real and imaginary parts of the 2 × 2 block C MS,RI ′ ΓΓ ′ , (Γ, Γ ′ = γ ν1 , γ ν2 ), for the four renormalization prescriptions: RI ′ 1 , RI ′ 2 , RI ′ 1 -bar, RI ′ 2 -bar.As shown in Fig. 5, the real part of the diagonal elements (first plot in the left column and second plot in the right column) has an almost linear (flat) dependence on y/a in the RI ′ i (RI ′ i -bar) schemes due to the presence (absence) of the pinch-pole singularity.As expected, the convergent behavior of RI ′ i -bar appears for large values of y/a (≳ 6).There are no significant differences in the real diagonal elements between schemes with index 1 and 2. The real part of the nondiagonal elements (first plot in the second column and second plot in the first column of Fig. 5) also converges for large values of y/a (≳ 8) for all RI ′ -type schemes.Note that the nondiagonal elements are identical between RI ′ i and RI ′ i -bar schemes (for the same i) at one loop.Differences between schemes with index 1 and 2 are now visible: the real nondiagonal elements for the RI ′ 2 scheme have smaller absolute values compared to RI ′ 1 , which are closer to zero.The relative size of the real nondiagonal elements compared to the diagonal ones is ∼ 20%.The imaginary part of the conversion matrix elements (the four plots in the last two rows of Fig. 5) gives a much milder contribution compared to the real one (≲ 10%).The imaginary (diagonal and nondiagonal) elements for RI ′ i -bar coincide with the corresponding elements for RI ′ i (for the same i) at one loop.A plateau is observed at large values of y/a; a more stable behavior is seen for the nondiagonal elements.As in the case of real nondiagonal elements, noticeable distinctions can also be spotted in the imaginary parts between schemes with index 1 and 2. However, there is no consistent pattern regarding which scheme yields smaller contributions, as it varies for each element.
Similar conclusions are extracted from Fig. 6 by considering the dependence of the conversion matrix elements on z/a.Here, we observe a convergent dependence for z/a ≳ 6 for all conversion matrix elements, except for the imaginary diagonal parts (the four plots in the last two rows of Fig. 6).Now, the real diagonal parts have flat behavior in both RI ′ i and RI ′ i -bar schemes since the pinch-pole singularity does not arise for large values of z/a.Examining the dependence of the conversion matrix elements on (y − y ′ )/a in Fig. 7, we conclude that a more flat behavior is observed for smaller values of (y − y ′ )/a, while for larger values the elements decrease (in most cases) rapidly.An almost linear dependence on (y − y ′ )/a with a negative slope is obtained in the real diagonal parts of the RI ′ i schemes (see the first plot in the first column and the second plot in the second column of Fig. 7), coming from the pinch-pole divergent term.
In summary, the conversion matrices exhibit significant contributions from the real diagonal components, whereas the imaginary diagonal and nondiagonal elements make comparatively milder yet perceptible contributions.Safer conclusions by comparing the different types of renormalization schemes can be obtained when combining the conversion matrices with the nonperturbative data.

IV.1. Green's functions
In the lattice calculation, we employ the Wilson/clover fermion action (see [130]) and a family of gluon Symanzik improved actions [131] of the form: Re Tr {1 − U rect. } , (IV.1.1) where U plaq. and U rect. are the standard 4-link "plaquette", and 6-link "2×1 rectangle" Wilson loops.We selected three of the most common choices of the Symanzik coefficients c i , called Wilson, Tree-Level Symanzik, and Iwasaki gluon actions, as shown in Table II.Since we consider mass-independent renormalization schemes, we set the quark mass equal to zero.Consequently, the results from this study are also applicable to the twisted mass fermions [132] in the chiral limit.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.).
The results for the one-loop lattice bare Green's functions of the asymmetric staple operators Λ LR Γ are presented below in terms of the MS-renormalized Green's functions, derived by the corresponding calculation in DR.The methodology for calculating the one-loop momentum integrals on the lattice is described in Refs.[105,108].Here,  we cite results for a general Wilson-line lattice operator with n cusps [(n + 1) segments] and no self-intersections, as calculated in our previous study considering symmetric staples [108].We confirm that the formula constructed in the latter publication gives the correct result for the difference between the bare lattice Green's functions and the MS-renormalized Green's functions, also in the case of asymmetric staple operators studied in this work.Thus, for a general Wilson-line lattice operator we have: where r is the 4-vector that connects the two endpoints of the Wilson line, ℓ j is the length of the j th straight-line segment and l ≡ n+1 j=1 ℓ j is the total length of the Wilson line, c SW is the clover coefficient in the fermion action, and νi (ν f ) is the direction of the Wilson line in the initial (final) endpoint.P 2 = 0.02401318111946489(1) [133]    Conclusions from this calculation are summarized below: 1. Linear divergence: On the lattice, there is a linear divergence, which depends on the length of the Wilson line and the gluon action that is employed.This divergence comes from diagrams contributing to the Wilson-line self-energy.At one loop, the contributing diagram is only d 4 .
2. Logarithmic divergences: In both regularizations (DR and lattice), there are end-point logarithmic divergences, coming from diagrams d 2 and d 3 , and cusp and contact logarithmic divergences, coming from diagram d 4 .The coefficients in front of these divergences are regularization-independent.
3. Operator mixing: Diagrams d 2 and d 3 give rise (upon summation) to the Dirac structure Γ ′ = (Γ / ν i + / ν f Γ)/2, which differs from the tree-level structure Γ of the operator.This indicates that operator mixing is present: In order to remove this additional structure, we need to renormalize the Wilson-line operators O Γ and O Γ ′ as a doublet by introducing a 2 × 2 mixing renormalization matrix.However, as concluded by symmetries, the employment of 4 × 4 mixing matrices for renormalizing quadruplets of asymmetric staple-shaped operators is expected to be required at higher loops.The mixing contributions at one loop depend solely on the direction of the Wilson line entering the endpoints, regardless of the shape of the Wilson line.Also, the coefficient α 2 + α 3 c SW in front of the structure Γ ′ depends on r and c SW ; in particular, α 2 vanishes when r = 0. Thus, the one-loop mixing contributions originate from the chirality-breaking parts of the fermion action.As concluded by symmetries, these specific contributions are expected to be absent when a chiral-fermion action is employed.Hence, the exact shape of the Wilson line does not affect the one-loop continuum and lattice Green's functions of the Wilson-line operators differently; the only additional contributions on the lattice depend on the total length, the number of cusps, and the direction of the Wilson line entering the endpoints.
In order to investigate RI ′ -bar schemes on the lattice, we have also calculated the one-loop bare Green's function of the Wilson loop ⟨L LR (z, y + y ′ )⟩, given below in terms of the corresponding MS-renormalized Green's function: ⟨L MS (z, y + y ′ )⟩ can be read from Eq. (III.1.44)by removing 1/ε terms.The result agrees with [129] in the case of Wilson gluon action.As expected, the linearly divergent term (1/a) of Eq. (IV.1.3)cancels the linearly divergent term of Eq. (IV.1.2) when calculating Λ LR Γ (q, z, y, y ′ ) (Eq. (II.4.10)).

IV.2. Renormalization matrices
The lattice renormalization matrices of the Wilson-line operators O Γ in the MS scheme (Z MS,LR ΓΓ ′ ) can be extracted from Eq. (IV.1.2) by imposing that the terms in the curly bracket are canceled when renormalizing both operator and external fermion fields in the lattice Green's function.
The lattice renormalization matrices for the modified asymmetric staple-shaped Wilson-line operators O Γ in the MS scheme have the same form as (IV.2.2).Their elements are given by: The linear divergence is now absent from the renormalization matrices.The nondiagonal elements are identical to (IV.2.2).

V. SUMMARY AND FUTURE PLANS
In this work, we present an extensive and comprehensive study of the renormalization of nonlocal quark bilinear operators featuring an asymmetric staple-shaped Wilson line.This project is motivated by the increased interest in studying TMDPDFs from lattice QCD using novel approaches, such as large momentum effective theory and shortdistance factorization, which require matrix elements of the operators under study.More details can be found in the TMD Handbook [134].
The analysis is based on a one-loop perturbative calculation of Green's functions of such operators in both lattice and continuum (dimensional) regularizations.Based on our results, we identify the mixing pattern of these operators and propose renormalization prescriptions applicable to perturbative and nonperturbative data.More precisely, we discuss RI ′ -type conditions by using different projectors that effectively address power and logarithmic divergences, as well as the finite mixing among operators with different Dirac structures.We have systematically analyzed the mixing patterns of these operators, leveraging symmetry arguments for both chiral and nonchiral fermions.We also introduce a variant of the RI ′ prescription, which removes the pinch-pole singularities inherent in staple operators of infinite length by incorporating rectangular Wilson loops.This strategy also eliminates residual power divergences.Another novel aspect of this work is the extraction of the conversion factors to the MS scheme using the results in dimensional regularization.Our calculations are performed for arbitrary values of the renormalization momentum scale and the spatial dimensions of the staple.This ensures their applicability across a broad spectrum of nonperturbative investigations that may use the results of this work.

1 .
Translational symmetry: The operators are covariant under translations in Euclidean space.Similarly to the case of local operators O(x) which cannot mix with translated versions of themselves (O(y)), mixing among nonlocal operators involving different paths and shapes for the Wilson line joining the fermion-antifermion pair cannot occur[98,100,103].

d 1 d 2 d 3 d 4 FIG. 2 :
FIG.2: One-loop Feynman diagrams contributing to the Green's functions of the asymmetric staple-shaped operator with external fermions.The straight (wavy) lines represent fermions (gluons).The operator insertion is denoted by a filled rectangle.

6 FIG. 7 :
FIG. 7: Real and imaginary parts of the conversion matrix elements C MS,RI ′ ΓΓ ′
and α i are numerical constants which depend on the gluon action and the Wilson parameter of the fermion action r; their values are given in Table III for the Wilson, Tree-level Symanzik and Iwasaki gluon actions and for r = 1.The first two terms in the curly brackets of Eq. (IV.1.2) come from the sum of Feynman diagrams d 2 and d 3 , while the last term comes from diagram d 4 .The expression for the case under study (asymmetric staple operator with two cusps) can be extracted by setting n = 2, l = |z| + |y| + |y ′ |, νi = sgn(y) ν2 , and νf = −sgn(y ′ ) ν2 .The corresponding MS-renormalized Green's function Λ MS Γ can be read from Eqs. (III.1.1 -III.1.41)by removing 1/ε terms.

4 .
Finite contributions: Diagram d 1 is identical in DR and lattice regularization (up to discretization effects O(a)) giving no contribution to the difference δΛ Γ .Finite contributions of diagrams d 2 and d 3 stem only from the endpoints of the Wilson line.Any parts of a segment that do not include the endpoints give finite contributions which differ between DR and lattice regularization only by discretization effects.In diagram d 4 , the finite contributions come from the cusps and the straight-line segments of the Wilson line; they depend on both the number of cusps (n) and the number of segments (n + 1).

TABLE II :
Selected sets of values for the Symanzik coefficients.In all cases, they satisfy c 0 + 8c 1 = 1.

TABLE III :
Numerical values of the coefficients α 1 − α 6 appearing in the difference δΛ Γ of Eq. (IV.1.2) for r = 1.A systematic error is quoted coming from the numerical integration over loop momenta.