Nonperturbative renormalization of asymmetric staple-shaped operators in twisted mass lattice QCD

Staple-shaped Wilson line operators are necessary for the study of transverse momentum-dependent parton distribution functions (TMDPDFs) in lattice QCD and beyond. In this work, we study the renormalization of such operators in the general case of an asymmetric staple. We analyze the mixing pattern of these operators using their symmetry properties, where we find that the possible mixing is restricted within groups of four operators. We then present numerical results using the regularization independent momentum subtraction (RI/MOM) scheme to study the importance of mixing using one operator in particular, the $\gamma_0$ operator. Based on these results, we consider the short distance ratio (SDR) scheme, which is desirable in the absence of mixing. Finally, we investigate a variant of the RI/MOM scheme, where the renormalization factors are computed at short distances.


I. INTRODUCTION
Collinear parton distribution functions (PDFs) probe the hadron structure from the perspective of the spin and longitudinal momentum distributions of the quarks and gluons that make up the hadron.The vast amount of work to determine these distributions over the last five decades, both theoretically and experimentally, has enormously expanded our view of the hadron structure, the proton in particular [1].And yet, these developments are mostly limited to probing the one-dimensional structure of the proton.In order to have a wider understanding of the proton structure, we also need to understand how the momentum and the spin are distributed in the transverse plane.For that, we need to measure and compute the generalized parton distributions and the transverse-momentum-dependent PDFs (TMDPDFs), the latter being the main focus of this present work.Although there has been an effort to obtain TMDPDFs from phenomenological fits to experimental data [2][3][4][5][6][7], their accuracy is still far from being at the same level as the collinear PDFs.This status will change in the coming years with new data coming from Jefferson Lab [8] and from the future Electron-Ion-Collider to be built at Brookhaven National Lab [9].It is, thus, of great importance to extract TMDPDFs from first principles calculations, namely, lattice QCD.In the last 10 years, there has been an enormous advance in computing the collinear PDFs using lattice QCD [10][11][12][13][14][15][16][17][18].By comparison, the computation of TMDPDFs is still in its infancy [19][20][21][22], although progress has been made in the computation of key elements required to build the TMDPDFS, like the soft functions [23,24] and the Collins-Soper kernel [25][26][27].Such computations have been restricted to ratios of TMDPDFs [28][29][30][31], and only recently, a first full lattice QCD calculation of the TMDPDFs themselves has been presented [32].One fundamental difficulty in these calculations is to have control over the renormalization procedure, which is more involved than in the case of collinear PDFs.In particular, the staple-shaped gauge link that enters in the evaluation of TMDPDFs has three types of divergences: i) linear divergence coming from the Wilson line, which connects the quark fields and which depends on the length of the staple-shaped link; ii) logarithmic divergences coming form the endpoints of the staple link, similarly to the case of straight gauge links; and iii) logarithmic divergences coming from the presence of cusps in the staple.In addition, in the limit of infinite-length staple, L, (which is the case of interest), pinch-pole singularities arise as positive powers of L, coming from the gluon exchange between the transverse segments of the staple.Moreover, as in the case of straight gauge links, staple-shaped operators of different Dirac structures can and will mix on the lattice among certain groups (when chiral symmetry is broken), as dictated by discrete symmetries.However, in this case, the mixing pattern of the operators employed in lattice regularization can be significantly more involved than in the case of the straight Wilson line.A first study within lattice perturbation theory to one-loop for the case of the symmetric staple [33] showed mixing between specific pairs of Dirac structures.This mixing cannot be avoided when one is interested in matching bare lattice Green's functions to the MS scheme (directly or indirectly through an intermediate scheme).The mixing depended solely on the direction of the staple link entering the endpoints of the staple, regardless of the shape of the staple.This implies that the same mixing pattern occurs also for asymmetric staple shapes.In Ref. [34], the mixing pattern of these staple-shaped operators1 has been studied using symmetry considerations.It was found that more mixing is present than observed in Ref. [33].This demonstrates that one-loop perturbation theory cannot fully reveal the mixing of the staple-shaped operators, unlike the case of straight Wilson-line operators [35], and a higher-loop computation is needed to confirm the additional mixing patterns found by symmetry arguments.The authors of Refs.[23,27,36] consider a maximal RI-type prescription, in which all 16 independent non-local Wilsonline quark bilinear operators are chosen to mix to eliminate all possible mixing effects.While Ref. [27] identified nonzero contributions in several off-diagonal elements of the renormalization matrix, Refs.[23,36] found negligible contributions, at least at small transverse separations, by setting specific momentum components to be zero.However, we emphasize that not all contributions are necessary for addressing the "unavoidable" mixing among the asymmetric staple-shaped operators on the lattice.Most off-diagonal elements in the 16 × 16 renormalization matrix are nonzero due to the non-minimal choice of renormalization conditions and not due to the unavoidable mixing.In this sense, it is preferable to construct a minimal intermediate scheme, keeping only the mixing sets that are needed for matching the bare lattice Green's functions to the corresponding MS-renormalized Green's functions (as obtained in dimensional regularization).In our study, we consider such a minimal scheme by using symmetry arguments to restrict the operators allowed to mix.Improving renormalization schemes on the lattice eliminates finite lattice spacing errors, which can come from different Dirac structures in Green's functions under consideration.A way of removing artifacts from all Dirac structures is to consider a wider mixing pattern, where higher dimensional operators multiplied by the appropriate power of the lattice spacing can also mix with the operators under study.This mixing is only present for finite values of the lattice spacing, while it vanishes when taking the continuum limit.The higher-dimensional operators will be higher twist since (by Lorentz invariance) they must have the same spin (with the operators under study), but their dimension will be higher.The unwanted effects of finite lattice spacing errors and higher-twist contributions are not considered in the present study; they will be addressed in a future publication.On the other hand, it has been shown [37] that the linear divergence in the lattice spacing a is not fully eliminated when the RI/MOM scheme is used to renormalize a straight Wilson line of length z.Ref. [36] has shown that this residual linear divergence remains in the case of the staple-shaped operator.In this scenario, an alternative approach that one can adopt is the so-called ratio scheme as proposed for the quasi-PDF case [11,38,39].In this approach, one subtracts the Ultra-Violet (UV) divergences by taking the ratio with a suitable object at a fixed short distance, where perturbation theory applies.Different choices of suitable objects to be used in the ratio have been proposed in Ref. [40].The authors of Refs.[32,36] use ratios of the matrix elements of the operator under study.As the divergences of the staple-shaped operator are independent of the longitudinal momentum of the state, one can use the matrix elements computed at different values of the momentum in order to cancel the divergences.Usually, matrix elements at zero momentum are chosen for the denominator, and such scheme has been named [32] short distance ratio (SDR) scheme.
As for the remaining divergences associated with the asymmetric staple-shaped link, one can cancel them by taking an appropriate ratio with the vacuum expectation value of a rectangular Wilson loop [41].However, we would like to stress that the SDR scheme is valid when operator mixing is absent or negligible.We check that in the case under study the mixing is indeed negligible and thus, one can employ the ratio scheme.Another option, as proposed in Ref. [34], is to employ RI/MOM in the spirit of the SDR scheme by fixing the dimensions of the staple at a short perturbative range (we will call this scheme RI-short).The paper is organized as follows.Section II discusses how TMDs can be accessed on a Euclidean lattice using the quasi-distribution approach.In Section III, we study operator mixing using symmetry arguments.Section IV presents our lattice setup.The size of the mixing is examined in Section V.In Sections V A, V B and V C we employ RI/MOM, SDR and RI-short schemes, respectively, by following the procedures described in Ref. [32,34,36].In Section VI, we present the renormalized beam functions in the three renormalization schemes discussed here, and then show the corresponding results in the MS scheme.The conclusions are presented in Section VII.

