Collins-Soper Kernel for TMD Evolution from Lattice QCD

The Collins-Soper kernel relates transverse momentum-dependent parton distribution functions (TMDPDFs) at different energy scales. For small parton transverse momentum $q_T\sim \Lambda_\text{QCD}$, this kernel is non-perturbative and can only be determined through experiment or first-principles calculations. This work presents the first exploratory determination of the Collins-Soper kernel using the lattice formulation of Quantum Chromodynamics. In a quenched calculation, the $N_f=0$ kernel is determined at scales in the range 250 MeV $<q_T<2$ GeV, and an analysis of the remaining systematic uncertainties is undertaken.


I. INTRODUCTION
Understanding the structure of matter has been a defining goal of physics for centuries. In the modern context, a primary objective is imaging the threedimensional spatial and momentum structure of the proton, and of other hadrons. Some important aspects of this structure related to the transverse momentum of quarks and gluons in a hadron state are encoded in transverse-momentum-dependent parton distribution functions (TMDPDFs) [1][2][3]. Experimentally, these quantities can be constrained for the proton by Drell-Yan production and semi-inclusive deep inelastic scattering (SIDIS) of electrons o↵ protons; the best current constraints are achieved via global fits to experimental data [4][5][6][7][8][9][10][11][12][13][14][15], with improvements expected in the coming years from measurements at COMPASS [16], the Thomas Je↵erson National Accelerator Facility [17], RHIC [18], and an Electron-Ion Collider [19]. Additional experimental information from dihadron production in e + e collisions at Belle, and new ways of looking at hadrons inside jets [20,21], may also help constrain these fits.
Key to global fits of TMDPDFs is the ability to relate these distributions determined in di↵erent processes, including those at di↵erent scales. That is, for a TMDPDF f TMD i x,b T , µ 0 , ⇣ 0 , defined for a parton of flavor i with longitudinal momentum fraction x, transverse displacementb T (the Fourier conjugate of the transverse momentumq T ), virtuality scale µ 0 , and hadron momentum scale ⇣ 0 which is related to the hard scale of the scattering process, it is critical to understand its evolution to scales (µ, ⇣): where b T = |b T |. The first exponential in this equation governs the µ-evolution of the TMDPDF, which is perturbative for scales {µ 0 , µ} ⇤ QCD . The evolution in ⇣ governed by the second exponential, however, is encoded in the Collins-Soper kernel 1 i ⇣ (µ, b T ), which is inherently non-perturbative for q T ⇠ b 1 T ⇠ ⇤ QCD , even for µ ⇤ QCD . Experimentally, the Collins-Soper kernel can be extracted by simultaneous global fits with the TMD-PDF, and recent global analyses show some discrepancy in determinations of the kernel in the region q T  500 MeV [22]. It would greatly improve systematic control if the Collins-Soper kernel could be independently determined from first-principles QCD calculations, and taken as input for global fits of experimental data.
Since TMDPDFs are defined in terms of light-cone correlation functions, they are challenging to calculate directly in the lattice formulation of QCD on a discrete Euclidean spacetime, which is the only known systematically improvable first-principles approach to nonperturbative QCD. Nevertheless, e↵orts to calculate aspects of TMD physics from equal-time correlation functions in boosted hadron states have been made in Refs. [23][24][25][26][27], and the large-momentum e↵ective theory (LaMET) framework [28,29] provides a promising pathway towards the determination of TMDPDFs by matching these matrix elements to the desired light-cone correlation functions at large hadron momentum [30][31][32][33][34][35][36][37]. In particular, it was recently shown in Refs. [32,33] how this approach may be used to extract the Collins-Soper kernel nonperturbatively from computations of matrix elements of nonlocal quark bilinear operators with staple-shaped Wilson lines. Here, this approach is implemented numerically for the first time, in a proof-of-principle calculation in quenched QCD. The Collins-Soper kernel is extracted at a range of q T scales, including in the non-perturbative region.
using lattice QCD and LaMET. Section III details the quenched lattice QCD calculation undertaken here, including discussion of the systematic uncertainties in the calculation, while Sec. IV outlines the requirements for a fully-controlled calculation of the Collins-Soper kernel to be achieved by this method.

