Nonperturbative renormalization of staple-shaped Wilson line operators in lattice QCD

Quark bilinear operators with staple-shaped Wilson lines are used to study transverse-momentumdependent parton distribution functions (TMDPDFs) from lattice quantum chromodynamics (QCD). Here, the renormalization factors for the isovector operators, including all mixings between operators with different Dirac structures, are computed nonperturbatively in the regularizationindependent momentum subtraction scheme for the first time. This study is undertaken in quenched QCD with three different lattice spacings. With Wilson flow applied to the gauge fields in the calculations, the operator mixing pattern due to chiral symmetry breaking with the lattice regularization is found to be significantly different from that predicted by one-loop lattice perturbation theory calculations. These results constitute a critical step towards the systematic extraction of TMDPDFs from lattice QCD.


I. INTRODUCTION
Building a quantitative description of the structure of the proton in terms of its fundamental parton constituents is a defining goal of hadronic physics research. A key aspect of this structure is encoded in transverse-momentum-dependent parton distribution functions (TMDPDFs), which describe the intrinsic transverse momentum of partons in the proton [1][2][3]. When the transverse momentum of the parton, q T , is in the perturbative region of quantum chromodynamics (QCD), i.e. q T Λ QCD , the TMDPDFs can be obtained perturbatively in terms of collinear parton distribution functions (PDFs) [3,4]. In contrast, when q T ∼ Λ QCD , the TMDPDFs are intrinsically nonperturbative, and constraining these fundamental aspects of proton structure remains a challenging problem for both theory and experiment.
Lattice QCD studies of TMDPDFs involve the nonperturbative calculation of hadron matrix elements of nonlocal bilinear operators with staple-shaped Wilson lines. These matrix elements are referred to as unsubtracted quasi TMDPDFs [26,27] or quasi beam functions [28,29]. An important component of such calculations is the renormalization of the bare quasi beam functions and their matching to the MS scheme. Bare quasi beam functions display both logarithmic and linear ultraviolet (UV) divergences. The linear divergences originate from the self-energies of the Wilson lines, and can be absorbed into exponential factors [32][33][34][35][36][37][38][39]. For hadronic matrix elements of quark bilinear operators with straight Wilson lines, which define the quasi PDFs, renormalization has been extensively studied in both perturbative and nonperturbative schemes [40][41][42][43][44][45][46][47][48][49][50] and multiplicative renormalizability in coordinate space has been proven to all orders in continuum perturbation theory [46,[51][52][53][54]. Similarly, it is expected that quark bilinear operators with staple-shaped Wilson lines are also multiplicatively renormalizable [23,30,55], such that they can be renormalized nonperturbatively via the regularization-independent momentum subtraction scheme (RI /MOM). The matching from RI /MOM to MS can then be calculated analytically in the continuum theory with dimensional regularization, which is free from linear divergences. The one-loop matching coefficient has been calculated for operators with zero longitudinal separation of the quark fields [55], which are relevant in the study of the x-moments of the TMDPDFs [19][20][21][22][23], and also for operators with quark fields separated longitudinally, which determine the xdependence of the TMDPDFs [30].
Here, the RI /MOM renormalization of quasi beam functions is studied numerically in quenched QCD with improved Wilson valence fermions. Due to the explicit breaking of Lorentz and chiral symmetries in the calculation, the multiplicatively-renormalizable quark bilinear arXiv:1911.00800v1 [hep-lat] 3 Nov 2019 operators with straight or staple-shaped Wilson lines mix with others with different Dirac structures [43,46,55,56]; the complete 16×16 operator mixing matrix for stapleshaped operators with all possible Dirac structures is therefore computed here for the first time. This study is undertaken at lattice spacings of 0.04, 0.06, and 0.08 fm, and with a single lattice volume, L ∼ 2 fm. This enables a detailed analysis of the lattice-spacing dependence of the mixing patterns induced in the lattice theory.
To increase the signal-to-noise ratio of the calculation, the renormalization factors are computed using gauge field configurations for which the gauge links have been subject to 100 steps of Wilson flow to flow-time t = 1.0 [57]. With the flow-time fixed in lattice units, Wilson flow corresponds to a smearing prescription whose effects vanish in the continuum limit [58]. A subset of the calculations were also repeated without flow applied to the gauge fields. Typically, this smearing prescription significantly reduces mixing between different operator structures. This modifies the mixing patterns such that the dominant mixings are not necessarily those predicted by one-loop lattice perturbation theory with the unflowed action. Calculations are predominantly undertaken with a light quark mass corresponding to a pion mass of m π ∼ 1.2 GeV. On the coarsest lattice, calculations with m π ∼ 340 MeV are also undertaken. While naively one might expect the heavy quark mass to enhance operator mixing due to the chiral-symmetry breaking of the mass terms in the fermion action, little massdependence is observed in the results.
Finally, the subset of the nonperturbative RI /MOM renormalization factors required to calculate renormalized quark-bilinear operators with Dirac structure γ 4 are matched to the MS scheme using matching coefficients computed in one-loop perturbation theory [30], and their lattice-spacing dependence is studied. This completes the nonperturbative renormalization prescription for quasi beam functions needed to study TMDPDFs from lattice QCD.

II. QUASI TMDPDFS
TMDPDFs, which are relevant to scattering processes such as Drell-Yan and SIDIS, can be expressed in terms of beam functions (also referred to as unsubtracted TMD-PDFs) which describe the incoming colinear partons in the scattering process, and soft functions, which encode the effects of soft gluon radiation by partons. Beam functions are defined as hadron matrix elements of quark bilinear operators with staple-shaped Wilson lines extended along the light-cone direction, while the soft functions are defined as the vacuum matrix elements of Wilson loops extended along the incoming and outgoing light-cone directions. Since they are defined on the lightcone, neither the beam nor soft functions can be directly calculated with lattice QCD formulated in Euclidean space. Constraints on TMDPDFs from lattice QCD, however, are possible via the LaMET approach [24,25]. The principle of LaMET is to approximate light-cone parton distributions by static quasi distributions, defined in terms of Euclidean matrix elements which can be calculated nonperturbatively in highly boosted hadron states. At large hadron momentum, quasi distributions are then matched to light-cone parton distributions perturbatively. To calculate TMDPDFs, quasi TMDPDFs have been constructed in terms of quasi beam and quasi soft functions [26][27][28][29]. Due to the complication of the quasi soft function, 1 the relation between quasi TMD-PDFs and TMDPDFs is expected to be nonperturbative; the explicit form of this relation was presented in Refs. [28,29]. It was also shown in those works that contributions from the soft sector, which do not depend on the hadron state, cancel in certain ratios of TMDPDFs and in the corresponding ratios of quasi TMDPDFs. For this reason, physical observables defined by ratios of TMDPDFs can be determined from lattice QCD calculations of quasi beam functions alone. For example, the Collins-Soper kernel can be obtained from ratios of quasi beam functions at different hadron momenta [28,29]. Similarly, ratios of the x-moments of TMDPDFs can also be determined with lattice QCD [19][20][21][22][23].
Precisely, quasi beam functions are calculated as matrix elements of quark bilinear operators with stapleshaped Wilson lines, in position space: Here, h(P ) denotes a boosted hadron state with fourmomentum P µ . The staple-shaped Wilson-line operator O q Γ (b µ , 0, η) is built as a bilinear of quark flavor q, with a Wilson line of staple length η in theẑ direction connecting endpoints separated by b µ = b z + b T , where T denotes a direction transverse toẑ: This operator, which is depicted graphically in Fig. 1, is constructed from spatial Wilson lines that are defined as Fourier transforms with respect to b z of the bare quasi beam function for different Dirac structures Γ ∈{I,γ µ ,γ 5 ,γ µ γ 5 ,σ µν } define quasi TMDPDFs with different spin structures. Including the quasi soft factor T < l a t e x i t s h a 1 _ b a s e 6 h T x 6 k / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S A T X x n W / n c L W 9 s 7 u X n G / d H B 4 d H x S P j 3 r 6 D h V D N s s F r H q B V S j 4 B L b h h u B v U Q h j Q K B 3 W B 6 v / C 7 T 6 g 0 j 2 X L z B L 0 I z q W P O S M G i s 1 W 8 N y x a 2 6 S 5 B N 4 u W k A j k a w / L X Y B S z N E J p m K B a 9 z 0 3 M X 5 G l e F M 4 L w 0 S D U m l E 3 p G P u W S h q h 9 r P l o X N y Z Z U R C W N l S x q y V H 9 P Z D T S e h Y F t j O i Z q L X v Y X 4 n 9 d P T V j z M y 6 T 1 K B k q 0 V h K o i J y e J r M u I K m R E z S y h T 3 N 5 K 2 I Q q y o z N p m R D 8 N Z f 3 i S d m 6 r n V r 3 m b a V e y + M o w g V c w j V 4 c A d 1 e I A G t I E B w j O 8 w p v z 6 L w 4 7 8 7 H q r X g 5 D P n 8 A f O 5 w + s p 4 z O < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 h T x 6 k / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S A T X x n W / n c L W 9 s 7 u X n G / d H B 4 d H x S P j 3 r 6 D h V D N s s F r H q B V S j 4 B L b h h u B v U Q h j Q K B 3 W B 6 v / C 7 T 6 g 0 j 2 X L z B L 0 I z q W P O S M G i s 1 W 8 N y x a 2 6 S 5 B N 4 u W k A j k a w / L X Y B S z N E J p m K B a 9 z 0 3 M X 5 G l e F M 4 L w 0 S D U m l E 3 p G P u W S h q h 9 r P l o X N y Z Z U R C W N l S x q y V H 9 P Z D T S e h Y F t j O i Z q L X v Y X 4 n 9 d P T V j z M y 6 T 1 K B k q 0 V h K o i J y e J r M u I K m R E z S y h T 3 N 5 K 2 I Q q y o z N p m R D 8 N Z f 3 i S d m 6 r n V r 3 m b a V e y + M o w g V c w j V 4 c A d 1 e I A G t I E B w j O 8 w p v z 6 L w 4 7 8 7 H q r X g 5 D P n 8 A f O 5 w + s p 4 z O < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 h T x 6 k / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S A T X x n W / n c L W 9 s 7 u X n G / d H B 4 d H x S P j 3 r 6 D h V D N s s F r H q B V S j 4 B L b h h u B v U Q h j Q K B 3 W B 6 v / C 7 T 6 g 0 j 2 X L z B L 0 I z q W P O S M G i s 1 W 8 N y x a 2 6 S 5 B N 4 u W k A j k a w / L X Y B S z N E J p m K B a 9 z 0 3 M X 5 G l e F M 4 L w 0 S D U m l E 3 p G P u W S h q h 9 r P l o X N y Z Z U R C W N l S x q y V H 9 P Z D T S e h Y F t j O i Z q L X v Y X 4 n 9 d P T V j z M y 6 T 1 K B k q 0 V h K o i J y e J r M u I K m R E z S y h T 3 N 5 K 2 I Q q y o z N p m R D 8 N Z f 3 i S d m 6 r n V r 3 m b a V e y + M o w g V c w j V 4 c A d 1 e I A G t I E B w j O 8 w p v z 6 L w 4 7 8 7 H q r X g 5 D P n 8 A f O 5 w + s p 4 z O < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 h T x 6 k / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S A T X x n W / n c L W 9 s 7 u X n G / d H B 4 d H x S P j 3 r 6 D h V D N s s F r H q B V S j 4 B L b h h u B v U Q h j Q K B 3 W B 6 v / C 7 T 6 g 0 j 2 X L z B L 0 I z q W P O S M G i s 1 W 8 N y x a 2 6 S 5 B N 4 u W k A j k a w / L X Y B S z N E J p m K B a 9 z 0 3 M X 5 G l e F M 4 L w 0 S D U m l E 3 p G P u W S h q h 9 r P l o X N y Z Z U R C W N l S x q y V H 9 P Z D T S e h Y F t j O i Z q L X v Y X 4 n 9 d P T V j z M y 6 T 1 K B k q 0 V h K o i J y e J r M u I K m R E z S y h T 3 N 5 K 2 I Q q y o z N p m R D 8 N Z f 3 i S d m 6 r n V r 3 m b a V e y + M o w g V c w j V 4 c A d 1 e I A G t I E B w j O 8 w p v z 6 L w 4 7 8 7 H q r X g 5 D P n 8 A f O 5 w + s p 4 z O < / l a t e x i t > z < l a t e x i t s h a 1 _ b a s e 6 4 = " J 9 J D d H F B s 2 k 0 q 4 n H K r k 0 8 d e d 1 l g = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E s M e C F 4 8 t 2 A 9 o Q 9 l s J + 3 a z S b s b o Q a + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 / Y 3 C p u l 3 Z 2 9 / Y P 7 P J h S 8 e p I r R J Y h 6 r T o g 1 5 U z S J j D g t J M o i k X I a T s c X c / 8 9 g N V m s X y D s Y J D Q Q e S B Y x g s F I P b s 8 u f d F e u 5 T w P 4 Q Q z a Z 9 u y K W 3 X n c F a J l 5 M K y t H o 2 V 9 + P y a p o B I I x 1 p 3 P T e B I M M K G O F 0 W v J T T R N M R n h A u 4 Z K L K g O s v n p U + f U K H 0 n i p U p C c 5 c / T 2 R Y a H 1 W I S m U 2 A Y 6 m V v J v 7 n d V O I a k H G Z J I C l W S x K E q 5 A 7 E z y 8 H p M 0 U J 8 L E h m C h m b n X I E C t M w K R V M i F 4 y y + v k t Z F 1 X O r 3 u 1 l p V 7 L 4 y i i Y 3 S C z p C H r l A d 3 a A G a i K C H t E z e k V v 1 p P 1 Y r 1 b H 4 v W g p X P H K E / s D 5 / A J i 5 l C 0 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / Y 3 C p u l 3 Z 2 9 / Y P 7 P J h S 8 e p I r R J Y h 6 r T o g 1 5 U z S J j D g t J M o i k X I a T s c X c / 8 9 g N V m s X y D s Y J D Q Q e S B Y x g s F I P b s 8 u f d F e u 5 T w P 4 Q Q z a Z 9 u y K W 3 X n c F a J l 5 M K y t H o 2 V 9 + P y a p o B I I x 1 p 3 P T e B I M M K G O F 0 W v J T T R N M R n h A u 4 Z K L K g O s v n p U + f U K H 0 n i p U p C c 5 c / T 2 R Y a H 1 W I S m U 2 A Y 6 m V v J v 7 n d V O I a k H G Z J I C l W S x K E q 5 A 7 E z y 8 H p M 0 U J 8 L E h m C h m b n X I E C t M w K R V M i F 4 y y + v k t Z F 1 X O r 3 u 1 l p V 7 L 4 y i i Y 3 S C z p C H r l A d 3 a A G a i K C H t E z e k V v 1 p P 1 Y r 1 b H 4 v W g p X P H K E / s D 5 / A J i 5 l C 0 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / Y 3 C p u l 3 Z 2 9 / Y P 7 P J h S 8 e p I r R J Y h 6 r T o g 1 5 U z S J j D g t J M o i k X I a T s c X c / 8 9 g N V m s X y D s Y J D Q Q e S B Y x g s F I P b s 8 u f d F e u 5 T w P 4 Q Q z a Z 9 u y K W 3 X n c F a J l 5 M K y t H o 2 V 9 + P y a p o B I I x 1 p 3 P T e B I M M K G O F 0 W v J T T R N M R n h A u 4 Z K L K g O s v n p U + f U K H 0 n i p U p C c 5 c / T 2 R Y a H 1 W I S m U 2 A Y 6 m V v J v 7 n d V O I a k H G Z J I C l W S x K E q 5 A 7 E z y 8 H p M 0 U J 8 L E h m C h m b n X I E C t M w K R V M i F 4 y y + v k t Z F 1 X O r 3 u 1 l p V 7 L 4 y i i Y 3 S C z p C H r l A d 3 a A G a i K C H t E z e k V v 1 p P 1 Y r 1 b H 4 v W g p X P H K E / s D 5 / A J i 5 l C 0 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / Y 3 C p u l 3 Z 2 9 / Y P 7 P J h S 8 e p I r R J Y h 6 r T o g 1 5 U z S J j D g t J M o i k X I a T s c X c / 8 9 g N V m s X y D s Y J D Q Q e S B Y x g s F I P b s 8 u f d F e u 5 T w P 4 Q Q z a Z 9 u y K W 3 X n c F a J l 5 M K y t H o 2 V 9 + P y a p o B I I x 1 p 3 P T e B I M M

l a t e x i t s h a 1 _ b a s e 6 4 = " O O h Y N I 0 n Y 2 s o h p f 2 V H L c S i 8 a j P 4 = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E s M e C F 4 8 t 9 A v a U D b b S b t 2 s w m 7 G 6 G E / g I v H h T x 6 k / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S A T X x n W / n c L W 9 s 7 u X n G / d H B 4 d H x S P j 3 r 6 D h V D N s s F r H q B V S j 4 B L b h h u B v U Q h j Q K B 3 W B 6 v / C 7 T 6 g 0 j 2 X L z B L 0 I z q W P O S M G i s 1 W 8 N y x a 2 6 S 5 B N 4 u W k A j k a w / L X Y B S z N E J p m K B a 9 z 0 3 M X 5 G l e F M 4 L w 0 S D U m l E 3 p G P u W S h q h 9 r P l o X N y Z Z U R C W N l S x q y V H 9 P Z D T S e h Y F t j O i Z q L X v Y X 4 n 9 d P T V j z M y 6 T 1 K B k q 0 V h K o i J y e J r M u I K m R E z S y h T 3 N 5 K 2 I Q q y o z N p m R D 8 N Z f 3 i S d m 6 r n V r 3 m b a V e y + M o w g V c w j V 4 c A d 1 e I A G t I E B w j O 8 w p v z 6 L w 4 7 8 7 H q r X g 5 D P n 8 A f O 5 w + s p 4 z O < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O O h Y N I 0 n Y 2 s o h p f 2 V H L c S i 8 a j P 4 = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E s M e C F 4 8 t 9 A v a U D b b S b t 2 s w m 7 G 6 G E / g I v H h T x 6 k / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S A T X x n W / n c L W 9 s 7 u X n G / d H B 4 d H x S P j 3 r 6 D h V D N s s F r H q B V S j 4 B L b h h u B v U Q h j Q K B 3 W B
6 v / C 7 T 6 g 0 j 2 X L z B L 0 I z q W P O S M G i s 1 W 8 N y x a 2 6 S 5 B N 4 u W k A j k a w / L X Y B S z N E J p m K B a 9 z 0 3 M X 5 G l e F M 4 L w 0 S D U m l E 3 p G P u W S h q h 9 r P l o X N y Z Z U R C W N l S x q y V H 9 P Z D T S e h Y F t j O i Z q L X v Y X 4 n 9 d P T V j z M y 6 T 1 K B k q 0 V h K o i J y e J r M u I K m R E z S y h T 3 N 5 K 2 I Q q y o z N p m R D 8 N Z f 3 i S d m 6 r n V r 3 m b a V e y + M o w g V c w j V 4 c A d 1 e I A G t I E B w j O 8 w p v z 6 L w 4 7 8 7 H q r X g 5 D P n 8 A f O 5 w + s p 4 z O < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 6 v / C 7 T 6 g 0 j 2 X L z B L 0 I z q W P O S M G i s 1 W 8 N y x a 2 6 S 5 B N 4 u W k A j k a w / L X Y B S z N E J p m K B a 9 z 0 3 M X 5 G l e F M 4 L w 0 S D U m l E 3 p G P u W S h q h 9 r P l o X N y Z Z U R C W N l S x q y V H 9 P Z D T S e h Y F t j O i Z q L X v Y X 4 n 9 d P T V j z M y 6 T 1 K B k q 0 V h K o i J y e J r M u I K m R E z S y h T 3 N 5 K 2 I Q q y o z N p m R D 8 N Z f 3 i S d m 6 r n V r 3 m b a V e y + M o w g V c w j V 4 c A d 1 e I A G t I E B w j O 8 w p v z 6 L w 4 7 8 7 H q r X g 5 D P n 8 A f O 5 w + s p 4 z O < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 s Y I D z D K 7 w 5 D 8 6 L 8 + 5 8 L F s L T j 5 z C n / g f P 4 A 5 j + M 9 A = = < / l a t e x i t > z µ + ⌘ẑ < l a t e x i t s h a 1 _ b a s e 6 4 = " l 7 B Q G Q e A C p h D / Y 3 C p u l 3 Z 2 9 / Y P 7 P J h S 8 e p I r R J Y h 6 r T o g 1 5 U z S J j D g t J M o i k X I a T s c X c / 8 9 g N V m s X y D s Y J D Q Q e S B Y x g s F I P b s 8 u f d F e u 5 T w P 4 Q Q z a Z 9 u y K W 3 X n c F a J l 5 M K y t H o 2 V 9 + P y a p o B I I x 1 p 3 P T e B I M M K G O F 0 W v J T T R N M R n h A u 4 Z K L K g O s v n p U + f U K H 0 n i p U p C c 5 c / T 2 R Y a H 1 W I S m U 2 A Y 6 m V v J v 7 n d V O I a k H G Z J I C l W S x K E q 5 A 7 E z y 8 H p M 0 U J 8 L E h m C h m b n X I E C t M w K R V M i F 4 y y + v k t Z F 1 X O r 3 u 1 l p V 7 L 4 y i i Y 3 S C z p C H r l A d 3 a A G a i K C H t E z e k V v 1 p P 1 Y r 1 b H 4 v W g p X P H K E / s D 5 / A J i 5 l C 0 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " l 7 B Q G Q e A C p h D x U i i 4 + o U O Z a 8 O u Y = " > A A A B + n i c b V B N S 8 N A E N 3 U r 1 q / U j 1 6 C R Z B E E o i g j 0 W v H i s Y D + g i W W z 3 b R L d z d h d 6 K 0 s T / F i w d F v P p L v P l v 3 L Y 5 a O u D g c d 7 M 8 z M C x P O N L j u t 1 V Y W 9 / Y 3 C p u l 3 Z 2 9 / Y P 7 P J h S 8 e p I r R J Y h 6 r T o g 1 5 U z S J j D g t J M o i k X I a T s c X c / 8 9 g N V m s X y D s Y J D Q Q e S B Y x g s F I P b s 8 u f d F e u 5 T w P 4 Q Q z a Z 9 u y K W 3 X n c F a J l 5 M K y t H o 2 V 9 + P y a p o B I I x 1 p 3 P T e B I M M K G O F 0 W v J T T R N M R n h A u 4 Z K L K g O s v n p U + f U K H 0 n i p U p C c 5 c / T 2 R Y a H 1 W I S m U 2 A Y 6 m V v J v 7 n d V O I a k H G Z J I C l W S x K E q 5 A 7 E z y 8 H p M 0 U J 8 L E h m C h m b n X I E C t M w K R V M i F 4 y y + v k t Z F 1 X O r 3 u 1 l p V 7 L 4 y i i Y 3 S C z p C H r l A d 3 a A G a i K C H t E z e k V v 1 p P 1 Y r 1 b H 4 v W g p X P H K E / s D 5 / A J i 5 l C 0 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " l 7 B Q G Q e A C p h D x U i i 4 + o U O Z a 8 O u Y = " > A A A B + n i c b V B N S 8 N A E N 3 U r 1 q / U j 1 6 C R Z B E E o i g j 0 W v H i s Y D + g i W W z 3 b R L d z d h d 6 K 0 s T / F i w d F v P p L v P l v 3 L Y 5 a O u D g c d 7 M 8 z M C x P O N L j u t 1 V Y W 9 / Y 3 C p u l 3 Z 2 9 / Y P 7 P J h S 8 e p I r R J Y h 6 r T o g 1 5 U z S J j D g t J M o i k X I a T s c X c / 8 9 g N V m s X y D s Y J D Q Q e S B Y x g s F I P b s 8 u f d F e u 5 T w P 4 Q Q z a Z 9 u y K W 3 X n c F a J l 5 M K y t H o 2 V 9 + P y a p o B I I x 1 p 3 P T e B I M M K G O F 0 W v J T T R N M R n h A u 4 Z K L K g O s v n p U + f U K H 0 n i p U p C c 5 c / T 2 R Y a H 1 W I S m U 2 A Y 6 m V v J v 7 n d V O I a k H G Z J I C l W S x K E q 5 A 7 E z y 8 H p M 0 U J 8 L E h m C h m b n X I E C t M w K R V M i F 4 y y + v k t Z F 1 X O r 3 u 1 l p V 7 L 4 y i i Y 3 S C z p C H r l A d 3 a A G a i K C H t E z e k V v 1 p P 1 Y r 1 b H 4 v W g p X P H K E / s D 5 / A J i 5 l C 0 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " l 7 B Q G Q e A C p h D x U i i 4 + o U O Z a 8 O u Y = " > A A A B + n i c b V B N S 8 N A E N 3 U r 1 q / U j 1 6 C R Z B E E o i g j 0 W v H i s Y D + g i W W z 3 b R L d z d h d 6 K 0 s T / F i w d F v P p L v P l v 3 L Y 5 a O u D g c d 7 M 8 z M C x P O N L j u t 1 V Y W 9 / Y 3 C p u l 3 Z 2 9 / Y P 7 P J h S 8 e p I r R J Y h 6 r T o g 1 5 U z S J j D g t J M o i k X I a T s c X c / 8 9 g N V m s X y D s Y J D Q Q e S B Y x g s F I P b s 8 u f d F e u 5 T w P 4 Q Q z a Z 9 u y K W 3 X n c F a J l 5 M K y t H o 2 V 9 + P y a p o B I I x 1 p 3 P T e B I M M K G O F 0 W v J T T R N M R n h A u 4 Z K L K g O s v n p U + f U K H 0 n i p U p C c 5 c / T 2 R Y a H 1 W I S m U 2 A Y 6 m V v J v 7 n d V O I a k H G Z J I C l W S x K E q 5 A 7 E z y 8 H p M 0 U J 8 L E h m C h m b n X I E C t M w K R V M i F 4 y y + v k t Z F 1 X O r 3 u 1 l p V 7 L 4 y i i Y 3 S C z p C H r l A d 3 a A G a i K C H t E z e k V v 1 p P 1 Y r 1 b H 4 v W g p X P H K E / s D 5 / A J i 5 l C 0 = < / l a t e x i t > q(z µ + b µ ) < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 c c A E j Q Z C V K o p Q 9 C H 3 Q + l 7 6 s t e U = " > A A A C A n i c b V B L S w M x G M z W V 6 2 v V U / i J V i E i l B 2 i 6 D H g h e P F e w D u m v J p t k 2 N M m u S V a o S / H i X / H i Q R G v / g p v / h u z 7 R 6 0 d S B h m P m G 5 J s g Z l R p x / m 2 C k v L K 6 t r x f X S x u b W 9 o 6 9 u 9 d S U S I x a e K I R b I T I E U Y F a S p q W a k E 0 u C e M B I O x h d Z n 7 7 n k h F I 3 G j x z H x O R o I G l K M t J F 6 9 o E X G T t L p 3 e T y s O t x 5 P T I L t P e n b Z q T p T w E X i 5 q Q M c j R 6 9 p f X j 3 D C i d C Y I a W 6 r h N r P 0 V S U 8 z I p O Q l i s Q I j 9 C A d A 0 V i B P l p 9 M V J v D Y K H 0 Y R t I c o e F U / Z 1 I E V d q z A M z y Z E e q n k v E / / z u o k O L / y U i j j R R O D Z Q 2 H C o I 5 g 1 g f s U 0 m w Z m N D E J b U / B X i I Z I I a 9 N a y Z T g z q + 8 S F q 1 q u t U 3 e u z c r 2 W 1 1 E E h + A I V I A L z k E d X I E G a A I M H s E z e A V v 1 p P 1 Y r 1 b H 7 P R g p V n 9 s E f W J 8 / U 5 O X T g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 c c A E j Q Z C V K o p Q 9 C H 3 Q + l 7 6 s t e U = " > A A A C A n i c b V B L S w M x G M z W V 6 2 v V U / i J V i E i l B 2 i 6 D H g h e P F e w D u m v J p t k 2 N M m u S V a o S / H i X / H i Q R G v / g p v / h u z 7 R 6 0 d S B h m P m G 5 J s g Z l R p x / m 2 C k v L K 6 t r x f X S x u b W 9 o 6 9 u 9 d S U S I x a e K I R b I T I E U Y F a S p q W a k E 0 u C e M B I O x h d Z n 7 7 n k h F I 3 G j x z H x O R o I G l K M t J F 6 9 o E X G T t L p 3 e T y s O t x 5 P T I L t P e n b Z q T p T w E X i 5 q Q M c j R 6 9 p f X j 3 D C i d C Y I a W 6 r h N r P 0 V S U 8 z I p O Q l i s Q I j 9 C A d A 0 V i B P l p 9 M V J v D Y K H 0 Y R t I c o e F U / Z 1 I E V d q z A M z y Z E e q n k v E / / z u o k O L / y U i j j R R O D Z Q 2 H C o I 5 g 1 g f s U 0 m w Z m N D E J b U / B X i I Z I I a 9 N a y Z T g z q + 8 S F q 1 q u t U 3 e u z c r 2 W 1 1 E E h + A I V I A L z k E d X I E G a A I M H s E z e A V v 1 p P 1 Y r 1 b H 7 P R g p V n 9 s E f W J 8 / U 5 O X T g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 c c A E j Q Z C V K o p Q 9 C H 3 Q + l 7 6 s t e U = " > A A A C A n i c b V B L S w M x G M z W V 6 2 v V U / i J V i E i l B 2 i 6 D H g h e P F e w D u m v J p t k 2 N M m u S V a o S / H i X / H i Q R G v / g p v / h u z 7 R 6 0 d S B h m P m G 5 J s g Z l R p x / m 2 C k v L K 6 t r x f X S x u b W 9 o 6 9 u 9 d S U S I x a e K I R b I T I E U Y F a S p q W a k E 0 u C e M B I O x h d Z n 7 7 n k h F I 3 G j x z H x O R o I G l K M t J F 6 9 o E X G T t L p 3 e T y s O t x 5 P T I L t P e n b Z q T p T w E X i 5 q Q M c j R 6 9 p f X j 3 D C i d C Y I a W 6 r h N r P 0 V S U 8 z I p O Q l i s Q I j 9 C A d A 0 V i B P l p 9 M V J v D Y K H 0 Y R t I c o e F U / Z 1 I E V d q z A M z y Z E e q n k v E / / z u o k O L / y U i j j R R O D Z Q 2 H C o I 5 g 1 g f s U 0 m w Z m N D E J b U / B X i I Z I I a 9 N a y Z T g z q + 8 S F q 1 q u t U 3 e u z c r 2 W 1 1 E E h + A I V I A L z k E d X I E G a A I M H s E z e A V v 1 p P 1 Y r 1 b H 7 P R g p V n 9 s E f W J 8 / U 5 O X T g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 c c A E j Q Z C V K o p Q 9 C H 3 Q + l 7 6 s t e U = " > A A A C A n i c b V B L S w M x G M z W V 6 2 v V U / i J V i E i l B 2 i 6 D H g h e P F e w D u m v J p t k 2 N M m u S V a o S / H i X / H i Q R G v / g p v / h u z 7 R 6 0 d S B h m P m G 5 J s g Z l R p x / m 2 C k v L K 6 t r x f X S x u b W 9 o 6 9 u 9 d S U S I x a e K I R b I T I E U Y F a S p q W a k E 0 u C e M B I O x h d Z n 7 7 n k h F I 3 G j x z H x O R o I G l K M t J F 6 9 o E X G T t L p 3 e T y s O t x 5 P T I L t P e n b Z q T p T w E X i 5 q Q M c j R 6 9 p f X j 3 D C i d C Y I a W 6 r h N r P 0 V S U 8 z I p O Q l i s Q I j 9 C A d A 0 V i B P l p 9 M V J v D Y K H 0 Y R t I c o e F U / Z 1 I E V d q z A M z y Z E e q n k v E / / z u o k O L / y U i j j R R O D Z Q 2 H C o I 5 g 1 g f s U 0 m w Z m N D E J b U / B X i I Z I I a 9 N a y Z T g z q + 8 S F q 1 q u t U 3 e u z c r 2 W 1 1 E E h + A I V I A L z k E d X I E G a A I M H s E z e A V v 1 p P 1 Y r 1 b H 7 P R g p V n 9 s E f W J 8 / U 5 O X T g = = < / l a t e x i t > q(z µ ) < l a t e x i t s h a 1 _ b a s e 6 4 = " r b 1 x w f C m 1 J X r S 3 V k x B b 4 u t l r X l Q = " > A A A B 7 3 i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C R a h X s p u K e i x 4 M V j B f s B 7 V q y a b Y N T b L b J C v U p X / C i w d F v P p 3 v P l v T N s 9 a O u D g c d 7 M 8 z M C 2 L O t H H d b y e 3 s b m 1 v Z P f L e z t H x w e F Y 9 P W j p K F K F N E v F I d Q K s K W e S N g 0 z n H Z i R b E I O G 0 H 4 5 u 5 3 3 6 k S r N I 3 p t p T H 2 B h 5 K F j G B j p c 6 k / P T Q E 8 l l v 1 h y K + 4 C a J 1 4 G S l B h k a / + N U b R C Q R V B r C s d Z d z 4 2 N n 2 J l G O F 0 V u g l m s a Y j P G Q d i 2 V W F D t p 4 t 7 Z + j C K g M U R s q W N G i h / p 5 I s d B    ∆ q S [26][27][28][29], the quasi TMDPDF in the MS scheme is defined as where Z MS ΓΓ (µ, b z ) renormalizes the quasi TMDPDF and matches it to the MS scheme at scale µ. The details of the definition and other properties of the quasi soft factor∆ q S are omitted here; since it only depends on b T , it will cancel in ratios of quasi TMDPDFs, as formed in key applications including calculations to extract the Collins-Soper kernel [19-23, 28, 29].
Both the quasi beam function and soft factor have linear UV divergences that are proportional to the total length of the Wilson line in the associated operator, and as such the quasi soft factor also acts as a counterterm to cancel the linear divergences in the bare quasi beam function. Nevertheless, there is still a remaining linear divergence proportional to |b z |, as well as other logarithmic divergences, which are renormalized by Z MS ΓΓ . The renormalization factor Z MS ΓΓ can be separated into pieces which renormalize the quasi beam function and soft factor, Z MS O ΓΓ and Z MS S , respectively: Since Z MS S (µ, b T , η) is independent of b z , it will also be canceled in ratios of quasi TMDPDFs at η and b T , thus leaving only the quasi beam functions to be renormalized for key applications. Taking into account mixing among O q Γ (b µ , z µ , η) with different Dirac structures, Z MS O ΓΓ is a 16×16 matrix that can be computed nonperturbatively via the RI /MOM prescription [59,60] with a perturbative matching to the MS scheme [30], as detailed in the next section.

III. NONPERTURBATIVE RENORMALIZATION
The bare staple-shaped Wilson line operator, Eq. (2), and hence the bare quasi beam function, can be renor-malized via the nonperturbative RI /MOM prescription [59,60]. In this approach, a renormalization constant is defined to relate the bare and tree-level amputated Green's functions for a given operator in a gaugefixed quark or gluon state at a fixed scale. A perturbative matching calculated in continuum perturbation theory then relates the resulting RI /MOM renormalized operator to the MS scheme. For a lattice operator O latt Γ , which implicitly depends on the staple extent η, the displacement between the staple endpoints b µ , and the lattice spacing a, this renormalization can be expressed as a matrix equation accounting for mixing of operators with different Dirac structures Γ: where p R denotes the matching scale introduced in the RI /MOM scheme. The determination of the nonperturbative RI /MOM renormalization matrix Z O ΓΓ (µ, p R ) is calculated at finite orders in perturbation theory; for the operators considered here, only one-loop results have been computed [30], so the cancellation is incomplete. Moreover, at finite a there are lattice artifacts that have p R dependence. In the pseudoscalar case (Γ = γ 5 ), Z RI/MOM O γ 5 Γ (p R , a) additionally develops a nonperturbative Goldstone boson pole that depends on p R [61].
For the operators considered here, the one-loop corrections in the diagonal terms R MS OΓΓ (µ, p R ) (i.e., with Γ = Γ ) are significantly larger than one [30], indicating that the perturbative series does not converge well. These large one-loop corrections can be cancelled by combining R MS O ΓΓ (µ, p R ) with the quasi soft factor∆ q S in the MS scheme, thus rendering the matching coefficient close to one [30] without affecting ratios of quasi TMDPDFs or the extraction of key physics results determined by such ratios. In the analysis presented here, higher-order perturbative corrections are neglected, and remnant p R dependence is treated as a discretization effect leading to systematic uncertainty in the results discussed in Sec. IV.
relating the bare and tree-level values of the operator's amputated Green's function in an off-shell quark state in the Landau gauge: where G OΓ denotes the Green's function for operator O Γ with Dirac structure Γ, which implicitly depends on the staple extent η and displacement between staple endpoints b µ , and S(p) is the quark propagator projected to momentum p. All quantities appearing on the left-handside of Eq. (7) implicitly depend on the lattice spacing; this dependence is suppressed in the following discussion. In Eq. (7), p 2 R acts a non-perturbative renormalization scale; however, since the operator O q Γ is nonlocal and frame dependent, the magnitude of p µ R alone is not sufficient to specify the renormalization condition. Different directions in p µ R amount to different renormalization schemes, which are related by finite renormalization factors. As a result, Z RI /MOM O ΓΓ (p R ) depends on p µ R rather than only its magnitude.
In a calculation with lattice volume V = L 3 ×T and lattice spacing a, the non-amputated quark-quark Green's function with one insertion of the operator O Γ is calculated as using the quark propagator The quark wavefunction renormalization Z q is defined via computed as Z q (p R ) = Tr i λ γ λ sin(ap λ )S −1 (p) 12 λ sin 2 (ap λ ) In terms of the projected vertex function the RI /MOM condition in Eq. (7), for an operator O Γ with endpoints separated by b µ , can be expressed as which yields an expression for the matrix of renormalization factors at p R : (18)

B. Conversion to the MS scheme
Since the renormalized matrix element in the RI /MOM scheme is independent of the UV regulator, it differs from the result in the continuum limit only by discretization effects at finite lattice spacing. The RI /MOM matrix element can thus be matched to the MS scheme in continuum perturbation theory, and then extrapolated to the continuum limit using nonperturbative calculations at different values of a.
Elements of the matrix of matching coefficients R MS O ΓΓ (µ, p R ) in Eq. (6) have been calculated at one-loop order in continuum perturbation theory with dimensional regularization (D = 4 − 2 ) for operators O Γ with both b z = 0 [55] and b z = 0 [30]. This matching matrix can be expressed as whereZ RI /MOM O;ΓΓ (p R , µ, ) is the perturbatively-computed RI /MOM renormalization factor for the quasi beam function, defined in Eq. (18). The factor Z MS ( , µ) is gauge-invariant and universal for all Dirac structures Γ [30,55]: where c f = 4/3. For Γ = γ λ , the matching coefficient R MS γ λ Γ has been calculated for all projectors Γ at one-loop order [30]. The results are summarized here for completeness: where Z (1) q (p R , µ) = 0 in the Landau gauge, and the subtraction of 1/ poles is implied. The explicit expression for the one-loop projected vertex functions V (1) γ λ ,γ ρ /γ ρ γ5 can be found in Ref. [30].
Defined in this way, the numerical values of the matching coefficients for the parameters of typical lattice QCD studies are much larger than one, which is due to η/b T terms that correspond to the rapidity divergences in the TMDPDF [27,29]. In Eq. (4), the quasi TMDPDF is defined with a quasi soft factor∆ q S = 1/ S q which cancels the linear power divergences as well as the η/b T dependence in the quasi beam function; redefining the matching coefficient to include the quasi soft factor removes this divergence and yields a matching coefficient close to one 2 [30]: In the numerical investigation presented in the following section, the "bent" quasi soft factor defined in Ref. [29,30] is adopted for this redefinition: Since the quasi soft factor only depends on the operator staple geometry in terms of b T and η, its inclusion will not change the p R or b z -dependence of the matching coefficient and therefore will not affect results for ratios of MS quasi beam functions at fixed b T and η.

IV. NUMERICAL INVESTIGATION
The RI /MOM renormalization of quark bilinear operators with staple-shaped Wilson lines is studied on three quenched QCD ensembles, detailed in Table I. These ensembles are tuned to have lattice spacings of 0.04, 0.06, and 0.08 fm, and a common lattice volume, L ∼ 2 fm. This enables study of the lattice-spacing dependence of renormalization factors and operator mixing patterns in 2 Note that unlike the proposal in Ref. [30], here the redefined matching coefficientR MS ΓΓ is not expanded as a series in αs. This maintains the numerical equivalence to MS matching when calculating ratios of quasi TMDPDFs. Moreover, with Eq. (24) expanded in αs as in Ref. [30], the quasi soft factor matching term only contributes to diagonal entriesR MS ΓΓ , which leads to significantly enhanced operator mixing in MS results that is not present in RI /MOM results.   [63,64]. β values were chosen in Ref. [65] to maintain a fixed physical volume, and n cfg con-  Table I), i.e., to almost half the lattice extent, and with staple widths and asymmetries b T and b z ranging from −η to η, for the complete 16×16 matrix of Dirac structures Γ, Γ . The gauge link fields used in the calculation have been subjected to Wilson flow to flow-time t = 1.0 [57], to enhance the signal-to-noise ratio in the numerical results; 3 a study of the impact of this smearing prescription on mixing patterns is given in App B. Valence quark propagators are computed with the tree-level O(a) improved Wilson clover fermion action [62] with κ values as given in Table I; these choices correspond to a pion mass of 1.2 GeV on each ensemble. For the E 24 ensemble, propagators corresponding to a pion mass of 340 MeV are also computed to enable a study of the mass-dependence of the renormalization patterns. Ten different momenta of the quark state are considered, tabulated in Tab. II. While the dependence of the RI /MOM renormalization on the matching scales p R and p z R would be cancelled by an allorders matching to the MS scheme, residual dependence on these scales remains with a matching calculated perturbatively to one-loop order. Studying various momenta at a range of scales p 2 R from 5.7 to 28 GeV 2 and p z R from 1.3 to 2.6 GeV allows an assessment of the systematic uncertainties in this matching.    mixing patterns revealed in the matrix of nonperturbative RI /MOM factors with the patterns predicted by perturbation theory, which have been studied in the special cases of local operators, straight Wilson-line operators, and symmetric staple-shaped Wilson line operators, provides an indication of the important nonperturbative mixings for each operator.
Figs. 2-7 display graphically the 16 × 16 matrices of RI /MOM renormalization factors for all Dirac structures Γ and projectors P, for a range of operators with different staple widths and asymmetries b T and b z , defined in Eq. (2). In each case, percentage mixings relative to the average diagonal element are displayed, defined as: where to illustrate the importance of mixings the maximum over momenta p R is taken over the ten momenta tabulated in Tab. II. Due to the off-shell nature of the quark in the Green's functions and the noncovariance of the operator O q Γ (b µ , 0, η) itself, there can be contributions from additional Dirac structures involving p µ R and b µ to the vertex function of O q Γ (b µ , 0, η), which do not break chiral symmetry and are also seen in continuum perturbation theory [59]. In lattice calculations, due to the breaking of chiral symmetry from the UV regularization, there are additional operator mixings that were predicted by one-loop lattice perturbation theory and symmetry arguments [43,46,55,56]. In addition, the interplay of these two mechanisms can generate new chiral-symmetry-breaking mixings, for example between γ 0 and 1. On each figure (except for Figs. 3 and 7), the chiral-symmetry-breaking mixing patterns predicted by one-loop lattice perturbation theory or symmetry arguments are highlighted for comparison with the numerical results.
In general, operators with longer Wilson lines are seen to suffer from greater mixing effects than operators with shorter Wilson lines; this is shown explicitly for the straight Wilson line operators in Fig. 3. Typically, the mixings predicted by lattice perturbation theory are found to be significant nonperturbatively, but in many cases other chiral-symmetry-preserving mixings are found to be equally, or more, important. The patterns of mixings computed on the three ensembles with different lattice spacings are consistent for each operator shape, with the relative magnitude of off-diagonal mixings relatively larger on the finer ensembles, as shown for the straight Wilson line case in Fig. 5.
A subset of calculations on the E 24 ensemble were repeated without Wilson flow applied to the gauge fields or Dirac operator. As outlined in Appendix B, Wilson flow generically reduces operator mixing, and in partic- From left to right, panels show results for operators with extent b z /a = 11 calculated on the ensembles E24, E32, E48, with progressively finer lattice spacing a. White circles indicate the mixings predicted by one-loop lattice perturbation theory and symmetry arguments [43,46,56]. White circles indicate the mixings predicted by one-loop lattice perturbation theory [55]. ular reduces some off-diagonal elements of the renormalization matrix significantly more than others, such that one-loop lattice perturbation theory (with an unflowed action) describes the unflowed mixing pattern somewhat better than the flowed mixing pattern.