II. QUASI-TMDPDFS
A first-principle computation of TMDPDFs in the context of LaMET is more involved than collinear PDFs.The main obstacle is to correctly subtract the so-called rapidity divergences, which are associated with the soft gluon radiation of the infinitely long Wilson lines present in the staple-shaped operator.This subtraction is made through the use of a soft function [42,43].Different ways how to perform this subtraction in the context of LaMET can be found in Ref. [44].In general, the TMD soft function involves two opposite light-like directions, and this makes a lattice calculation significantly challenging.In Refs.[21,41], the authors define a rapidity-independent reduced soft function S r and show that it can be extracted from a form factor and the quasi-TMD wave function of a light meson.Using S r , the rapidity scheme-independent TMDPDF f T M D (x, b, µ, ζ) can be written as [15,41] where f (x, b, µ, ζ z ) is the so-called quasi-TMDPDF for a nucleon with mass M , x is the longitudinal momentum fraction, and b is a separation transverse to direction of the momentum P z carried by the nucleon.The scale µ defines the renormalization scale and ζ z = (2xP z ) 2 is the Collins-Soper scale of the quasi-TMDPDF, with ζ being the scale for the light-cone correlation.The factor H ζz µ 2 is the perturbative matching kernel that connects the TMDPDFs to the quasi-TMDPDFs, and K(b, µ) is the Collins-Soper kernel.The renormalization scheme in Eq. ( 1) is left unspecified; usually, it is computed in the MS scheme.The quasi-TMDPDF on the lattice is defined as [20,44] f with B Γ (z, b, µ, P z ) the renormalized beam function.It is obtained from the bare beam function, B 0,Γ (z, b, L, P z ; 1/a), defined as the matrix element of the non-local staple-shaped quark bilinear operator, O Γ (z, b, L), of length L: where N (P z ) is a nucleon state with momentum boost of (0, 0, P z ), a is the lattice spacing, W(z, b, L; 1/a) is the staple-shaped Wilson line, and ψ(z) is the standard up and down quark doublet.In practice, we compute the flavor non-singlet combination u − d, and hence the operator also includes a Pauli τ 3 matrix in flavor space.For the unpolarized TMDPDF, Γ can be either γ 0 or γ 3 or, in general, a combination of the two.In this work, we show results for Γ = γ 0 , since we observe a better signal for this operator compared to Γ = γ 3 .γ 0 was also used for the quasi-PDF case, in which case mixing is absent, while γ 3 mixed with 1.In the case of the staple, however, no such advantage exists, as will be shown in Section III.In general, W(z, b, L; 1/a) is given by: where we dropped the dependence on 1/a to keep the notation compact.The arguments of the Wilson lines on the right hand side (RHS) of the equation above are defined in the following way: the indexes z and ⊥ in W z and W ⊥ refers to the direction ẑ of the boost and a direction perpendicular to it, respectively.The transverse direction could be in the x or ŷ direction.The first argument of both W z and W ⊥ refers to the initial point of the Wilson line, while the second argument is the displacement of the Wilson line in the direction of the indexes (z or ⊥).Mathematically, the expression for a Wilson line starting at the point ⃗ x, with a displacement L in the ẑ direction, is given by: with P the path order operator and A the gluon field.The shape of the asymmetric staple is shown in Fig. 1.
The shape of the asymmetric staple defined in the operator of the quasi-beam function