II. COLLINS-SOPER KERNEL FROM LATTICE QCD
In Refs. [32,33] a method was proposed to determine the quark Collins-Soper kernel using lattice QCD and LaMET. Precisely, it was shown that q ⇣ (µ, b T ) can be extracted from a ratio of nonsinglet quasi TMDPDFsf TMD ns at di↵erent momenta, which are defined using equal-time correlation functions within hadron states at large momentum in the z-direction: , (2) up to power corrections which are discussed further below. In this expression, P z i ⇤ QCD are the z-component of the hadron momenta and C TMD ns is a perturbative matching coe cient that has been obtained at one-loop order [32,33]. The quasi TMDPDFf TMD ns , defined below, approximates the physical TMDPDF involving light-like paths, as detailed in Ref. [33], and complications involving matching in the soft sector [31,33,35,36] are eliminated in the ratio that gives the Collins-Soper kernel. Similar constructions have been used in calculations of ratios of x-moments of TMDPDFs from lattice QCD [23][24][25][26][27].
The unpolarized quasi TMDPDF is defined in terms of a quasi beam functionB i and a quasi soft factor˜ S [30][31][32][33], both of which are calculable in lattice QCD: where a denotes the lattice spacing, the subscript i is the flavor index, and summation over Dirac structures is implied. This summation accounts for the operator mixings among di↵erent Dirac structures in lattice QCD calculations defined on a hypercubic space-time lattice [38][39][40].
Additional mixing with gluon operators, not shown in Eq. (3), cancels in the flavor nonsinglet combination used in Eq. (2), which is defined asf TMD ns =f TMD u f TMD d . BothB i and˜ S include logarithmic (⇠ ln a) and linear (⇠ 1/a) ultraviolet divergences, with the latter proportional to the total lengths of the Wilson lines. Both functions also include contributions diverging linearly as 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 6 K g L b K b A Z 6 V V v L v 7 n d R M T X v s p k 3 F i q C T L R W H C k Y n Q / H k 0 Y I o S w 6 e W Y K K Y v R W R E V a Y G B t R w Y b g r b 6 8 T l r V i u d W v L t a q V 7 N 4 s j D G Z x D G T y 4 g j r c Q g O a Q I D D M 7 z C m z N x X p x 3 5 2 P Z m n O y m V P 4 A + f z B 4 f e j 5 I = < / l a t e x i t > < 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 6 K g L b K b A Z 6 V V v L v 7 n d R M T X v s p k 3 F i q C T L R W H C k Y n Q / H k 0 Y I o S w 6 e W Y K K Y v R W R E V a Y G B t R w Y b g r b 6 8 T l r V i u d W v L t a q V 7 N 4 s j D G Z x D G T y 4 g j r c Q g O a Q I D D M 7 z C m z N x X p x 3 5 2 P Z m n O y m V P 4 A + f z B 4 f e j 5 I = < / l a t e x i t > < 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 6 K g L b K b A Z 6 V V v L v 7 n d R M T X v s p k 3 F i q C T L R W H C k Y n Q / H k 0 Y I o S w 6 e W Y K K Y v R W R E V a Y G B t R w Y b g r b 6 8 T l r V i u d W v L t a q V 7 N 4 s j D G Z x D G T y 4 g j r c Q g O a Q I D D M 7 z C m z N x X p x 3 5 2 P Z m n O y m V P 4 A + f z B 4 f e j 5 I = < / l a t e x i t > < 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 6 K g L b K b A Z 6 V V v L v 7 n d R M T X v s p k 3 F i q C T L R W H C k Y n Q / H k 0 Y I o S w 6 e W Y K K Y v R W R E V a Y G B t R w Y b g r b 6 8 T l r V i u d W v L t a q V 7 N 4 s j D G Z x D G T y 4 g j r c Q g O a Q I D D M 7 z C m z N x X p x 3 5 2 P Z m n O y m V P 4 A + f z B 4 f e j 5 I = < / l a t e x i t > q i (z µ + b µ ) < l a t e x i t s h a 1 _ b a s e 6 4 = " o T d L z m 9 5 P L r n g o Y n a 2 K E l O M Q I T k = " > A A A C B H i c b V B L S w M x G M z W V 6 2 v V Y + 9 L B a h I p R d E f R Y 9 O K x g n 1 A t y 7 Z N N u G 5 r E m W a E u P X j x r 3 j x o I h X f 4 Q 3 / 4 3 Z d g / a O p A w z H x D 8 k 0 Y U 6 K 0 6 3 5 b h a X l l d W 1 4 n p p Y 3 N r e 8 f e 3 W s p k U i E m 0 h Q I T s h V J g S j p u a a I o 7 s c S Q h R S 3 w 9 F l 5 r f v s V R E 8 B s 9 j n G P w Q E n E U F Q G y m w y 7 4 w d p Z O 7 y Y B q T 7 c + i w 5 D r P 7 K L A r b s 2 d w l k k X k 4 q I E c j s L / 8 v k A J w 1 w j C p X q e m 6 s e y m U m i C K J y U / U T i G a A Q H u G s o h w y r X j p d Y u I c G q X v R E K a w 7 U z V X 8 n U s i U G r P Q T D K o h 2 r e y 8 T / v G 6 i o / N e S n i c a M z R 7 K E o o Y 4 W T t a I 0 y c S I 0 3 H h k A k i f m r g 4 Z Q Q q R N b y V T g j e / 8 i J p n d Q 8 t + Z d n 1 b q F 3 k d R V A G B 6 A K P H A G 6 u A K N E A T I P A I n s E r e L O e r B f r 3 f q Y j R a s P L M P / s D 6 / A H m 0 5 g 6 < / 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 T d L z m 9 5 P L r n g o Y n a 2 K E l O M Q I T k = " > A A A C B H i c b V B L S w M x G M z W V 6 2 v V Y + 9 L B a h I p R d E f R Y 9 O K x g n 1 A t y 7 Z N N u G 5 r E m W a E u P X j x r 3 j x o I h X f 4 Q 3 / 4 3 Z d g / a O p A w z H x D 8 k 0 Y U 6 K 0 6 3 5 b h a X l l d W 1 4 n p p Y 3 N r e 8 f e 3 W s p k U i E m 0 h Q I T s h V J g S j p u a a I o 7 s c S Q h R S 3 w 9 F l 5 r f v s V R E 8 B s 9 j n G P w Q E n E U F Q G y m w y 7 4 w d p Z O 7 y Y B q T 7 c + i w 5 D r P 7 K L A r b s 2 d w l k k X k 4 q I E c j s L / 8 v k A J w 1 w j C p X q e m 6 s e y m U m i C K J y U / U T i G a A Q H u G s o h w y r X j p d Y u I c G q X v R E K a w 7 U z V X 8 n U s i U G r P Q T D K o h 2 r e y 8 T / v G 6 i o / N e S n i c a M z R 7 K E o o Y 4 W T t a I 0 y c S I 0 3 H h k A k i f m r g 4 Z Q Q q R N b y V T g j e / 8 i J p n d Q 8 t + Z d n 1 b q F 3 k d R V A G B 6 A K P H A G 6 u A K N E A T I P A I n s E r e L O e r B f r 3 f q Y j R a s P L M P / s D 6 / A H m 0 5 g 6 < / 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 T d L z m 9 5 P L r n g o Y n a 2 K E l O M Q I T k = " > A A A C B H i c b V B L S w M x G M z W V 6 2 v V Y + 9 L B a h I p R d E f R Y 9 O K x g n 1 A t y 7 Z N N u G 5 r E m W a E u P X j x r 3 j x o I h X f 4 Q 3 / 4 3 Z d g / a O p A w z H x D 8 k 0 Y U 6 K 0 6 3 5 b h a X l l d W 1 4 n p p Y 3 N r e 8 f e 3 W s p k U i E m 0 h Q I T s h V J g S j p u a a I o 7 s c S Q h R S 3 w 9 F l 5 r f v s V R E 8 B s 9 j n G P w Q E n E U F Q G y m w y 7 4 w d p Z O 7 y Y B q T 7 c + i w 5 D r P 7 K L A r b s 2 d w l k k X k 4 q I E c j s L / 8 v k A J w 1 w j C p X q e m 6 s e y m U m i C K J y U / U T i G a A Q H u G s o h w y r X j p d Y u I c G q X v R E K a w 7 U z V X 8 n U s i U G r P Q T D K o h 2 r e y 8 T / v G 6 i o / N e S n i c a M z R 7 K E o o Y 4 W T t a I 0 y c S I 0 3 H h k A k i f m r g 4 Z Q Q q R N b y V T g j e / 8 i J p n d Q 8 t + Z d n 1 b q F 3 k d R V A G B 6 A K P H A G 6 u A K N E A T I P A I n s E r e L O e r B f r 3 f q Y j R a s P L M P / s D 6 / A H m 0 5 g 6 < / 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 T d L z m 9 5 P L r n g o Y n a 2 K E l O M Q I T k = " > A A A C B H i c b V B L S w M x G M z W V 6 2 v V Y + 9 L B a h I p R d E f R Y 9 O K x g n 1 A t y 7 Z N N u G 5 r E m W a E u P X j x r 3 j x o I h X f 4 Q 3 / 4 3 Z d g / a O p A w z H x D 8 k 0 Y U 6 K 0 6 3 5 b h a X l l d W 1 4 n p p Y 3 N r e 8 f e 3 W s p k U i E m 0 h Q I T s h V J g S j p u a a I o 7 s c S Q h R S 3 w 9 F l 5 r f v s V R E 8 B s 9 j n G P w Q E n E U F Q G y m w y 7 4 w d p Z O 7 y Y B q T 7 c + i w 5 D r P 7 K L A r b s 2 d w l k k X k 4 q I E c j s L / 8 v k A J w 1 w j C p X q e m 6 s e y m U m i C K J y U / U T i G a A Q H u G s o h w y r X j p d Y u I c G q X v R E K a w 7 U z V X 8 n U s i U G r P Q T D K o h 2 r e y 8 T / v G 6 i o / N e S n i c a M z R 7 K E o o Y 4 W T t a I 0 y c S I 0 3 H h k A k i f m r g 4 Z Q Q q R N b y V T g j e / 8 i J p n d Q 8 t + Z d n 1 b q F 3 k d R V A G B 6 A K P H A G 6 u A K N E A T I P A I n s E r e L O e r B f r 3 f q Y j R a s P L M P / s D 6 / A H m 0 5 g 6 < / l a t e x i t > 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 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 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 > 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 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K L S P J b 3 Z p q g H 9 G R 5 C F n 1 F i p + T Q o V 9 y q u w B Z J 1 5 O K p C j M S h / 9 Y c x S y O U h g m q d c 9 z E + N n V B n O B M 5 K / V R j Q t m E j r B n q a Q R a j 9 b H D o j F 1 Y Z k j B W t q Q h C / X 3 R E Y j r a d R Y D s j a s Z 6 1 Z u L / 3 m 9 1 I Q 1 P + M y S Q 1 K t l w U p o K Y m M y / J k O u k B k x t Y Q y x e 2 t h I 2 p o s z Y b E o 2 B G / 1 5 X X S v q p 6 b t V r X l f q t T y O I p z B O V y C B z d Q h z t o Q A 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 > < 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 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K L S P J b 3 Z p q g H 9 G R 5 C F n 1 F i p + T Q o V 9 y q u w B Z J 1 5 O K p C j M S h / 9 Y c x S y O U h g m q d c 9 z E + N n V B n O B M 5 K / V R j Q t m E j r B n q a Q R a j 9 b H D o j F 1 Y Z k j B W t q Q h C / X 3 R E Y j r a d R Y D s j a s Z 6 1 Z u L / 3 m 9 1 I Q 1 P + M y S Q 1 K t l w U p o K Y m M y / J k O u k B k x t Y Q y x e 2 t h I 2 p o s z Y b E o 2 B G / 1 5 X X S v q p 6 b t V r X l f q t T y O I p z B O V y C B z d Q h z t o Q A 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 > < 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 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K L S P J b 3 Z p q g H 9 G R 5 C F n 1 F i p + T Q o V 9 y q u w B Z J 1 5 O K p C j M S h / 9 Y c x S y O U h g m q d c 9 z E + N n V B n O B M 5 K / V R j Q t m E j r B n q a Q R a j 9 b H D o j F 1 Y Z k j B W t q Q h C / X 3 R E Y j r a d R Y D s j a s Z 6 1 Z u L / 3 m 9 1 I Q 1 P + M y S Q 1 K t l w U p o K Y m M y / J k O u k B k x t Y Q y x e 2 t h I 2 p o s z Y b E o 2 B G / 1 5 X X S v q p 6 b t V r X l f q t T y O I p z B O V y C B z d Q h z t o Q A 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 > < 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 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K L S P J b 3 Z p q g H 9 G R 5 C F n 1 F i p + T Q o V 9 y q u w B Z J 1 5 O K p C j M S h / 9 Y c x S y O U h g m q d c 9 z E + N n V B n O B M 5 K / V R j Q t m E j r B n q a Q R a j 9 b H D o j F 1 Y Z k j B W t q Q h C / X 3 R E Y j r a d R Y D s j a s Z 6 1 Z u L / 3 m 9 1 I Q 1 P + M y S Q 1 K t l w U p o K Y m M y / J k O u k B k x t Y Q y x e 2 t h I 2 p o s z Y b E o 2 B G / 1 5 X X S v q p 6 b t V r X l f q t T y O I p z B O V y C B z d Q h z t o Q A 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 = " / 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 = " p 3 x A c p N + 1 A 5 k g o 0 d 9 t D Q B B K N D y M = " > A A
F e X c + 5 q 0 5 J 5 s 5 h D 9 w P n 8 A D I u Q f 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 = " p 3 x A c p N + 1 A 5 k g o 0 F e X c + 5 q 0 5 J 5 s 5 h D 9 w P n 8 A D I u Q f 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 = " p 3 x A c p N + 1 A 5 k g o 0 F e X c + 5 q 0 5 J 5 s 5 h D 9 w P n 8 A D I u Q f 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 = " p 3 x A c p N + 1 A 5 k g o 0 ⇠ ⌘/b T in the limit ⌘ ! 1 [33]. The ⌘/a and b T /a divergences, as well as renormalizes the remaining linear (⇠ b z /a) and logarithmic divergences in the quasi TMDPDF and matches it to the quasi TMDPDF with Dirac structure 4 (where '4' indexes the temporal direction) in the MS scheme at scale µ [34,38,39]. An alternate choice is to consider the quasi TMDPDF with Dirac structure 3 ; both 4 and 3 can be boosted onto + and thus define quasi TMDPDFs which can be matched to the spin-independent TMDPDF in the infinite-momentum limit.
Quasi beam function: The quasi beam functions in Eq. (3) are defined as matrix elements of quark bilinear operators with staple-shaped Wilson lines: Here h(P z ) denotes a boosted hadron state with four-momentum P µ = (0, 0, P z , E . It is convenient to define a dimensionless 'bare' nonsinglet beam function: The (4) is defined as a quark bilinear with a staple-shaped Wilson line, depicted in Fig. 1: where f W (⌘; b µ ; z µ ) is a spatial Wilson line of staple length ⌘ in theẽ z direction connecting endpoints separated by b µ = (b T , b z , 0). Here T denotes a direction transverse tõ e z , and all spatial Wilson lines are defined as Quasi soft factor: The quasi soft factor˜ S (b T , a, ⌘) can be computed as the vacuum matrix element of a closed spatial Wilson loop, whose definition and properties are detailed in Refs. [30][31][32][33]. This factor cancels in the ratios of quasi TMDPDFs which define the Collins-Soper evolution kernel by Eq. (2), and will thus not be discussed further here.
Renormalization factor: The renormalization factor Z MS 4 can be separated into two parts which renormalize the quasi beam function and soft factor respectively, denoted by Z MS O 4 and Z MS S : Both Z MS O 4 and Z MS S include linear power divergences proportional to ⌘/a and b T /a that cancel between the two terms, such that the complete renormalization factor Z MS 4 is independent of ⌘ and b T . Z MS O 4 can be computed nonperturbatively using the regularization independent momentum subtraction (RI 0 /MOM) scheme, with a perturbative matching to the MS scheme via a multiplicative factor R MS O 4 as described in Refs. [34,38]. In this approach, Z MS O 4 can be expressed as where Z is the RI 0 /MOM renormalization factor and p R denotes the matching scale introduced in the RI 0 /MOM scheme. At all orders in perturbation theory, the scheme conversion factor R MS O 4 cancels the depen- on p R and on the direction ofb T (up to discretization artefacts).
The authors have previously calculated Z , and thereby Z MS O 4 , by this approach in a quenched lattice QCD study [39]; those results are used for the numerical study in this work. The renormalization factor Z MS S does not need to be evaluated for a computation of the Collins-Soper kernel, as detailed in the following subsection.
Collins-Soper kernel: In the ratio of quasi TMDPDFs which gives the Collins-Soper kernel in Eq. (2),˜ S and its renormalization factor Z MS S , which do not depend on b z , cancel between the numerator and denominator. As a result, q ⇣ (µ, b T ) can be expressed in terms of the quasi beam function and its renormalization only, at the cost of introducing power-law divergences in ⌘ and b T separately in the numerator and denominator (divergences which were canceled by the quasi soft factor and its renormalization in the original expression for the kernel). Moreover, to ensure that the renormalization and matching between RI 0 /MOM and MS is performed in the perturbative region, the scale b T must be taken to be much smaller than ⇤ 1 QCD , a condition which does not permit an extraction of the Collins-Soper kernel at b T values in the nonperturbative region. A perturbative renormalization matching scale b T = b R T ⌧ ⇤ 1 QCD in Eq. (8) can, however, be defined by exploiting the b T -independence of Z MS 4 (µ, b z , a), as described in Ref. [34]. In this approach, for the choice = 4 in Eq. (3), the Collins-Soper kernel can be expressed as where a modified MS-renormalized quasi beam function B MS has been defined as Here, the additional factorR has been introduced into the modified MS-renormalized quasi beam function to compensate for the power-law divergences ⇠ |b T b R T |/a which would otherwise a↵ect both the numerator and denominator of Eq. (10): In this definition, fixed choices ofp R ,p 0 R , and of the directions ofb T andb R T , are taken. Since the factorR is independent of b z , and thus cancels between the numerator and denominator of Eq. (10), the specific choice of definition will not a↵ect the determination of the Collins-Soper kernel 2 . In the numerical study in this work, an average overb T andb R T orientations, and over several choices of p R andp 0 R , is performed in the same manner detailed in Appendix C in the numerator and denominator ofR. 3 2 The definition ofR used here di↵ers from that in Ref. [34] by the omission of the quasi soft factor and by allowingp 0 R to be di↵erent fromp R . 3 In the numerical study presented here, a set of ten momenta p R with p 2 R ranging from 5.7 to 28 GeV 2 , as described in Ref. [39], are used to constructR.   [41,42]. The lattice spacing a is determined from an analysis of scale setting in Ref. [43], and the lattice geometry parameters L and T are specified in units of a. For operator structures with Dirac index = 4 , n cfg configurations are analyzed, with nsrc source locations chosen on each. For other operator Dirac structures 6 = 4 , a subset with 25 configurations is analyzed, with 1 source location computed on each.
Several observations are pertinent to the computation of the Collins-Soper evolution kernel by Eq. (10). First, since the kernel is independent of the external state [32], one may calculate the quasi beam functions in the state with the best signal-to-noise properties in a lattice QCD calculation, e.g., for the pion. In a quenched calculation, a heavier-than-physical valence quark mass can be chosen for the same reason. Moreover, since although the kernel is state-independent, the power-corrections to the kernel are not, and so variation of the choice of external state, and external state momenta, provides a test of systematic e↵ects in a numerical calculation. Second, the Collins-Soper kernel does not depend on the longitudinal momentum fraction x or on the hadron momenta P z i , at O (b T /⌘, 1/(b T P z )). Although the truncation in the b zspace Fourier integral will induce oscillatory behavior in x-space, varying these parameters provides insight into these additional systematic uncertainties.
An alternative approach to extracting the Collins-Soper kernel by transforming the product of the matching coe cient and MS quasi beam function in Eq. (10) into a convolution integral in b z -space was advocated in Ref. [34]. Appendix E provides an investigation of this approach and finds that it su↵ers from significant systematic uncertainties.

III. LATTICE QCD STUDY
The Collins-Soper evolution kernel is computed by Eq. (10) in a lattice QCD calculation using a single quenched ensemble, detailed in Table I. The calculation is undertaken on gauge fields that have been subjected to Wilson flow to flow-time t = 1.0 [44], in order to increase the signal-to-noise ratio of the numerical results, and gauge-fixed to Landau gauge, in order to permit the use of gauge non-invariant quark wall sources. Quasi beam functions are constructed for a pion external state using valence quark propagators that are computed with the tree-level O(a) improved Wilson clover fermion action [45] and a  value that corresponds to a heavy pion mass of 1.207 (3) GeV. This choice may be made without introducing systematic bias, since the Collins-Soper kernel is independent of state. Three external state momenta are studied,P = P zẽ z with P z = n z 2⇡/L for  14) for pion states with momenta |P | = n z 2⇡/L. Shaded bands display the result of single-exponential fits to the two-point correlation functions for each non-zero momentum, and a twoexponential fit at zero momentum; the number of states in each fit is chosen to maximize an information criterion as described in the text, and the fit ranges shown correspond to the highest-weight fits in the weighted average over successful two-point function fits as discussed in Appendix A.
n z 2 {2, 3, 4}, corresponding to P z 2 {1.29, 1.94, 2.58} GeV, allowing the kernel to be computed from three di↵erent momentum ratios. To improve the overlap of boosted pion interpolating operators onto their respective ground states and improve statistical precision, a combination of wall sources and momentum-smeared sinks [46] are used to construct two-point and three-point correlation functions.
Bare quasi beam functions B bare (b z ,b T , a, ⌘, P z ) are extracted for non-local quark bilinear operators (Eq. (6)) with Wilson line staple geometries defined by staple extents ⌘ ranging between 0.6 and 0.8 fm (⌘/a 2 {10, 12, 14}), and with staple widths and asymmetries corresponding to |b T | and b z ranging from (⌘ a) to (⌘ a). In order for the mixing contributions to Eq. (10) to be consistently included, bare quasi beam functions are computed for all Dirac operator structures . As detailed in the caption of Table I, however, lower statistics are used for operators with Dirac structures 6 = 4 , whose contributions to the Collins-Soper kernel are suppressed by the renormalization factors. Previously, the 16-dimensional vector of MS renormalization factors Z MS O 4 0 (µ, b z ,b T , a, ⌘) was computed for the same ensemble and operator parameters as studied here [39], and those results are used in this work.
The two-point correlation function for the pion, projected to a given three-momentumP , is defined as: where ZP denotes the combination of overlap fac-tors for the source and sink interpolation operators and the ellipsis in Eq. (13) denotes contributions from higher excitations, which are exponentially suppressed for large t and discussed further in Appendix A. Wall-source interpolating operators ⇡P ,W (t) = u(t,P /2) 5 d(t,P /2) are used as sources for correlation functions, where momentum projected quark fields are defined by q(t,P ) = Px e iP ·x q(x, t) for q = {u, d}. Momentum-smeared interpolating operators ⇡P ,S (x, t) = u S(P /2) (x, t) 5 d S(P /2) (x, t) are used as sinks, where q S(P ) (x, t) are quasi local smeared quark fields obtained through iterative application of the Gaussian momentum-smearing operator defined in Ref. [46]. In particular, 50 steps of iterative momentum-smearing with smearing radius " = 0.25, as defined in Ref. [46], are used to construct momentum-smeared sinks for each momentum corresponding to n z 2 {2, 3, 4}. An e↵ective energy function that asymptotes to EP can be defined from the two-point correlation function by Two-point correlation functions for the three momenta which are considered here are displayed in Fig. 2. The extracted energies are slightly smaller than those obtained with the continuum dispersion relation, with relative deviations from EP = q m 2 ⇡ + |P | 2 ranging from 1.5(4)% for n z = 2 to 3.5(4)% for n z = 4. These deviations are consistent with the expected size of lattice artifacts which are neglected in this exploratory work.