B. Renormalization results
The row of the MS renormalization matrix for bare quasi beam functions with operator Dirac structure Γ = γ 4 is sufficient to determine MS-renormalized matrix elements of O MS γ4 , given bare matrix elements O latt Γ for all 16 choices of Γ. These MS renormalization factors are defined from the nonperturbatively computed RI /MOM factor and the perturbative one-loop matching by where the left hand side is independent of the choice of p R up to discretization effects, nonperturbative effects that vanish at asymptotically large p 2 R , and neglected two-loop perturbative matching corrections. Here, Z MS Oγ 4 Γ (µ, p R ) implicitly includes the quasi soft factor included iñ R MS γ4Γ (µ, p R ) (and thus differs from Z MS Oγ 4 Γ (µ, p R ) defined in Eq. (5) by terms which cancel in suitable ratios of renormalized TMDPDFs), and both Z MS Oγ 4 Γ (µ, p R ) and Z RI /MOM O Γ Γ (p R ) implicitly depend on a. This renormalization factor is computed for each choice of Γ with each of the 10 p R shown in Table II, for staple-shaped operators with −η < b T < η, −η ≤ b z ≤ η, for three values of η on each ensemble shown in Table I.
To determine Z MS Oγ 4 Γ from numerical results at different choices of p R , one could fit the data to a model of the discretization effects in the renormalization matrix. However, statistical noise in the nonlocal operator renormalization grows exponentially with the length of the Wilson line; in the present study it is not possible to constrain discretization effects from the 10 momenta used for all but the smallest nonlocal operator separations. In particular Bayes and Akaike information criteria prefer constant fits to more complicated fit forms including the leading discretization artifacts in the data (the functional form of these effects is made explicit in Appendix A). Moreover, the covariance matrices for nonlocal operators are not reliably estimated from the current data.
Rather than performing uncorrelated fits to correlated results, weighted averages are used to remove residual p R (pR) and Z MS Oγ 4 γ 4 (µ, pR) for the E32 ensemble with η/a = 10, b z /a = 3, bT /a = −4, µ = 2 GeV, are displayed as orange circles and blue squares, respectively. Results at ten choices of pR given in Table II are  dependence, where the weights are chosen to sum to unity and to be proportional to the inverse variance of the result for each momentum: The central value of this weighted average is identical to the central value of an uncorrelated fit and ensures that the fit is constrained most heavily by the most precise data. The inverse variance of this weighted average is the average inverse variance of the data, while the inverse variance of an uncorrelated χ 2 -minimization fit is equal to the same quantity times the number of data points. Uncorrelated fits to correlated data therefore lead to a spurious reduction in the uncertainty of the fit result that is avoided by Eq. (28). The systematic uncertainty term in Eq. (28) is included to reflect the uncertainty arising from unresolved discretization and nonperturbative effects. The resulting systematic error on Z MS is < 15% in all cases; for all but the largest Wilson line extents the systematic uncertainty is 2%. Similar results hold  for Z MS Oγ 4 Γ with Γ = γ 4 apart from cases where Z MS Oγ 4 Γ is consistent with zero. Fig. 8 shows a representative example of this weighted averaging procedure for an asymmetric staple operator with η/a = 10, b T /a = 3, and b z /a = 4, computed on ensemble E 32 . For this example and in general, the MS renormalization constant Z MS Oγ 4 Γ (µ, p R ) is more consistent with a constant and has smaller systematic uncertainties in a weighted average than Z RI /MOM Oγ 4 Γ (p R ), which indicates that one-loop matching accounts for some of the p R -dependence of the bare vertex function. Results for operators with displacements in the x−z and y−z planes, where x and y are the directions transverse to the staple extent η, are fit independently and found to be consistent within uncertainties, and the renormalization constants for operators of different shapes are found to be relatively smooth functions of the staple geometry parameterized by b z , b T , and η. Samples of these results are shown for the diagonal renormalization constant Z MS Oγ 4 γ 4 (µ) in Fig. 9. Here and throughout, µ = 2 GeV is used as a reference scale. The off-diagonal terms Z MS Oγ 4 Γ (µ) with Γ = γ 4 describing operator mixing indicate that such mixing is a percent-level effect for operators with small Wilson lines, but grows to become a 5 − 10% effect for the largest Wilson lines studied. A representative set of off-diagonal mixing results are shown in Fig. 10.

□ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △
In order to study the quark mass dependence of Z MS , calculations on the E 24 ensemble are repeated using a second quark mass corresponding to m π ∼ 340 MeV. For all Z MS Oγ 4 Γ , results for m π ∼ 340 MeV are found to be consistent within uncertainties with those calculated using m π ∼ 1.2 GeV, as shown in Fig. 11. This suggests that the large quark mass used in this work does not significantly affect results for Z MS Oγ 4 Γ . Before averaging over momentum, statistically significant differences between m π ∼ 1.2 GeV and m π ∼ 340 MeV results can be seen at the smallest momenta considered here, which is consistent with expectations that renormalization factors include nonperturbative quark mass effects proportional to m q qq /p 4 that vanish at large momentum [66][67][68][69][70][71][72][73]. After averaging over momentum, results with m π ∼ 1.2 GeV and m π ∼ 340 MeV are consistent within combined   with fixed b T and b z has an approximately exponential dependence on η as shown in Fig. 12. This η-dependence should cancel corresponding η-dependence in bare matrix elements, resulting in approximately ηindependent MS renormalized matrix elements. Similar exponential dependence on the Wilson line extent is seen in the b z and b T dependence of Z MS Oγ 4 γ 4 in Fig. 9 for large b T and large −b z . At smaller values of b T and −b z , additional structure beyond simple exponential dependence on the Wilson line length is demonstrated by the curvature visible in Fig. 9.
Ratios of off-diagonal and diagonal elements of the renormalization matrices Z MS Oγ 4 Γ / Z MS Oγ 4 γ 4 are seen to have mild a-dependence, as shown in Fig. 14. This is consistent with general arguments that the a → 0 divergence structure of Z MS Oγ 4 Γ does not depend on Γ discussed in Ref. [19].