III. OPERATOR MIXING THROUGH SYMMETRY
In order to use symmetry arguments to determine which operators mix, we first organize the staple-shaped operator using all possible ways to connect the fermion fields.We show in Fig. 2 all the distinct possibilities of connecting a quark and an antiquark field located at a spatial separation of (0, b, z), using an asymmetric staple-shaped Wilson line.
Since we build the staples only in the yz-plane, in the following discussion and in Fig. 2, we denote the coordinate space (t, x, y, z) as only (y, z).Without loss of generality, we can fix the quark field (q) at (0, 0).Then, the only possibilities for the position of the antiquark field (q) are (b, z), (b, −z), (−b, z) and (−b, −z).Here, b and z are strictly non-negative.We define the staple-shaped Wilson line connecting q at (0, 0) and q at (b, z) as O ++ (Γ).The two plus signs denote the direction of separation of q from q.The first sign is for the transverse (ŷ) and the second sign for the longitudinal (ẑ) axes.The Γ denotes the insertion gamma matrix as defined in Eq. 3. In a similar fashion, the staple-shaped operator in the cases of q being located at (b, −z), We analyze the symmetry properties using generalized parity (P 1,2 F α ) and time reversal (T 1,2 F α ) transformations with discrete flavor rotation, as well as charge conjugation (C), for the fermion fields in the twisted mass basis, χ(x α , x), and with the gauge link, U (x α , x; α), in some direction α.The symmetry transformations are defined in Appendix A. The operators O ±± (Γ), O ±± c (Γ) do not have definite properties with respect to the symmetries of the lattice action.Instead, one needs to consider their linear combinations, defined by The only combinations that have definite symmetry properties are: FIG. 2: All distinct possibilities for an asymmetric staple-shaped Wilson line connecting a quark and an antiquark field spatially separated by ±bŷ and ±z ẑ (b, z ≥ 0).
i.e. the combinations with all signs i, j, k, l reversed are equivalent from the point of view of symmetry transformations, with irrelevant global phase.
As an example, we look at the symmetry properties of γ 0 and γ 0 γ 3 given in Tables I and II, respectively, restricting ourselves here to the flavor non-singlet case, u − d (τ 3 matrix in flavor space).The combinations of operators that mix are those which have all signs equal in the nine rows of Tables I and II.For example, the symmetry properties of (+ + ++) c for Γ = γ 0 (second column of Table I) are identical to the ones of the combination (+ − −+) c for Γ = γ 0 γ 3 (last column of Table II).Thus, we conclude that the following mixings occur: where the first combination in a pair pertains to Γ = γ 0 and the second one to Γ = γ 0 γ 3 .Additional mixings appear with Γ = γ 0 γ 2 and Γ = γ 5 γ 1 , thus forming a quadruple of operators that mix, {γ 0 , γ 0 γ 2 , γ 0 γ 3 , γ 5 γ 1 }.These are the only relevant mixings for the present work, which concerns Γ = γ 0 for unpolarized TMDs.
For the general case, the symmetry properties for all Γ's are summarized in Table IV in Appendix A. They imply that the possible mixing is between Γ and {Γγ 2 , Γγ 3 , Γγ 2 γ 3 } in the asymmetric staples case.In turn, symmetric staples restrict the mixing by eliminating always one member of the mixing quadruple, i.e.Γ = γ 0,1,5 mix with Γγ 2 and Γγ 3 , while Γ = γ 2,3 form a triple {γ 2 , γ 3 , γ 2 γ 3 } of operators that mix.We point out that different mixing patterns have been considered in recent studies of other groups, including mixing among all 16 operators of different Dirac structures [27], or at least mixing in pairs (Γ, Γγ 3 ) [34].In our work, we choose to consider the minimal set of staple-shaped operators (of the same dimension) that are allowed to mix by the above-mentioned C, P, and T symmetries.We have not considered mixing with higher-dimensional operators allowed by Lorentz symmetry (see, e.g., Ref. [45] for the straight Wilson-line case), since it is power suppressed and not relevant when one takes the continuum limit a → 0. Also, in contrast to Ref. [34], in our analysis, the staple line has been chosen to be in a 2-dimensional (D), and not 3D, plane in Euclidean space formed by the transverse (ŷ) and longitudinal (ẑ) directions.In this respect, we end up with a basis of 8 (instead of 16) operators, which are eigenstates of C, P and T transformations.We also note that calculations in one-loop lattice perturbation theory [33] show a smaller mixing pattern compared to the symmetries; however, this cannot guarantee a reduced mixing in higher loops.