A. Quasi beam functions
Bare quasi beam functions can be computed from three-point correlation functions with insertions of the non-local quark bilinear operators O i (b µ , z µ , ⌘), defined in Eq. (6). For the special case where pion momenta are taken only in the z-direction (i.e., consistent with the definition of quasi beam functions in Eq. (4)), three-point correlation functions are defined as A ratio of three-and two-point correlation functions then enables the bare isovector quasi beam functions of Eq. (5) to be extracted: Constraining the bare quasi beam functions B bare from ratios of two-and three-point functions R for all staple geometries (specified by {⌘, b µ }), all Dirac structures , and all momenta P z , considered in this work, requires fits for a very large number of operators (35,660) to be performed. These fits are automated using a fit procedure discussed in Appendix A. An example of the result of these fits, for = 4 , and specific choices of b T and ⌘, is given in Fig. 3(a); a second example figure holding b z fixed, but showing all Dirac structures, is shown in Fig. 3(b). Additional examples of the real and imaginary parts of the extracted bare quasi beam functions are displayed in Appendix B.
The bare quasi beam functions obtained by Eq. (16) are renormalized to the MS scheme by Eq. (11), using renormalization factors Z MS O 4 which were computed for the same ensemble and operators as studied here in Ref. [39]. The fractional contributions to the renormalized quasi beam function from the bare quasi beam functions with di↵erent Dirac structures is shown in Fig. 3(c), for a particular choice of parameters. The size of these contributions is observed to grow with increasing b T and with increasing (⌘ b z ); while the relative magnitudes of bare beam functions with di↵erent do not vary significantly with these parameters, the relative importance of the o↵-diagonal renormalization factors varies significantly as discussed in Ref. [39]. Across the parameters studied, the combined contributions from mixing to the renormalized quasi beam function with Dirac structure = 4 are at the 5%-25% level. Calculations of the quasi beam functions, and the Collins-Soper evolution kernel, to better than this precision thus require bare quasi beam functions to be computed for several Dirac structures .
The functional dependence of the renormalized quasi beam function B MS 4 (µ, b z ,b T , a, ⌘, b R T , P z ) (defined in Eq. (11)) is shown in Fig. 4. The factorR(b T , b R T , a, ⌘) (Eq. (12)) was included in the definition of the renormalized quasi beam function to cancel the dependence of the bare beam function on ⌘ and on b R T . It is clear that over choices of b R T within the perturbative region, this dependence is indeed removed to better than the statistical uncertainties of this study. A weighted average of the renormalized quasi beam function over these parameters, as well as over di↵erent directions ofb T , is thus taken as detailed in Appendix C to define averaged quasi beam functions B shown in Appendix C. These are the key results used to extract the Collins-Soper kernel, as discussed in the next subsection.