V. SUMMARY
In this work, the nonperturbative RI /MOM renormalization of staple-shaped Wilson line operators, as relevant to lattice QCD studies of transverse-momentumdependent parton distribution functions, is investigated for the first time. The renormalization factors are computed nonperturbatively for a basis of nonlocal quark bilinear operators with a wide range of transverse and longitudinal separations in quenched QCD with three different lattice spacings, namely 0.04, 0.06, and 0.08 fm, and a single lattice volume, L ∼ 2 fm. Renormalization factors are found to depend exponentially on the length of the Wilson line in lattice units, as expected from per-turbation theory, although additional dependence on the shape of the Wilson line is clearly visible. Quark mass dependence is found to be negligible compared to uncertainties from statistical noise and lattice artifacts, for quark masses corresponding to m π = {0.4, 1.2} GeV.
Mixing between nonlocal quark bilinears with different Dirac operator structures is observed to be larger than the corresponding mixing between local quark bilinears; this mixing can not be neglected in studies of nonlocal quark bilinears targeting precision better than the 10% level. These operator mixing effects show mild lattice spacing dependence, with the effects typically found to be larger both for finer discretization scales and for operators built from Wilson lines with longer staple extents. While the mixing patterns predicted by one-loop lattice perturbation theory are observed, additional mixing effects that are as, or even more, significant than the predicted mixings are also often present. The results of this work allow bare matrix elements for a basis of non-local quark bilinear operators with staple-shaped Wilson lines to be renormalized, with the mixing between operators with different Dirac structures fully accounted for. This completes a critical step towards the systematic extraction of TMDPDFs, and also TMD distribution amplitudes, from lattice QCD. considered here can be expressed as [73] Z MS Oγ 4 Γ (µ, p) = Z MS Oγ 4 Γ (µ) + c 1pz + c 2p 2 + c 3p 2 t + c 4p [4] p 2 + c 5p 2 ln(p 2 ) + d 1 p 2 + . . . , where p [4] = 4 µ=1 p 4 µ and the parameters c i are adependent constants that can be extracted from fits to numerical data. On the right-hand side of this expression, momenta have been replaced with the momentum variable that arises in a discrete Fourier transform of the lattice action, namelỹ p µ ≡ 1 a sin(ap µ ).
Local operator renormalization factors and Z q have additional symmetry constraints leading to c 1 = 0 and c 3 = 0 up to negligible symmetry-breaking effects from the different extent of the lattice space and time directions. This leads to the functional form Z MS q (µ, p) = Z MS q (µ) + c 2p 2 + c 4p [4] p 2 + c 5p 2 ln(p 2 ) + d 1 p 2 + . . . .