IV. LATTICE SETUP
For the lattice simulation, we use an N f = 2 + 1 + 1 clover-improved twisted mass fermion ensemble of size 24 3 × 48, produced by the Extended Twisted Mass Collaboration (ETMC) [46].Our study is done using N conf configurations with N src source positions for each configuration.To increase statistics, the boost is taken in all three directions, and both positive and negative.For each such direction of boost, the staple is then constructed in both the remaining transverse directions.This gives us 12 measurements (6 boost directions × 2 transverse directions) for each source position.The details of the lattice simulation are summarized in Table III.

TABLE III:
We give the parameters of the lattice ensemble and measurements used in the calculation.a is the lattice spacing, µ l is the bare twisted light quark mass, m π the pion mass, N conf the number of configurations, N src the number of source positions and N meas the total number of measurements.
The bare matrix element for the quasi-beam function is calculated through a ratio of a 3-point to a 2-point function, where τ is the insertion time of the operator O Γ and t s defines the source-sink time separation.P is the parity projector.In this work, we show results for a single source-sink separation of t s = 10a and a longitudinal momentum of P z = 6π/24a ∼ 1.7 GeV.The 3-point function is constructed for the isovector combination u − d, by inserting a τ 3 in flavor space.This choice ensures the elimination of the disconnected contributions and only connected diagrams need to be calculated.
Momentum smearing [47] is applied to improve the signal for large boosts.It is observed that applying stout smearing to the gauge links used in construction of the staple also reduces the statistical errors.Here, we have applied 5 steps of stout smearing to the staple-shaped Wilson line.
In Fig. 3 we show the bare matrix elements as a function of z at a fixed b/a = 1 for different values of L. These and the following results for the quasi-beam function have been symmetrized using the relation To ensure that this property holds in our lattice simulation, we have used the staple-shaped operator O ++ on the left hand side of the above equation and O −− c on the right hand side.Our results show a clear decay of the magnitude of the beam function as L increases, as expected [15].In Fig. 4 we show the dependence of the bare beam function on b.This is expected to decay as exp (−LV (b)) [15].Assuming a linear function for the potential V (b), we also show a fit of exp (c 0 + c 1 b) to the lattice data at two different values of z/a and a fixed L/a = 10.