B. Collins-Soper kernel
Computing the Collins-Soper evolution kernel by Eq. (10) requires taking the Fourier transform of the MSrenormalized quasi beam function with respect to b z . It is clear from the results shown in Fig. 5, however, that with the parameter ranges explored in this study the Fourier transform will su↵er from significant truncation e↵ects since the quasi beam function is not yet consistent with zero within uncertainties at the largest b z values that are used. For this reason, models are used to fit the P z b z -dependence of the lattice data for the quasi beam function before the Fourier transform is taken to evaluate the Collins-Soper kernel. The results obtained by taking discrete Fourier transforms instead are shown in Appendix D, and a discussion of what will be required for future calculations to achieve robust results for the Collins-Soper kernel without this modeling step is presented in Sec. IV.
Two models of the P z b z -dependence of the quasi beam functions are considered, based on Hermite and Bernstein polynomial bases. The models are constructed to yield xindependent Collins-Soper kernels, as would be expected in the absence of systematic artifacts, assuming the leading order value for the perturbative matching coe cient, i.e., C TMD ns = 1. Including higher orders in the matching factor while guaranteeing an x-independent kernel would require more complicated functional forms to be fit to the quasi beam functions to compensate for the x-dependence of the matching. It is expected that the matching uncertainties are small relative to the systematic uncertainties inherent in introducing models for the P z b z -dependence of the quasi beam functions, and these e↵ects are thus neglected in this work. While a modelindependent result for the Collins-Soper kernel cannot be achieved from the data presented here, the comparison between results obtained using the two di↵erent models considered nevertheless provides some indication of the severity of the model-dependence, and the quality of fits to these functional forms not including power corrections also provides a measure of their importance.
The first functional form which is fit to the MSrenormalized quasi beam function is where H n (x) is the n-th Hermite polynomial. The fit parameter ! is taken to be complex, while the other free parameters are real. Allowing Im(!) 6 = 0 allows the T and ⌘ to the renormalized quasi beam function as a function of b z and P z (at the fixed a of the calculation), as described in the text. (a) Fourier transform of F Herm N (P z , b z P z ; {a k }, , !, ) with respect to b z P z to be complex, and correspondingly enables F Herm N (P z , b z P z ; {a k }, , !, ) to be an asymmetric function of b z P z . The real and imaginary parts of the quasi beam function are symmetric and antisymmetric functions of b z respectively in the ⌘ ! 1 limit; however, the numerical results presented in this work show significant departures from these expectations, particularly for large b T , as shown in Fig. 5(b). The observed asymmetry could arise from finite-volume e↵ects: e↵ective field theory calculations [49] have demonstrated that finite-volume e↵ects for pion matrix elements of nonlocal operators with separation`generically take the form e m⇡(L `) . In this work, one therefore expects b zdependent finite-volume e↵ects of the form e m⇡(L ⌘+b z ) as well as additional b z independent finite-volume effects. In addition, exponential dependence on b z could arise from an imperfect cancellation between power-lawdivergent lattice artifacts in B bare (b z ,b T , a, ⌘, P z ) and . Taking Im(!) 6 = 0 al-lows the fit form in Eq. (17) to include exponential dependence on b z and is found to significant improve the quality of fits to the numerical results with large b T & 0.5 fm.
The second model considered assumes that the Fourier transform of the quasi beam function has compact support on the interval 0 < x < 1 [30,31,33], which is expected to become valid for large P z , and takes the form where B r,N 1 , for r 2 {0, . . . N 1} are the N Bernstein basis polynomials of degree N 1 normalized as in Ref. [50], and asymmetry in b z is accommodated by taking Im(a r ) 6 = 0.
Using either functional form, F Herm  (18). That is, the resulting Collins-Soper kernel is independent of x by construction. The full procedure by which each functional form is fit to the numerical results for the quasi beam function is described in Appendix C, and examples of the resulting fits are shown both in Fig. 6 and in Appendix C. Briefly, the fits are undertaken simultaneously at all P z and b z values for a given b T , and an information criterion is used to choose the model truncation N for each fit. While both models fit the quasi beam function well within the range of P z b z values constrained by the lattice data (with an average 2 /N dof over all fits of 0.9, tabulated in Appendix C), it is clear from Fig. 6 that they correspond to substantially di↵erent models outside this range.
The Collins-Soper kernel determined from each set of model fits is shown in Fig. 7. The results obtained using the two model forms, i.e., the Hermite polynomial model, in which the quasi beam function has support on 1 < x < 1, and the Bernstein polynomial model, with support on 0 < x < 1, are consistent. This encouragingly suggests that q,MS ⇣ is well-constrained by the numerical results at the P z and b z values of this calculation, and that the model-dependence introduced in the Fourier transform is relatively mild. Perturbative results for the 0-flavor Collins-Soper kernel [47,48] are  also shown in Fig. 7 for comparison. 4 It is noteworthy that the lattice QCD results for q,MS ⇣ obtained here are consistent with perturbative calculations of the 0-flavor Collins-Soper kernel [47,48] in the region b T ⇠ 0.2 fm. For b T < 0.2 fm, the results di↵er significantly from the perturbative calculation, which is likely due to power corrections to the lattice QCD results of the form 1/(b T P z ), which have not been estimated here.
Although the Collins-Soper kernel shown in Fig. 7 has been obtained in a 0-flavor calculation, it can also be compared qualitatively with the results of fits to experimental data, in which several di↵erent parametrizations of the nonperturbative behavior of the kernel at large b T have been used. In early fits to Drell-Yan data [4][5][6], the Collins-Soper kernel was parametrized as a quadratic function in b T in the nonperturbative region. It was later found in Refs. [53,54], however, that these fits cannot describe SIDIS data. More recently, it has been argued that q,MS ⇣ should approach a constant as b T ! 1; phenomenological fits to Drell-Yan data under this assumption suggest that this constant is ⇠ 0.6 [55]. Finally, in a recent fit to both SIDIS and Drell-Yan data [14], the kernel was parametrized to behave linearly at large b T . The results of this numerical study are qualitatively consistent with constant or linear behavior of the kernel in the nonperturbative region. Once lattice QCD results with controlled systematic uncertainties are available, it will be possible to test these and other phenomenological expectations for the large-b T behavior of the Collins-Soper kernel with QCD predictions, and begin incorporating lattice QCD constraints in phenomenological analyses. The requirements for a fully controlled lattice QCD determination of the Collins-Soper kernel are discussed in the following section.