(A3)
Results for Z MS q (µ, p) are fit to Eq. (A3) for each ensemble. Fit results for c 4 are used to remove rotationally non-invariant lattice artifacts as Z MS q (µ, p) = Z MS q (µ, p) − c 4p [4] p 2 .
(A4) Fig. 15 shows numerical results forZ MS q (µ = 2 GeV, p) as well as the best fit to Eq. (A3) for each ensemble studied here. Results for c 2 , c 4 , and c 5 are consistent within uncertainties for all three ensembles, as expected. As 0) with extent η/a = 9, bT /a = 3, respectively. The bare quark mass is tuned separated for flowed and unflowed gauge fields in order to give mπ ∼ 1.2 GeV in both cases. White circles indicate the mixings predicted by one-loop lattice perturbation theory and symmetry arguments [43,46,55,56]. discussed in Sec. IV B, for the nonlocal operator renormalization constants Z MS Oγ 4 Γ (µ, p), lattice artifacts cannot be clearly resolved, and simple constant fits are preferred over fits to Eq. (A1) by information criteria. For quark bilinear operators with straight Wilson lines, computed without Wilson flow, the mixings predicted by one-loop lattice perturbation theory are also the largest mixings nonperturbatively. With Wilson flow, these mixings are reduced significantly and become smaller than mixings that are not predicted by one-loop lattice perturbation theory. For symmetric staple-shaped Wilson line operators (b z = 0) without Wilson flow, mixings between operators with Dirac structures Γ and Γ ∈ {Γ,ẑ /} dominate over those predicted by one-loop lattice perturbation theory [55]. With Wilson flow applied to the gauge fields, all mixings are significantly reduced. It will be interesting to see whether the flowed mixing patterns are postdicted by flowed one-loop lattice perturbation theory.