V. NON-PERTURBATIVE RENORMALIZATION
As in the case of the straight Wilson line, there are logarithmic and linear divergences associated with the length z of the staple-shaped link.However, here there are two extra divergences: one, which is associated with the cusps of the staple, and one associated with the gluon exchange between the gauge links of length L → ∞, the so-called pinch-pole singularities.On the other hand, the finite transverse separation b mitigates possible discretization effects associated with small z separation, mainly in the non-perturbative region (large b), which is the main interest of a lattice calculation.In the following, we will analyze three different ways to carry out the renormalization, namely the RI/MOM scheme, the short distance ratio scheme, and a modified version of the RI/MOM scheme, where the renormalization constants are computed at short distances.

A. RI/MOM
The RI/MOM scheme [48] was first adapted for non-local operators employed in the quasi-distribution approach in Refs.[35,49].The RI/MOM renormalization constants Z RI O are defined by the condition where Λ Γ 0 is defined in terms of the amputated Green's function with S q the off-shell quark propagator.The Green's function is calculated as Because G Γ , and thus Λ Γ 0 , have the same divergences as B 0,Γ , all the divergences, in principle, cancel in the renormalization procedure.Z RI q is the quark wave function renormalization defined as Tr (S q (p; 1/a)) −1 S Born q (p) and the corresponding renormalized beam function is then Notice that the relation between the bare and renormalized beam functions involves, in principle, a 16 × 16 matrix, if one considers all possible mixing among the full set of Dirac structures.The case of interest in this work is Γ = γ 0 , computed using the lattice ensemble of Table III.As pointed out in Section III, γ 0 mixes with γ 0 γ 2 , γ 0 γ 3 , and γ 5 γ 1 .
The renormalization factors for the diagonal and off-diagonal cases are shown in Figs. 5 and 6, respectively.We set the renormalization scale 4-vector µ 0 equal to 2π 6+0.5 Nt , 3 Ns , 3 Ns , 3  Ns , where N s = 24a, N t = 2N s , and anti-periodic boundary conditions have been employed in the time direction.We have chosen an isotropic momentum in the spatial directions and "democratic" momentum, which obeys the criterion ρ sin 4 (ap ρ /2)/( ρ sin 2 (ap ρ /2)) 2 | p=µ0 < 0.3, in order to reduce Lorentz-non-invariant contributions in the vertex functions.Also, the selected momentum lies in a perturbative region, where perturbation theory is reliable and, at the same time, lattice artifacts are under control.A more detailed study using different values of µ 0 will be presented in a future extension of our study.Using this set-up, we observe that the contributions from the off-diagonal mixing terms are ≲ 4% of the diagonal contribution for all z when b = 1a, 2a.It is interesting to examine how the relative size of the off-diagonal to the diagonal renormalization factors changes as b increases, as the primary focus of lattice computation of TMDPDFs is the the non-perturbative region in the transverse separation b.We show in Fig. 7 the renormalization factors as a function of b for the diagonal and offdiagonal contributions.We see a contribution of ≲ 7 − 8% up to b = 6a for the off-diagonal terms.Although their size shows a possible tendency of growth with increasing b, this growth seems to be only moderate.We notice that the authors of Ref. [27] considered the entire set of Dirac structures for the operator mixing, and observe similar results to ours, if we take from Ref. [27] only the results from the operators that are allowed to mix with γ 0 , according to our symmetry arguments.For example, their result for the contribution from γ 5 γ 1 is much larger than that of γ 0 γ 2 and γ 0 γ 3 , and there is a steady increase going to larger b values.For the values of b we use in this work, the magnitude of mixing for these operators found in Ref. [27] is also comparable to our findings.However, considering the dominant diagonal renormalization factors, it seems that there are remaining divergences related to large valus of z and b, as also observed in Ref. [36].Combined with the fact that the bare matrix element decays exponentially with b [15], it is significantly hard to renormalize for large transverse separations using this scheme.However, in the small b region where we can extract information using RI/MOM, it can be confirmed that the mixing is negligible.As an example, we show in Fig. 8 the renormalized beam function for the b = 1a case, with +'s denoting the RI/MOM procedure using the full mixing matrix and with the ×'s denoting the case of only taking into account the diagonal contribution.
Considering that including the mixing has a negligible effect, one can safely ignore operator mixing for these values of b.