IV. OUTLOOK
This manuscript presents an exploratory calculation of the nonperturbative Collins-Soper kernel in quenched lattice QCD based on the method developed in Refs. [32][33][34]. In this approach, the kernel is computed from quasi beam functions defined from matrix elements of quark bilinear operators with staple-shaped Wilson lines in boosted hadron states. These beam functions are 4 In this work, ↵ [51] to lower scales using the four-loop function [52], integrating out bottom and charm quarks, and finally matching ↵ is calculated in Ref. [39]. This procedure gives the result ⇤ renormalized to the MS scheme via the RI 0 /MOM prescription, and a ratio of Fourier-transformed quasi beam functions at di↵erent hadron boost momenta determines the Collins-Soper kernel. In this calculation, the kernel is extracted over a range of scales b T 2 (0.1, 0.8) fm. The final results presented here rely on modeling the b z -space quasi beam functions to control truncation e↵ects in the Fourier transform. Nevertheless, the results are robust under the variation of models considered here.
For a future controlled and model-independent determination of the Collins-Soper kernel by this method, several improvements will be essential. Critically, larger lattice volumes must be studied; the overwhelming systematic in this calculation arises from modeling to facilitate the Fourier transform, which is necessary because of truncation e↵ects su↵ered due to the limited b z range over which quasi beam functions could be computed on the lattice volume used here. One measure of truncation effects is given by the model truncation error defined as where the truncated discrete Fourier transform (DFT) and untruncated Fourier transform are defined as Equation (19) is evaluated by applying both the DFT and untruncated Fourier transform to the best-fit model for large b T and prevent a reliable determination of the Collins-Soper kernel using a DFT approach. Results with significantly larger P z b z max , however, could be used to obtain a model independent prediction of the Collins-Soper kernel directly from a DFT of lattice QCD results (see Appendix D). For example, based on the results of the present study, quasi beam function calculations with (b z max , P z ) ⇠ (2.5 fm, 2.5 GeV) are likely to su↵er from truncation e↵ects of less than 5% for b T scales across the range of those studied in this work.
In addition to truncation artifacts, extractions of the Collins-Soper kernel by the method pursued here su↵er from power corrections at O (b T /⌘, 1/(b T P z )) [33,36]. These e↵ects could not be resolved by model fits in this study, and as such, the coe cients of these power corrections could not be constrained. Nevertheless, larger physical lattice volumes, as needed to reduce truncation artifacts, will simultaneously enable O (b T /⌘, 1/(P z ⌘)) e↵ects to be reduced by allowing larger staple extents ⌘ to be investigated at fixed b T . Larger boost momentum, again needed to control truncation e↵ects, will also simultaneously enable control over power corrections of O(1/(b T P z )). These artifacts make comparison of the Collins-Soper kernel extracted by the method pursued here with perturbative predictions, which are accurate in the region b T ⌧ ⇤ 1 QCD , challenging. Precise comparisons in this region will be an important test of systematics in the lattice QCD approach. The infinite volume and continuum limits must also ultimately be taken for a fully controlled result.
Future studies would also gain significantly by exploiting the state-independence of the Collins-Soper kernel to obtain multiple constraints on the kernel from the same calculation and thus test systematic e↵ects. An alternative, complementary, approach to extracting the Collins-Soper kernel from Eq. (2) was proposed in Ref. [37]. This strategy uses the Mellin moments of the expressions so that one only needs to calculate the quasi beam function or its derivatives at b z = 0, which reduces the computational cost and has the advantage that renormalization factors cancel in the ratio. This approach also, however, requires a nontrivial integration over the TMDPDF that is extracted from experiments over a limited kinematic range. Comparison of results of the two approaches will also be valuable in future calculations. Despite the significant challenges described above, the results presented here suggest that controlled firstprinciples calculations of the Collins-Soper kernel at nonperturbative scales as large as b T ⇠ 1 fm are tractable with current methods. Refs. [11,[13][14][15] indicate that such a calculation at 10% precision at scales b T 2 (0.2, 1.0) fm will be su cient to di↵erentiate di↵erent models of the Collins-Soper kernel and will thereby provide important input for fitting low-energy SIDIS data. Ultimately, larger values of b T , e.g. b T . 2 fm, will be important input for determinations of the TMDPDFs; this will be attainable with larger lattice volumes in the future.
viding the gauge field configurations used in this project. Calculations were performed using the Qlua [56] and Chroma [57]  Appendix A: Fits to three-and two-point functions As shown in Eq. (16), ratios of three-point and twopoint correlation functions asymptote to the bare quasi beam function in the double limit {⌧, t ⌧ } ! 1. Ratios computed at finite t and ⌧ , however, will include contributions from matrix elements in excited states. Twopoint correlation functions have the spectral representation where n labels QCD energy eigenstates with the quantum numbers of a pion with momentumP , E ñ P is the energy of state n, q Z ñ P = 2E ñ P h0| ⇡P |ni are overlap factors for the interpolating operator ⇡P onto state n, and thermal e↵ects arising from the finite Euclidean time extent T are included. Dependence of Z ñ P on the type of source/sink interpolating operators used (wall or momentum-smeared) is suppressed in Eq. (A1) and throughout this section. Fits of lattice QCD two-point correlation function results to Eq. (A1) can be used to extract the ground-state energies E 0 P shown in Fig. 2, as well as excited-state energies and overlap factors.
Three-point functions have an analogous spectral representation where n, m index energy eigenstates. Combing Eq. (A1) and Eq. (A2) and isolating the ground-state contributions yields a spectral representation for ratios of three-and two-point functions: After determining the spectral quantities appearing on the left-hand-side of Eq. (A3) from fits to lattice QCD results for C 2pt , where in practice the sum over states is truncated at n = N states as discussed below, the bare quasi beam functions and other parameters appearing on the right-hand-side of Eq. (A3) can be determined from fits to lattice QCD results for three-point to twopoint function ratios. Fitting directly to these ratios has the advantages that ground-state overlap factors cancel exactly between three-and two-point functions and that correlated ratios are determined more precisely than three-point functions alone. Including the additional factor on the left-hand-side of Eq. (A3), which depends only on energies and overlaps obtained in two-point function fits, removes the need to model excited-state contamination in the two-point function during fits to the ratio (which would require fitting several additional parameters entering 2 -minimization nonlinearly) without spoiling these correlations. Three-point correlation functions are computed for six source/sink separations t/a 2 {9, 12, 15, 18, 21, 24} and all operator insertion points 0 < ⌧ < t. Signal-to-noise ratios of two-point and three-point correlation functions are proportional to e (E 0 P m⇡)t , where m ⇡ is the pion mass, and for n z 3 the large-separation results are very noisy and so only results with t/a 2 {9, 12, 15, 18} are used in fits. Correlated 2 -minimization fits of twopoint functions to Eq. (A1), followed subsequently by fits to Eq. (A3), are performed in an automated manner as follows: • The minimum source/sink separation t min are varied over the range 2  t min  t max t plateau , where t max is chosen to be the largest t for which the signal-to-noise ratio of C 2pt (t, a,P ) is greater than a fixed value (a threshold of 2 is used in the results presented here) and t plateau is a free parameter specifying the minimum number of points in a fit (results presented here use t plateau = 3). The re-striction t min 2 is set by the degree of nonlocality in the lattice action. For every possible choice of t min within this range, correlated 2 -minimization fits to Eq. (A1) are performed using two-point function results with t min  t  t max . The two-point function fitting procedure is identical to that described in Appendix B of Ref. [58]. To summarize, one-state fits are performed first, followed by two-state fits. If the Akaike information criterion (AIC) [59] is improved su ciently by the addition of a second state to the fit (a threshold of AIC < 2N dof , where N dof is the number of degrees of freedom in the fit, is used in the final results), then a three-state fit is performed and the same criterion is used to judge whether the three-state fit is preferred. This procedure is repeated until adding additional states does not su ciently improve the fit, in order to select the optimal number of states to include in the fit for each t min . The best-fit parameters are determined using nonlinear optimization for the energies E ñ P , with linear systems of equations solved to determine Z ñ P at each iteration. Covariance matrices are determined using optimal shrinkage [60] as described in Refs. [58,61] in order to reduce finite-statistic e↵ects leading to poorly conditioned sample covariance matrices. Several checks on numerical 2 optimization described in Ref. [58] are then performed to verify the reliability of the fit. If these checks are passed, an acceptable two-point function model has been found for this choice of t min , and three-to two-point function ratios are subsequently analyzed using fits to Eq. (A3).
• The minimum source/operator and source/sink separations corresponding to ⌧ 2 [⌧ min , ⌧ max ] are varied over the ranges 2  ⌧ min  (t min t plateau )/2 and 2  t ⌧ max  (t min t plateau )/2. Threepoint to two-point ratios using all available t 2    Examples of fits to the ratio of three-and two-point functions R (t, ⌧, b µ , a, ⌘, P z ) (Eq. (16)), obtained as described in the text. Shaded bands of the same colors as the points show 68% bootstrap confidence intervals of the ⌧ and t-dependent fits from the fit range (specifically the choice of tmin, ⌧min, and ⌧max) that had the highest weight in the weighted average of successful fits. Gray horizontal bands show the total uncertainty on the bare quasi beam functions extracted from the fits, including the statistical uncertainty and the systematic uncertainty from variation of the results between di↵erent fit range choices.
a, ⌘ a], and all P z corresponding to n z 2 {2, 3, 4}, are performed using uncorrelated 2 -minimization 5 In order to determine the appropriate polynomial to use in Eqs. (17)(18), the AIC is employed. Fits to Eq. (17) with Hermite polynomials of degree N = {0, 1, 2} are performed, and minimum-AIC fits with b T /a 2 [1,13] are obtained with minimum 2 /N dof as given in Table II. Fits to Eq. (18) with Bernstein polynomials of degree N = {0, 1, 2, 3} are similarly performed, and minimum-AIC fits are obtained with minimum 2 /N dof as given in Table III 1  2  3  4  5  6  7  8  9  10  11  12  13  N  1  3  3  2  0  3  2  2  2  1  3  1    The instability in Fig. 12(b) can be understood as a consequence of the limited range of b z and ⌘ considered in this study. From Eq. (10), the ratio of quasi beam functions should stabilize in a certain x-region for a given momentum pair {P z 1 , P z 2 }. Given a limited range of b z in the Fourier transform, this stabilization can be expected to be robust for values of x ⇠ 0.5 only. One might naively choose the peak at around x ⇠ 0.5 as the central value, and take some variation around the peak to define the systematic uncertainties. However, since |b z |  (⌘ a) in this calculation, the Fourier transform will introduce an oscillatory term to the quasi beam function with frequency P z (⌘ a). The same issue has been encountered in lattice QCD calculations of collinear PDFs [67][68][69][70], and certain model assumptions have been suggested to avoid this challenge [71,72]. For P z = 6⇡/L, ⌘ = 12a, as in this calculation, the period of this oscillation is T = 2⇡/(P z (⌘ a)) ⇠ 1.0, and the oscillatory behavior is not apparent within the region 0 < x < 1, as shown in Fig. 12(a). The oscillatory behavior will persist in ratios of quasi beam functions at di↵erent momenta P z 1 , P z 2 , as an interference between oscillations with frequen-cies P z 1 (⌘ a) and P z 2 (⌘ a). As a result, the Collins-Soper kernel extracted from the peak can be shifted significantly, which adds an important systematic error to numerical calculations via this approach.
In future calculations with increased ranges of ⌘ or P z , such that there are more rapid oscillations of the DFTs of quasi beam functions within the range 0 < x < 1, this approach may nevertheless be used robustly. For example, if the frequency P z (⌘ a) were doubled, then the Collins-Soper kernel would oscillate around the true value for at least two cycles, which would allow it to be determined by taking an average of the central local maximum and minimum within the oscillating region. To illustrate this point, a toy model for the quasi beam functionB ns in x-space is constructed: To study the oscillatory behavior in this toy model, the inverse FT of the quasi beam functionB ns (x, b T , µ, P z ) is taken first. Then, a DFT of the truncated quasi beam functions back to x-space is performed for ⌘ = {12a, 24a, 36a}, and the Collins-Soper kernel is computed using Eq. (10). The results are shown in Fig. 13. It is apparent that the kernel su↵ers from oscillations due to the truncated DFT, and for ⌘ = 12a the shape of the curve is qualitatively similar to those in Fig. 12(b). Moreover, the peaks or local maximums around x = 0.5 are significantly shifted from the true value of the Collins-Soper kernel for all ⌘ choices. Nevertheless, for ⌘ = 24a and 36a, taking the average of the central peak and trough values provides a close approximation to the Collins-Soper kernel, as shown in Fig. 14. With more rapid oscillations, this averaging method will lead to more accurate results. Future calculations with larger lattices sizes and higher pion momenta will thus likely enable reliable determination of the Collins-Soper kernel with the DFT method, although the toy model results shown in Fig. 14  large ⌘ values may be required to achieve percent-level precision.
Appendix E: Alternate approach in position space In this appendix an alternate approach to extract the Collins-Soper kernel in b z -space is investigated, as suggested in Ref. [34]. By taking the FT of the matching kernel C ns (µ, xP z ), the Collins-Soper kernel can be expressed as q,FI ⇣ (µ, b T ) = 1 ln(P z 1 /P z 2 ) ⇥ ln R db zC ns (µ, y b z P z 1 , P z 1 )P z 1 B MS 4 (b z ,b T , µ, P z 1 ) R db zC ns (µ, y b z P z 2 , P z 2 )P z 2 B MS 4 (b z ,b T , µ, P ⇥ ln R db zC 0 ns (µ, y b z P z 1 , P z 2 )P z 1 B MS 4 (b z ,b T , µ, P z 1 ) R db zC 0 ns (µ, y b z P z 2 , P z 1 )P z 2 B MS 4 (b z ,b T , µ, P z 2 ) , wherē C 0 ns (µ, b z P z , P z ) ⌘ Z dx e ix(b z P z ) C ns (µ, xP z ) . (E4) Eqs. (E1) and (E3) are denoted as Form I and Form II, respectively; studying both forms enables a consistency check. Similar to the results obtained via the DFT method outlined in App. D, the Collins-Soper kernel obtained by either Form I or Form II should not depend on the value of y or on the momentum pair {P z 1 , P z 2 }, which provides another handle on the relevant systematic uncertainties.
Collins-Soper kernels extracted with the position-space approach for the toy model of Eq. (D1) are shown in Fig. 15. The two forms do not yield consistent answers, which indicates that they are not numerically equivalent and that the perturbative convergence is lost in either or both of the convolution integrals. Nevertheless, it is clear that with Form I the extracted Collins-Soper kernel does not stabilize to the correct result with increasing ⌘, while with Form II the ratio stabilizes around the true value for su ciently large ⌘. With su ciently large ⌘ it is possible that this approach may provide a reliable determination of the Collins-Soper kernel, although this will need to be investigated carefully in future work.
The results of applying the position space approach to the lattice QCD results in this study are shown in Fig. 16, which is compared to the Collins-Soper kernel extracted from the fits with the Bernstein polynomial model, dis-cussed in the main text. With the range of ⌘ values computed in the numerical study there is no plateau in the y-space Collins-Soper kernels, and the di↵erent choices of momentum pairs do not appear to converge. Although the result extracted at {P z 1 , P z 2 } = {3, 4} ⇥ 2⇡/L is consistent with the results extracted using the model fits at the minimum value, this consistency is not found at di↵erent b T values. Further investigation is needed to confirm whether the position-space approach via Form I or II can provide a valuable consistency check against the DFT approach with larger physical lattice volumes used in calculations.