B. The short distance ratio scheme
According to Figs. 6 and 7, mixing different operators can be neglected at least up to b/a ≲ 6, and we can assume that the renormalization of the staple-shaped link is multiplicative.This justifies the approach taken in Refs.[32,36].We first note that the vacuum expectation value of a rectangular Wilson loop Z E with sides 2L + z and b, is, by construction, a product of the staple-shaped gauge link, as defined in Eq. ( 4), and its reflection.Therefore, √ Z E has the same divergences as that of the staple-shaped gauge link.Thus, it should cancel the pinch-pole singularity associated with the length L of the staple, as well as the divergences associated with the cusps.Furthermore, as the sides of Z E are also dependent on the longitudinal displacement z, the exponential divergence associated with z present in the staple-shaped link must also be cancelled if an appropriated ratio is taken, namely the one proposed in Refs.[21,36]: To illustrate this point, we show in Fig. 9 the beam function as a function of L, for a fixed z = 2a and two values of b, before and after dividing by √ Z E .The cancellation of the divergences related to the length of the Wilson line is explicitly shown.This ratio, thus, takes care of the divergence associated with the length L and width b of the staple-shaped operator.
After dividing the beam function by the square root of Z E , the only remaining divergences are the UV divergences associated with the quark field and its end-points connecting to the gauge links.As such, they can be cancelled by taking ratios, as they have a multiplicative nature [26,50].Hence, one can define where, Because the remaining divergences are independent of the length of the Wilson line, one is free to choose z 0 and b 0 .
In order to connect these quantities to the MS scheme via a perturbative scheme conversion, z 0 and b 0 should be small enough for perturbation theory to be valid.However, the use of small values for both z 0 and b 0 can introduce sizable discretization errors in the renormalization factors, which can affect the validity of the SDR scheme.To address this issue, different approaches can be employed in order to reduce finite lattice-spacing errors from the nonperturbative data for all values of the staple lengths b/a and z/a.Ideally, the elimination of discretization errors requires calculations of physical matrix elements at different finite values of the lattice spacing a and an extrapolation a → 0. When data for multiple values of a are not available, a number of different approaches can be employed in order to reduce discretization errors at each lattice spacing.A standard method is to apply an improved discretization, in both the action and the operators under study, using the Symanzik-improvement program [51,52].Another way to reduce this kind of systematic error from a lattice calculation is to subtract one-loop artifacts employing lattice perturbation theory to all orders in a, from the non-perturbative vertex functions calculated in lattice simulations.
Our group has successfully applied this method to the renormalization of local quark bilinear operators [53][54][55], and more recently to the renormalization of non-local straight Wilson-line operators for quasi-PDFs [56].These studies have provided a useful feedback on the effectiveness of artifacts in the renormalization factors for different ranges of the scales entering the renormalization procedure.Since this is our first non-perturbative study considering non-local staple-shaped operators we do not consider finite-a errors, but we intend to apply the method of subtracting one-loop artifacts in a future extension of our study.The renormalized matrix elements are converted to the MS scheme using perturbation theory.We have computed the vertex, sail, and tadpole one-loop diagrams for external quark states with a general momentum p z .For p z → 0, we obtain which agrees with Eq. ( 6) of Ref. [36].Note that this factor equals B MS Γ (z = z 0 , b = b 0 , P z = 0, µ 0 ).Details of this calculation for a general external momentum will be presented in Ref. [57].The renormalized beam function in the MS scheme is then given by Results for B MS (z, b, P z ) are presented in Section VI, where we also discuss the cancellation of the z 0 , b 0 dependence in Eq. ( 19).

C. Short distance RI/MOM
As discussed in Section V A, the usual RI scheme may be problematic at large z and b as the magnitude of the Z RI factors grows exponentially.Also, as shown in Refs.[36,37], the usual RI scheme may still contain a residual linear divergence, which may not be properly canceled.On the other hand, as discussed in Section V B, the √ Z E factor cancels all divergences in z and b present in the staple.Hence, we can define a vertex function that is free of such divergences, Because the divergences related to the lengths of the Wilson line have been removed, we can compute the renormalization factors as in Section V B at some fixed z 0 , b 0 [34] The Z RI−short factor defined at fixed z 0 and b 0 is used to renormalize the bare B γ0 defined at arbitrary values of z and b.We have labeled this procedure as RI-short.In principle, the vertex function in the standard RI/MOM could also be modified by taking its ratio with √ Z E .This would reduce the growth of the Z-factors with increasing b.However, when combining Z RI with the bare B γ0 , the √ Z E factors cancel each other, and therefore √ Z E has no effect on the renormalized matrix elements.On the other hand, the √ Z E factor appearing in the vertex function of the RI-short scheme is defined at fixed values of z and b, contrary to the √ Z E factor appearing in the bare B γ0 .Hence the cancellation of the Z E factors do not happen in the RI-short scheme.As in Section V B, we choose the pair z 0 , b 0 to be in the perturbative region in order to make the perturbative conversion to the MS scheme reliable.Moreover, the study of possible lattice artifacts associated with the use of small values of z 0 /a and b 0 /a will be considered in future extensions of our study by using one-loop lattice perturbation theory.The corresponding renormalized beam function in this scheme is then given by For illustration, we show in 48 , 3  24 , 3 24 , 3 24 .We also show the off-diagonal renormalization factors, which, in principle, mix with γ 0 .As expected, they are also independent of L and can be omitted in the full calculation as their contribution is negligible as compared to the diagonal contribution.Finally, we convert B RI−short Γ to the MS scheme using one-loop perturbation theory, The conversion matrix Z MS,RI−short O is calculated in dimensional regularization for arbitrary values of the momentum scale µ 0 .Explicit expressions for all Γ, Γ ′ are given in Ref. [57].Similar perturbative studies can be found in Refs.[26,33].A comparison between the conversion factor in RI-short and the SDR scheme is shown in Fig. 11 for the case Γ = Γ ′ = γ 0 .We observe that the real parts are compatible between the two schemes, while the RI-short scheme shows a non-zero imaginary part, in contrast to the SDR scheme.This is a consequence of using a non-zero momentum scale µ 0 in the RI-short scheme, while SDR is defined at zero momentum.Note that specific choices of µ 0 in RI-short can lead to a vanishing imaginary part in the renormalization factors.In particular, by setting to zero the two momentum components parallel to the staple segments, the one-loop expression for the renormalization factors results into a vanishing imaginary part.However, such choice of momentum gives rise to unwanted Lorentz non-invariant contributions in the non-perturbative calculations.Thus, in our study we follow the common practice of employing democratic momentum with reduced Lorentz non-invariant contributions at the cost of introducing imaginary part in the renormalization factors.

VI. RESULTS
In this Section, we present the renormalized beam functions in the 3 renormalization schemes discussed in Section V, as well as the corresponding results in the MS scheme.In Fig. 12, we show the real and imaginary parts of B γ0 (z, b, µ 0 , P z ) as a function of z for two values of the transverse separation, b/a = 1, 2. For the real parts, the schemes agree within errors for |z/a| ⪅ 6.For |z/a| ⪆ 6, B SDR γ0 and B RI−short γ0 are consistent, while B RI γ0 becomes increasingly more negative with increasingly larger errors.This discrepant behaviour between the RI-MOM and the RI-short and SDR schemes indicates the existence of a residual linear divergence at large |z| in the standard RI/MOM, as observed in [36].For the imaginary parts, we observe a similar behavior.Namely, the SDR and the RI-short schemes are almost identical for any |z/a| for the two values of b/a considered.For the usual RI/MOM, we observe deviations not only for |z/a| ⪆ 6 but also for most of the |z/a| region with an increase in errors with increasing |z/a|.The conversion to the MS scheme is done using one-loop perturbation theory for all three schemes considered.For the two cases where the renormalization factors are fixed at a short perturbative scale, we use z 0 /a = b 0 /a = 1 when X=RI-short X=SDR (z, z 0 , b, b 0 , µ 0 , P z ) and B SDR Γ (z, z 0 , b, b 0 , P z ).Since the final results in the MS scheme should be independent of the values we choose for the z 0 , b 0 pair, we examine the stability of our results on the choice of the values for this pair.To this end, we show in Fig. 13 the product of the renormalization factors for each scheme with the corresponding conversion factors to the MS scheme, Z MS,SDR .These products should be independent on the choice of values of z 0 , b 0 as long as the one-loop perturbative calculation is a reliable approximation.We notice that these products are equal to the renormalization factors in the MS scheme in the absence of mixing.For the ratio scheme, there is a 10% − 20% correction when going from b 0 /a = 1 to b 0 /a = 2 for the three values of z 0 considered here, with the correction increasing for larger values of z 0 .This indicates a limitation on the use of one-loop perturbation theory to perform the conversion to the MS scheme already at a transverse separation as low as b/a = 2.As a consequence, O(α 2 s ) corrections cannot be disregarded.For the z 0 dependence, the choice of z 0 /a = 0 and z 0 /a = 1 is nearly equivalent within one sigma error.For z/a = 2, the discrepancy grows, especially for b/a = 2.For the RI-short scheme, the situation is not significantly changed.From this discussion, we conclude that our choice for the pair z 0 , b 0 is a likely safe region to fix the short perturbative scale.We show in Fig. 14 the renormalized MS matrix elements computed in the ratio and the RI-short schemes at P z ≃ 1.7 GeV.The imaginary part of the matrix elements are in full agreement, the same happening in the the real part, for large z values.In the small z region of the real part there is a tendency to discrepancy, as already observed in the intermediate schemes shown in Fig. 12, although they agree within one sigma error.

VII. CONCLUSION
In this work, we have addressed the non-perturbative renormalization of the asymmetric staple-shaped quark bilinear operators.We used symmetry arguments to restrict the possible mixing between the operators, and we found that mixing is allowed within sets of four operators.We then employed the RI/MOM scheme to estimate the importance of mixing coming from the non-diagonal terms, and we found that they can be neglected up to transverse separations of at least 6a.This result justifies the use of multiplicative renormalization for the asymmetric staple-shaped operator.Based on this conclusion, we computed the renormalization factors in the SDR scheme, and in the RI-short scheme, where the renormalization factors are computed at short distances.We found that both schemes are nearly equivalent over all the |z| region considered.Finally, we converted both schemes to the MS scheme and presented the corresponding results for the renormalized beam function.Given our previous computation of the Collins-Soper kernel and of the soft function [24], the next natural step will be to compute the TMDPDFs themselves over a wide range of b, and we plan to present results for them in the near future.

(A6)
The mixing properties of the staple-shaped operators is determined from symmetry arguments.Specifically, we observe the transformations of the objects (ijkl) c defined by Eq. ( 6).This gives rise to 4 relevant combinations with definite symmetry properties: (+ + ++) c , (+ − +−) c , (+ + −−) c , and (+ − −+) c .In Table IV, we summarize the symmetries of the 16 Dirac structures, where the four signs under any transformation denote the symmetry properties of the four relevant combinations, see the caption for more details.).The last row indicates the symmetry properties with respect to charge conjugation that depend on the sign c and are common to all combinations.The combinations that mix have the same signs in a given column (see e.g. the first column of γ 0 , the second column of γ 0 γ 2 , the third column of γ 5 γ 1 and the fourth column of γ 0 γ 3 , representing the quadruple {γ 0 , γ 0 γ 2 , γ 0 γ 3 , γ 5 γ 1 } of operators that mix relevant to this work; note that the mixing with γ 5 γ 1 involves the combination (+ + −−) −c with opposite relative sign in front of the charge-conjugated operators).
In the case of symmetric staple-shaped operators (z = 0), the mixing sets are reduced to: • {γ 2 , γ 3 , γ 2 γ 3 }, −b, z) and (−b, −z) can be defined as O +− (Γ), O −+ (Γ) and O −− (Γ) respectively.A visual representation of these operators is shown in Fig. 2 in black.Due to the asymmetric nature of the staple-shaped Wilson line, there are four more operators that are obtained from the charge conjugation of the above defined four.These are also shown in Fig. 2, but in red, with the charge-conjugated version of O ±± (Γ) denoted by O ±± c (Γ) and quark/antiquark fields having exchanged positions.For the symmetric case (z = 0), the charge-conjugated operators are redundant, since O ++ c with i, j, k, l = ±1 denoting the signs of O ++ , O −− , O +− , O −+ , and c = ±1 representing the relative sign of the charge-conjugated versions of O ±± .

10 FIG. 3 :
FIG.3: Real and imaginary parts of the bare matrix elements for a transverse separation of b/a = 1.

2 z a = 4 FIG. 4 :
FIG.4: Exponential decay of the bare matrix element with increasing b at L = 10a and z/a = 2 (red points) and z/a = 4 (blue points).The red and blue bands are the fits to the results at z/a = 2 and z/a = 4, respectively.

FIG. 9 :
FIG. 9: The effect of taking the ratio of the bare matrix elements with√ Z E at a fixed z/a = 2.

= 1 FIG. 10 :
FIG.10:The Z RI−short renormalization factor within the perturbative range b 0 /a = 1 as a function of the length L of the staple.

= 2 FIG. 11 :
FIG. 11: Conversion factor to the MS scheme (top: real part, bottom: imaginary part) for the RI-short and SDR schemes at b/a = 1 (left) and b/a = 2 (right).

FIG. 13 : 2 FIG. 14 :U
FIG.13: Dependence of the product between the MS conversion factors and the SDR and RI − short renormalization factors on the choice of the short perturbative scale.

TABLE I :
Symmetry properties of operators with Γ = γ 0 .The ± sign for P F /T F transformations denotes that a given combination of staple-shaped operators is symmetric/antisymmetric with respect to the symmetry transformation given in the first column.The last row indicates symmetry properties with respect to charge conjugation, which depend on the sign c, i.e. (• • ••) + combinations are symmetric and (• • ••) − antisymmetric with respect to C in this case.

TABLE II :
Symmetry properties of operators with Γ = γ 0 γ 3 .See the caption of Fig.I for explanation.