Collins-Soper Kernel for TMD Evolution from Lattice QCD
FFERMILAB-PUB-20-102-TMIT/CTP-5184
Collins-Soper Kernel for TMD Evolution from Lattice QCD
Phiala Shanahan, ∗ Michael Wagman,
1, 2, † and Yong Zhao
1, 3, ‡ Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA, USA 02139 Fermi National Accelerator Laboratory, Batavia, IL 60510, USA Physics Department, Brookhaven National Laboratory, Bldg. 510A, Upton, NY 11973, USA
The Collins-Soper kernel relates transverse momentum-dependent parton distribution functions(TMDPDFs) at different energy scales. For small parton transverse momentum q T ∼ Λ QCD , thiskernel is non-perturbative and can only be determined with controlled uncertainties through ex-periment or first-principles calculations. This work presents the first exploratory determination ofthe Collins-Soper kernel using the lattice formulation of Quantum Chromodynamics. In a quenchedcalculation, the N f = 0 kernel is determined at scales in the range 250 MeV < q T < I. INTRODUCTION
Understanding the structure of matter has been adefining goal of physics for centuries. In the mod-ern context, a primary objective is imaging the three-dimensional spatial and momentum structure of the pro-ton, and of other hadrons. Some important aspectsof this structure related to the transverse momentumof quarks and gluons in a hadron state are encodedin transverse-momentum-dependent parton distributionfunctions (TMDPDFs) [1–3]. Experimentally, thesequantities can be constrained for the proton by Drell-Yan production and semi-inclusive deep inelastic scat-tering (SIDIS) of electrons off protons; the best currentconstraints are achieved via global fits to experimentaldata [4–15], with improvements expected in the comingyears from measurements at COMPASS [16], the ThomasJefferson National Accelerator Facility [17], RHIC [18],and an Electron-Ion Collider [19]. Additional experi-mental information from dihadron production in e + e − collisions at Belle, and new ways of looking at hadronsinside jets [20, 21], may also help constrain these fits.Key to global fits of TMDPDFs is the ability to relatethese distributions determined in different processes, in-cluding those at different scales. That is, for a TMDPDF f TMD i (cid:0) x,(cid:126)b T , µ , ζ (cid:1) , defined for a parton of flavor i withlongitudinal momentum fraction x , transverse displace-ment (cid:126)b T (the Fourier conjugate of the transverse mo-mentum (cid:126)q T ), virtuality scale µ , and hadron momentumscale ζ which is related to the hard scale of the scatter-ing process, it is critical to understand its evolution toscales ( µ, ζ ): f TMD i (cid:16) x,(cid:126)b T , µ, ζ (cid:17) = f TMD i (cid:16) x,(cid:126)b T , µ , ζ (cid:17) × exp (cid:20)(cid:90) µµ d µ (cid:48) µ (cid:48) γ iµ ( µ (cid:48) , ζ ) (cid:21) exp (cid:20) γ iζ ( µ, b T ) ln ζζ (cid:21) , (1) ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] where b T = | (cid:126)b T | . The first exponential in this equationgoverns the µ -evolution of the TMDPDF, which is per-turbative for scales { µ , µ } (cid:29) Λ QCD . The evolution in ζ governed by the second exponential, however, is en-coded in the Collins-Soper kernel γ iζ ( µ, b T ), which is in-herently non-perturbative for q T ∼ b − T ∼ Λ QCD , even for µ (cid:29) Λ QCD . Experimentally, the Collins-Soper kernel canbe extracted by simultaneous global fits with the TMD-PDF, and recent global analyses show some discrepancyin determinations of the kernel in the region q T ≤ q T scales, including in the non-perturbativeregion.Section II outlines the procedure, developed inRefs. [32, 33], for constraining the Collins-Soper kernel The Collins-Soper kernel is also often denoted by K ( b T , µ ) [3],and it is defined as −D ( b T , µ ) in Ref. [12]. a r X i v : . [ h e p - l a t ] A ug using lattice QCD and LaMET. Section III details thequenched lattice QCD calculation undertaken here, in-cluding discussion of the systematic uncertainties in thecalculation, while Sec. IV outlines the requirements for afully-controlled calculation of the Collins-Soper kernel tobe achieved by this method. II. COLLINS-SOPER KERNEL FROM LATTICEQCD
In Refs. [32, 33] a method was proposed to determinethe quark Collins-Soper kernel using lattice QCD andLaMET. Precisely, it was shown that γ qζ ( µ, b T ) can be ex-tracted from a ratio of nonsinglet quasi TMDPDFs ˜ f TMDns at different momenta, which are defined using equal-timecorrelation functions within hadron states at large mo-mentum in the z -direction: γ qζ ( µ, b T ) = 1ln( P z /P z ) × ln C TMDns ( µ, xP z ) ˜ f TMDns ( x,(cid:126)b T , µ, P z ) C TMDns ( µ, xP z ) ˜ f TMDns ( x,(cid:126)b T , µ, P z ) , (2)up to power corrections which are discussed further be-low. In this expression, P zi (cid:29) Λ QCD are the z -componentof the hadron momenta and C TMDns is a perturbativematching coefficient that has been obtained at one-looporder [32, 33]. The quasi TMDPDF ˜ f TMDns , defined below,approximates the physical TMDPDF involving light-likepaths, as detailed in Ref. [33], and complications involv-ing matching in the soft sector [31, 33, 35, 36] are elim-inated in the ratio that gives the Collins-Soper kernel.Similar constructions have been used in calculations ofratios of x -moments of TMDPDFs from lattice QCD [23–27].The unpolarized quasi TMDPDF is defined in terms ofa quasi beam function ˜ B Γ i and a quasi soft factor ˜∆ S [30–33], both of which are calculable in lattice QCD:˜ f TMD i (cid:0) x,(cid:126)b T , µ, P z (cid:1) ≡ lim a → η →∞ (cid:90) d b z π e − i b z ( xP z ) Z MS γ Γ ( µ, b z , a ) × P z E (cid:126)P ˜ B Γ i (cid:0) b z ,(cid:126)b T , a, η, P z (cid:1) ˜∆ S ( b T , a, η ) , (3)where a denotes the lattice spacing, the subscript i is theflavor index, and summation over Dirac structures is im-plied. This summation accounts for the operator mixingsamong different Dirac structures in lattice QCD calcula-tions defined on a hypercubic space-time lattice [38–40].Additional mixing with gluon operators, not shown inEq. (3), cancels in the flavor nonsinglet combination usedin Eq. (2), which is defined as ˜ f TMDns = ˜ f TMD u − ˜ f TMD d .Both ˜ B Γ i and ˜∆ S include logarithmic ( ∼ ln a ) and lin-ear ( ∼ /a ) ultraviolet divergences, with the latter pro-portional to the total lengths of the Wilson lines. Bothfunctions also include contributions diverging linearly as T
The quasi beam functions inEq. (3) are defined as matrix elements of quark bilinearoperators with staple-shaped Wilson lines:˜ B Γ i ( b z ,(cid:126)b T , a, η, P z ) = (cid:68) h ( P z ) (cid:12)(cid:12) O i Γ ( b µ , , η ) (cid:12)(cid:12) h ( P z ) (cid:69) . (4)Here h ( P z ) denotes a boosted hadron state withfour-momentum P µ = (0 , , P z , E ( h ) (cid:126)P ), with E ( h ) (cid:126)P = (cid:113) (cid:126)P + m h and where m h is the mass of thehadron h . States are normalized as (cid:104) h ( (cid:126)P (cid:48) ) | h ( (cid:126)P ) (cid:105) =2 E ( h ) (cid:126)P (2 π ) δ (3) ( (cid:126)P − (cid:126)P (cid:48) ). It is convenient to define a di-mensionless ‘bare’ nonsinglet beam function: B bareΓ ( b z ,(cid:126)b T , a, η, P z ) = 12 E (cid:126)P (cid:16) ˜ B Γ u ( b z ,(cid:126)b T , a, η, P z ) − ˜ B Γ d ( b z ,(cid:126)b T , a, η, P z ) (cid:17) . (5)The operator O i Γ ( b µ , , η ) in Eq. (4) is defined as a quarkbilinear with a staple-shaped Wilson line, depicted inFig. 1: O i Γ ( b µ , z µ , η ) = ¯ q i ( z µ + b µ ) Γ2 W ˆ z ( z µ + b µ ; η − b z ) × W † T ( z µ + η ˆ z ; b T ) W † ˆ z ( z µ ; η ) q i ( z µ ) ≡ ¯ q i ( z µ + b µ ) Γ2 (cid:102) W ( η ; b µ ; z µ ) q i ( z µ ) , (6)where (cid:102) W ( η ; b µ ; z µ ) is a spatial Wilson line of staple length η in the (cid:126)e z direction connecting endpoints separated by b µ = ( (cid:126)b T , b z , T denotes a direction transverse to (cid:126)e z , and all spatial Wilson lines are defined as W ˆ α ( x µ ; η ) = P exp (cid:20) i g (cid:90) η d s A α ( x µ + s ˆ α ) (cid:21) . (7) Quasi soft factor:
The quasi soft factor ˜∆ S ( b T , a, η )can be computed as the vacuum matrix element of aclosed spatial Wilson loop, whose definition and prop-erties are detailed in Refs. [30–33]. This factor cancels inthe ratios of quasi TMDPDFs which define the Collins-Soper evolution kernel by Eq. (2), and will thus not bediscussed further here. Renormalization factor:
The renormalization factor Z MS γ Γ can be separated into two parts which renormal-ize the quasi beam function and soft factor respectively,denoted by Z MS O γ and Z MS S : Z MS γ Γ ( µ, b z , a ) = Z MS O γ ( µ, b z , b T , a, η ) Z MS S ( µ, b T , a, η ) . (8)Both Z MS O γ and Z MS S include linear power divergencesproportional to η/a and b T /a that cancel between thetwo terms, such that the complete renormalization factor Z MS γ Γ is independent of η and b T . Z MS O γ can be computednonperturbatively using the regularization independentmomentum subtraction (RI (cid:48) / MOM) scheme, with a per-turbative matching to the MS scheme via a multiplicativefactor R MS O γ as described in Refs. [34, 38]. In this ap-proach, Z MS O γ can be expressed as Z MS O γ ( µ, b z , b T , a, η ) = R MS O γ ( µ, p R , b z ,(cid:126)b T , η ) × Z RI (cid:48) / MOM O γ ( p R , b z ,(cid:126)b T , a, η ) , (9)where Z RI (cid:48) / MOM O γ is the RI (cid:48) / MOM renormalization fac-tor and p R denotes the matching scale introduced in theRI (cid:48) / MOM scheme. At all orders in perturbation theory,the scheme conversion factor R MS O γ cancels the depen-dence of Z RI (cid:48) / MOM O γ on p R and on the direction of (cid:126)b T (upto discretization artefacts).The authors have previously calculated Z RI (cid:48) / MOM O γ , andthereby Z MS O γ , by this approach in a quenched latticeQCD study [39]; those results are used for the numericalstudy in this work. The renormalization factor Z MS S doesnot 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 TMDPDFswhich gives the Collins-Soper kernel in Eq. (2), ˜∆ S andits renormalization factor Z MS S , which do not depend on b z , cancel between the numerator and denominator. Asa result, γ qζ ( µ, b T ) can be expressed in terms of the quasibeam function and its renormalization only, at the cost of introducing power-law divergences in η and b T separatelyin the numerator and denominator (divergences whichwere canceled by the quasi soft factor and its renormal-ization in the original expression for the kernel). More-over, to ensure that the renormalization and matchingbetween RI (cid:48) / MOM and MS is performed in the pertur-bative region, the scale b T must be taken to be muchsmaller than Λ − , a condition which does not permitan extraction of the Collins-Soper kernel at b T values inthe nonperturbative region. A perturbative renormal-ization matching scale b T = b RT (cid:28) Λ − in Eq. (8) can,however, be defined by exploiting the b T -independence of Z MS γ Γ ( µ, b z , a ), as described in Ref. [34]. In this approach,for the choice Γ = γ in Eq. (3), the Collins-Soper kernelcan be expressed as γ qζ ( µ, b T ) = 1ln( P z /P z ) ln (cid:34) C TMDns ( µ, xP z ) C TMDns ( µ, xP z ) × (cid:82) d b z e − ib z xP z P z lim a → η →∞ B MS γ ( µ, b z ,(cid:126)b T , a, η, P z ) (cid:82) d b z e − ib z xP z P z lim a → η →∞ B MS γ ( µ, b z ,(cid:126)b T , a, η, P z ) (cid:35) , (10)where a modified MS-renormalized quasi beam function B MSΓ has been defined as B MS γ ( µ, b z ,(cid:126)b T , a, η, P z ) = Z MS O γ ( µ, b z , b RT , a, η ) × ˜ R ( b T , b RT , a, η ) B bareΓ ( b z ,(cid:126)b T , a, η, P z ) . (11)Here, the additional factor ˜ R has been introduced intothe modified MS-renormalized quasi beam function tocompensate for the power-law divergences ∼ | b T − b RT | /a which would otherwise affect both the numerator anddenominator of Eq. (10):˜ R ( b T , b RT , a, η ) = Z RI (cid:48) / MOM O γ γ ( p R = ˜ p R , b z = 0 ,(cid:126)b T , a, η ) Z RI (cid:48) / MOM O γ γ ( p R = ˜ p (cid:48) R , b z = 0 ,(cid:126)b RT , a, η ) . (12)In this definition, fixed choices of ˜ p R , ˜ p (cid:48) R , and of the di-rections of (cid:126)b T and (cid:126)b RT , are taken. Since the factor ˜ R is in-dependent of b z , and thus cancels between the numeratorand denominator of Eq. (10), the specific choice of defini-tion will not affect the determination of the Collins-Soperkernel . In the numerical study in this work, an averageover (cid:126)b T and (cid:126)b RT orientations, and over several choices of˜ p R and ˜ p (cid:48) R , is performed in the same manner detailed inAppendix C in the numerator and denominator of ˜ R . The definition of ˜ R used here differs from that in Ref. [34] bythe omission of the quasi soft factor and by allowing ˜ p (cid:48) R to bedifferent from ˜ p R . In the numerical study presented here, a set of ten momenta p R with p R ranging from 5 . , as described in Ref. [39],are used to construct ˜ R . Label β a [fm] L × T κ n src n cfg E ×
64 0.1222 2 200TABLE I: The ensemble of quenched QCD gauge field con-figurations used in this work [41, 42]. The lattice spacing a isdetermined from an analysis of scale setting in Ref. [43], andthe lattice geometry parameters L and T are specified in unitsof a . For operator structures with Dirac index Γ = γ , n cfg configurations are analyzed, with n src source locations chosenon each. For other operator Dirac structures Γ (cid:54) = γ , a sub-set with 25 configurations is analyzed, with 1 source locationcomputed on each. Several observations are pertinent to the computationof 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 statewith the best signal-to-noise properties in a lattice QCDcalculation, e.g., for the pion. In a quenched calculation,a heavier-than-physical valence quark mass can be cho-sen for the same reason. Moreover, since although thekernel is state-independent, the power-corrections to thekernel are not, and so variation of the choice of exter-nal state, and external state momenta, provides a test ofsystematic effects in a numerical calculation. Second, theCollins-Soper kernel does not depend on the longitudinalmomentum fraction x or on the hadron momenta P zi , at O ( b T /η, / ( b T P z )). Although the truncation in the b z -space Fourier integral will induce oscillatory behavior in x -space, varying these parameters provides insight intothese additional systematic uncertainties.An alternative approach to extracting the Collins-Soper kernel by transforming the product of the match-ing coefficient and MS quasi beam function in Eq. (10)into a convolution integral in b z -space was advocated inRef. [34]. Appendix E provides an investigation of thisapproach and finds that it suffers from significant sys-tematic uncertainties. III. LATTICE QCD STUDY
The Collins-Soper evolution kernel is computed byEq. (10) in a lattice QCD calculation using a singlequenched ensemble, detailed in Table I. The calculationis undertaken on gauge fields that have been subjectedto Wilson flow to flow-time t = 1.0 [44], in order to in-crease the signal-to-noise ratio of the numerical results,and gauge-fixed to Landau gauge, in order to permit theuse of gauge non-invariant quark wall sources. Quasibeam functions are constructed for a pion external stateusing valence quark propagators that are computed withthe tree-level O ( a ) improved Wilson clover fermion ac-tion [45] and a κ value that corresponds to a heavy pionmass of 1.207(3) GeV. This choice may be made with-out introducing systematic bias, since the Collins-Soperkernel is independent of state. Three external state mo-menta are studied, (cid:126)P = P z (cid:126)e z with P z = n z π/L for □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ □ ○ △ FIG. 2: Effective energy function defined by Eq. (14) forpion states with momenta | (cid:126)P | = n z π/L . Shaded bands dis-play the result of single-exponential fits to the two-point cor-relation functions for each non-zero momentum, and a two-exponential fit at zero momentum; the number of states ineach fit is chosen to maximize an information criterion as de-scribed in the text, and the fit ranges shown correspond tothe highest-weight fits in the weighted average over successfultwo-point function fits as discussed in Appendix A. n z ∈ { , , } , corresponding to P z ∈ { . , . , . } GeV, allowing the kernel to be computed from threedifferent momentum ratios. To improve the overlapof boosted pion interpolating operators onto their re-spective ground states and improve statistical precision,a combination of wall sources and momentum-smearedsinks [46] are used to construct two-point and three-pointcorrelation functions.Bare quasi beam functions B bareΓ ( b z ,(cid:126)b T , a, η, P z ) areextracted for non-local quark bilinear operators (Eq. (6))with Wilson line staple geometries defined by stapleextents η ranging between 0.6 and 0.8 fm ( η/a ∈{ , , } ), and with staple widths and asymmetriescorresponding to | b T | and b z ranging from − ( η − a )to ( η − a ). In order for the mixing contributions toEq. (10) to be consistently included, bare quasi beamfunctions are computed for all Dirac operator structuresΓ. As detailed in the caption of Table I, however, lowerstatistics are used for operators with Dirac structuresΓ (cid:54) = γ , whose contributions to the Collins-Soper kernelare suppressed by the renormalization factors. Previ-ously, the 16-dimensional vector of MS renormalizationfactors Z MS O γ (cid:48) ( µ, b z ,(cid:126)b T , a, η ) was computed for the sameensemble and operator parameters as studied here [39],and those results are used in this work.The two-point correlation function for the pion, pro-jected to a given three-momentum (cid:126)P , is defined as: C ( t, (cid:126)P ) = (cid:88) (cid:126)x e i (cid:126)P · (cid:126)x (cid:104) | π (cid:126)P ,S ( (cid:126)x, t ) π † (cid:126)P ,W (0) | (cid:105) t (cid:29) −→ Z (cid:126)P aE (cid:126)P e − E (cid:126)P t + . . . , (13)where Z (cid:126)P denotes the combination of overlap fac-tors for the source and sink interpolation opera-tors and the ellipsis in Eq. (13) denotes contribu-tions from higher excitations, which are exponentiallysuppressed for large t and discussed further in Ap-pendix A. Wall-source interpolating operators π (cid:126)P ,W ( t ) = u ( t, (cid:126)P / γ d ( t, (cid:126)P /
2) are used as sources for correlationfunctions, where momentum projected quark fields aredefined by q ( t, (cid:126)P ) = (cid:80) (cid:126)x e i (cid:126)P · (cid:126)x q ( (cid:126)x, t ) for q = { u, d } .Momentum-smeared interpolating operators π (cid:126)P,S ( (cid:126)x, t ) = u S ( (cid:126)P / ( (cid:126)x, t ) γ d S ( (cid:126)P / ( (cid:126)x, t ) are used as sinks, where q S ( (cid:126)P ) ( (cid:126)x, t ) are quasi local smeared quark fields ob-tained through iterative application of the Gaussianmomentum-smearing operator defined in Ref. [46]. Inparticular, 50 steps of iterative momentum-smearingwith smearing radius ε = 0 .
25, as defined in Ref. [46],are used to construct momentum-smeared sinks for eachmomentum corresponding to n z ∈ { , , } . An effectiveenergy function that asymptotes to E (cid:126)P can be definedfrom the two-point correlation function by E eff (cid:126)P ( t ) = 1 a arccosh (cid:32) C ( t + a, (cid:126)P ) + C ( t − a, (cid:126)P )2 C ( t, (cid:126)P ) (cid:33) t (cid:29) −→ E (cid:126)P + . . . . (14)Two-point correlation functions for the three momentawhich are considered here are displayed in Fig. 2. The ex-tracted energies are slightly smaller than those obtainedwith the continuum dispersion relation, with relative de-viations from E (cid:126)P = (cid:113) m π + | (cid:126)P | ranging from 1 . n z = 2 to 3 . n z = 4. These deviations areconsistent with the expected size of lattice artifacts whichare neglected in this exploratory work. A. Quasi beam functions
Bare quasi beam functions can be computed fromthree-point correlation functions with insertions of thenon-local quark bilinear operators O i Γ ( b µ , z µ , η ), definedin Eq. (6). For the special case where pion momenta aretaken only in the z -direction (i.e., consistent with the def-inition of quasi beam functions in Eq. (4)), three-pointcorrelation functions are defined as C Γ ,i ( t, τ, b µ , a, η, (cid:126)P = P z (cid:126)e z )= (cid:88) (cid:126)x,(cid:126)z e i (cid:126)P · (cid:126)x (cid:104) | π (cid:126)P ,S ( (cid:126)x, t ) O i Γ ( b µ , ( (cid:126)z, τ ) , η ) π † (cid:126)P ,W (0) | (cid:105) t (cid:29) τ (cid:29) −−−−−→ Z (cid:126)P aE (cid:126)P e − E (cid:126)P t ˜ B Γ i ( b z ,(cid:126)b T , a, η, P z ) + . . . . (15)A ratio of three- and two-point correlation functions thenenables the bare isovector quasi beam functions of Eq. (5) □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - (a) Example of the computed bare quasi beam functions. ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ - - - (b) Comparison of the values of bare quasi beam functions fordifferent Dirac structures Γ, at the b z parameter indicated bythe orange dotted vertical line in subfigure (a). ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ - (c) Contribution to the renormalized quasi beam function B MS γ / ˜ R from each of the bare quasi beam functions shown insubfigure (b), as a fraction of the dominant contribution. Thelarge relative uncertainties result from the lower statistics usedto compute the off-diagonal renormalization factors and the barebeam functions with Dirac structures Γ (cid:54) = γ . FIG. 3: Examples of the extracted bare quasi beam functions B bareΓ ( b z ,(cid:126)b T = b T (cid:126)e x , a, η, P z = n z π/L ), defined in Eq. (5),for various parameter choices. Additional examples of thebare quasi beam functions with different parameter choicesare displayed in Appendix B. to be extracted: R Γ ( t, τ, b µ , a, η, P z )= C Γ ,u ( t, τ, b µ , a, η, P z (cid:126)e z ) − C Γ ,d ( t, τ, b µ , a, η, P z (cid:126)e z ) C ( t, P z (cid:126)e z ) t (cid:29) τ (cid:29) −−−−−→ B bareΓ ( b z ,(cid:126)b T , a, η, P z ) + . . . . (16)Constraining the bare quasi beam functions B bareΓ fromratios of two- and three- point functions R Γ for all staplegeometries (specified by { η, b µ } ), all Dirac structures Γ,and all momenta P z , considered in this work, requiresfits for a very large number of operators (35,660) to beperformed. These fits are automated using a fit proce-dure discussed in Appendix A. An example of the resultof these fits, for Γ = γ , and specific choices of b T and η , is given in Fig. 3(a); a second example figure hold-ing b z fixed, but showing all Dirac structures, is shownin Fig. 3(b). Additional examples of the real and imag-inary parts of the extracted bare quasi beam functionsare displayed in Appendix B.The bare quasi beam functions obtained by Eq. (16)are renormalized to the MS scheme by Eq. (11), us-ing renormalization factors Z MS O γ which were computedfor the same ensemble and operators as studied herein Ref. [39]. The fractional contributions to the renor-malized quasi beam function from the bare quasi beamfunctions with different Dirac structures Γ is shown inFig. 3(c), for a particular choice of parameters. The sizeof these contributions is observed to grow with increasing b T and with increasing ( η − b z ); while the relative magni-tudes of bare beam functions with different Γ do not varysignificantly with these parameters, the relative impor-tance of the off-diagonal renormalization factors variessignificantly as discussed in Ref. [39]. Across the param-eters studied, the combined contributions from mixing tothe renormalized quasi beam function with Dirac struc-ture Γ = γ are at the 5%–25% level. Calculations ofthe quasi beam functions, and the Collins-Soper evolu-tion kernel, to better than this precision thus require barequasi beam functions to be computed for several Diracstructures Γ.The functional dependence of the renormalized quasibeam function B MS γ ( µ, b z ,(cid:126)b T , a, η, b RT , P z ) (defined inEq. (11)) is shown in Fig. 4. The factor ˜ R ( (cid:126)b T , b RT , a, η )(Eq. (12)) was included in the definition of the renor-malized quasi beam function to cancel the dependenceof the bare beam function on η and on b RT . It is clearthat over choices of b RT within the perturbative region,this dependence is indeed removed to better than thestatistical uncertainties of this study. A weighted aver-age of the renormalized quasi beam function over theseparameters, as well as over different directions of (cid:126)b T , isthus taken as detailed in Appendix C to define averagedquasi beam functions B MS γ ( µ, b z , b T , a, P z ). Examples ofthe P z b z -dependence of the resulting quasi beam func-tions are shown in Fig. 5, and additional examples are shown in Appendix C. These are the key results used toextract the Collins-Soper kernel, as discussed in the nextsubsection. B. Collins-Soper kernel
Computing the Collins-Soper evolution kernel byEq. (10) requires taking the Fourier transform of the MS-renormalized quasi beam function with respect to b z . It isclear from the results shown in Fig. 5, however, that withthe parameter ranges explored in this study the Fouriertransform will suffer from significant truncation effectssince the quasi beam function is not yet consistent withzero within uncertainties at the largest b z values thatare used. For this reason, models are used to fit the P z b z -dependence of the lattice data for the quasi beamfunction before the Fourier transform is taken to eval-uate the Collins-Soper kernel. The results obtained bytaking discrete Fourier transforms instead are shown inAppendix D, and a discussion of what will be requiredfor future calculations to achieve robust results for theCollins-Soper kernel without this modeling step is pre-sented in Sec. IV.Two models of the P z b z -dependence of the quasi beamfunctions are considered, based on Hermite and Bernsteinpolynomial bases. The models are constructed to yield x -independent Collins-Soper kernels, as would be expectedin the absence of systematic artifacts, assuming the lead-ing order value for the perturbative matching coefficient,i.e., C TMDns = 1. Including higher orders in the match-ing factor while guaranteeing an x -independent kernelwould require more complicated functional forms to befit to the quasi beam functions to compensate for the x -dependence of the matching. It is expected that thematching uncertainties are small relative to the system-atic uncertainties inherent in introducing models for the P z b z -dependence of the quasi beam functions, and theseeffects are thus neglected in this work. While a model-independent result for the Collins-Soper kernel cannot beachieved from the data presented here, the comparisonbetween results obtained using the two different modelsconsidered nevertheless provides some indication of theseverity of the model-dependence, and the quality of fitsto these functional forms not including power correctionsalso provides a measure of their importance.The first functional form which is fit to the MS-renormalized quasi beam function is F Herm N ( P z , b z P z ; { a k } , γ, ω, σ )= N (cid:88) k =1 a k (cid:90) ∞−∞ dx e i ( b z P z ) x e − ( x − ω ) / σ ( P z x ) γ H k − ( x ) , (17)where H n ( x ) is the n -th Hermite polynomial. The fitparameter ω is taken to be complex, while the otherfree parameters are real. Allowing Im( ω ) (cid:54) = 0 allows the □ □ □○ ○ ○△ △ △ □ ○ △ (a) □ □ □○ ○ ○△ △ △ □ ○ △ (b) □ □ □ □○ ○ ○ ○ △ △ △ △ □ ○ △ (c) □ □ □ □ ○ ○ ○ ○ △ △ △ △ □ ○ △ (d) FIG. 4: Renormalized quasi beam function B MS γ ( µ, b z ,(cid:126)b T , a, η, b RT , P z ) in Eq. (11) (right column), and the same quantitydivided by the factor ˜ R ( b T , b RT , a, η ) in Eq. (12) (left column), similarly averaged, for various parameter choices. The horizontalshaded bands show the results of constant fits in b RT 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) □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○ ○ ○○○○○○○○ ○ ○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - - (b) FIG. 5: Averaged renormalized quasi beam function B MS γ ( µ, b z , b T , a, P z = n z π/L ) at small (a) and large (b) b T , afteraveraging over directions of (cid:126)b T , and weighted averaging over b RT and η , as detailed in Appendix C. Further examples of theaveraged renormalized quasi beam functions at different choices of b T are also given in Appendix C. △△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△ - - -
10 0 10 20 30 - - - (a) △△△△△ △△△△ △△△△△△△△△ △△△△△△△ △△□□□□□ □□□□ □□□□□□□□□ □□□□□□□ □□ - - -
10 0 10 20 30 - - (b) FIG. 6: Examples of fits to the averaged renormalized quasibeam functions B MS γ ( µ, b z , b T , a, P z ) using functional formsbased on Hermite and Bernstein polynomials (Eqs. (17-18)).Further examples of fits at different choices of the b T and P z parameters are shown in Appendix C. Fourier transform of F Herm N ( P z , b z P z ; { a k } , γ, ω, σ ) withrespect to b z P z to be complex, and correspondingly en-ables F Herm N ( P z , b z P z ; { a k } , γ, ω, σ ) to be an asymmetricfunction of b z P z . The real and imaginary parts of thequasi beam function are symmetric and antisymmetricfunctions of b z respectively in the η → ∞ limit; how-ever, the numerical results presented in this work showsignificant departures from these expectations, particu-larly for large b T , as shown in Fig. 5(b). The observedasymmetry could arise from finite-volume effects: effec-tive field theory calculations [49] have demonstrated thatfinite-volume effects for pion matrix elements of non-local operators with separation (cid:96) generically take theform e − m π ( L − (cid:96) ) . In this work, one therefore expects b z -dependent finite-volume effects of the form e − m π ( L − η + b z ) as well as additional b z independent finite-volume ef-fects. In addition, exponential dependence on b z couldarise from an imperfect cancellation between power-law-divergent lattice artifacts in B bareΓ ( b z ,(cid:126)b T , a, η, P z ) and Z MS O γ ( µ, b z ,(cid:126)b T , a, η ) ˜ R ( b T , b RT , a, η ). Taking Im( ω ) (cid:54) = 0 al- lows the fit form in Eq. (17) to include exponential depen-dence on b z and is found to significant improve the qualityof fits to the numerical results with large b T (cid:38) . < x < P z , and takes theform F Bern N ( P z , b z P z ; { a rn } , γ, A, B )= N − (cid:88) r =0 a r (cid:90) dx e i ( b z P z ) x x A (1 − x ) B ( P z x ) γ B r,N − ( x ) , (18)where B r,N − , for r ∈ { , . . . N − } are the N Bern-stein basis polynomials of degree N − b z is accommodated bytaking Im( a r ) (cid:54) = 0.Using either functional form, F Herm N or F Bern N , as amodel for B MS γ , and evaluating Eq. (10) with the tree-level matching factor C TMDns = 1, gives the result γ q, MS ζ = γ , where γ is the model parameter appearing in Eqs. (17)-(18). That is, the resulting Collins-Soper kernel is in-dependent of x by construction. The full procedure bywhich each functional form is fit to the numerical resultsfor the quasi beam function is described in Appendix C,and examples of the resulting fits are shown both in Fig. 6and in Appendix C. Briefly, the fits are undertaken simul-taneously at all P z and b z values for a given b T , and aninformation criterion is used to choose the model trunca-tion N for each fit. While both models fit the quasi beamfunction well within the range of P z b z values constrainedby the lattice data (with an average χ /N dof over all fitsof 0.9, tabulated in Appendix C), it is clear from Fig. 6that they correspond to substantially different modelsoutside this range.The Collins-Soper kernel determined from each set ofmodel fits is shown in Fig. 7. The results obtained us-ing the two model forms, i.e., the Hermite polynomialmodel, in which the quasi beam function has support on −∞ < x < ∞ , and the Bernstein polynomial model,with support on 0 < x <
1, are consistent. This en-couragingly suggests that γ q, MS ζ is well-constrained bythe numerical results at the P z and b z values of this cal-culation, and that the model-dependence introduced inthe Fourier transform is relatively mild. Perturbativeresults for the 0-flavor Collins-Soper kernel [47, 48] are FIG. 7: Collins-Soper evolution kernel obtained using fits to the renormalized quasi beam functions based on Hermite andBernstein polynomial bases (Eqs. (17-18)), computed as described in the text. The background shading density is proportionalto 1 / ( b T P z ) + b T /η , indicating regions of greater and lesser sensitivity to power corrections which are not included in theuncertainties presented. The black dotted, dashed and solid lines show perturbative results for the 0-flavor Collins-Soper kernelup to three-loop order [47, 48]. Perturbative results become singular at b T ∼ .
25 fm because they reach the Landau poleassociated with Λ MS ,N f =0QCD = 639 MeV. (a) (b) FIG. 8: Fractional truncation effects in the MS-renormalized quasi beam functions, defined by Eq. (19), evaluated at x = 0 . b T values shown. The red vertical line denotes the maximum b z used in this study; the vertical axis rangecorresponds to the P z range of this study. It is noteworthythat the lattice QCD results for γ q, MS ζ obtained here areconsistent with perturbative calculations of the 0-flavorCollins-Soper kernel [47, 48] in the region b T ∼ . b T < . / ( b T P z ),which have not been estimated here.Although the Collins-Soper kernel shown in Fig. 7 hasbeen obtained in a 0-flavor calculation, it can also becompared qualitatively with the results of fits to experi-mental data, in which several different parametrizationsof the nonperturbative behavior of the kernel at large b T have been used. In early fits to Drell-Yan data [4–6],the Collins-Soper kernel was parametrized as a quadraticfunction in b T in the nonperturbative region. It was laterfound in Refs. [53, 54], however, that these fits cannotdescribe SIDIS data. More recently, it has been arguedthat γ q, MS ζ should approach a constant as b T → ∞ ; phe-nomenological fits to Drell-Yan data under this assump-tion suggest that this constant is ∼ − . b T . The results of this numerical study are qualitativelyconsistent with constant or linear behavior of the kernelin the nonperturbative region. Once lattice QCD resultswith controlled systematic uncertainties are available, itwill be possible to test these and other phenomenologi-cal expectations for the large- b T behavior of the Collins-Soper kernel with QCD predictions, and begin incorpo-rating lattice QCD constraints in phenomenological anal-yses. The requirements for a fully controlled lattice QCDdetermination of the Collins-Soper kernel are discussed inthe following section. IV. OUTLOOK
This manuscript presents an exploratory calculationof the nonperturbative Collins-Soper kernel in quenchedlattice QCD based on the method developed in Refs. [32–34]. In this approach, the kernel is computed fromquasi beam functions defined from matrix elements ofquark bilinear operators with staple-shaped Wilson linesin boosted hadron states. These beam functions are In this work, α MS ,N f =0 s is determined by evolving α MS ,N f =5 s ( µ = M Z ) from Ref. [51] to lower scales usingthe four-loop β function [52], integrating out bottom and charmquarks, and finally matching α MS ,N f =0 s ( µ ) to α MS ,N f =3 s ( µ ) atthe scale µ = 2 GeV where Z MS O ΓΓ (cid:48) ( µ, b z , b T , a, η ) is calculatedin Ref. [39]. This procedure gives the result Λ MS ,N f =0QCD = 639MeV, which determines α MS ,N f =0 s ( µ ) at all µ ; throughout thiswork α MS ,N f =0 s ( µ = 2 GeV) = 0 . renormalized to the MS scheme via the RI (cid:48) / MOM pre-scription, and a ratio of Fourier-transformed quasi beamfunctions at different hadron boost momenta determinesthe Collins-Soper kernel. In this calculation, the kernelis extracted over a range of scales b T ∈ (0 . , .
8) fm. Thefinal results presented here rely on modeling the b z -spacequasi beam functions to control truncation effects in theFourier transform. Nevertheless, the results are robustunder the variation of models considered here.For a future controlled and model-independent deter-mination of the Collins-Soper kernel by this method, sev-eral improvements will be essential. Critically, larger lat-tice volumes must be studied; the overwhelming system-atic in this calculation arises from modeling to facilitatethe Fourier transform, which is necessary because of trun-cation effects suffered due to the limited b z range overwhich quasi beam functions could be computed on thelattice volume used here. One measure of truncation ef-fects is given by the model truncation error defined as δ trunc B MS = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) DFT b z max (cid:2) F Herm N (cid:3) ( P z , x )FT ∞ (cid:2) F Herm N (cid:3) ( P z , x ) − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (19)where the truncated discrete Fourier transform (DFT)and untruncated Fourier transform are defined asDFT b z max (cid:2) F Herm N ( P z , b z P z ) (cid:3) ( P z , x ) (20a)= P z b z max (cid:88) b z = − b z max e − ixP z b z F Herm N ( b z P z ) , FT ∞ (cid:2) F Herm N ( b z P z ) (cid:3) ( P z , x ) (20b)= (cid:90) ∞−∞ P z db z π e − ixP z b z F Herm N ( b z P z ) . Equation (19) is evaluated by applying both the DFTand untruncated Fourier transform to the best-fit model F Herm N ( P z , b z P z ; { a k } , γ, ω, σ ) obtained by fitting the lat-tice QCD results at fixed b T , as described in Sec. III.This provides an estimate, based on the model, of theeffects of extending b z max to larger values of η than thoseused for the numerical calculations in this work. Fig-ure 8 shows δ trunc B MS ( x ) for x = 0 . b T over the range of P z values used in this work, andat values of b z max both comparable to the value usedhere and considerably larger. For the value of b z max used here, truncation effects in the DFT results are O (1)for large b T and prevent a reliable determination of theCollins-Soper kernel using a DFT approach. Results withsignificantly larger P z b z max , however, could be used toobtain 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 ofthe present study, quasi beam function calculations with( b z max , P z ) ∼ (2 . , . b T scales across therange of those studied in this work.In addition to truncation artifacts, extractions of theCollins-Soper kernel by the method pursued here suffer1from power corrections at O ( b T /η, / ( b T P z )) [33, 36].These effects could not be resolved by model fits in thisstudy, and as such, the coefficients of these power cor-rections could not be constrained. Nevertheless, largerphysical lattice volumes, as needed to reduce truncationartifacts, will simultaneously enable O ( b T /η, / ( P z η ))effects to be reduced by allowing larger staple extents η to be investigated at fixed b T . Larger boost momen-tum, again needed to control truncation effects, will alsosimultaneously enable control over power corrections of O (1 / ( b T P z )). These artifacts make comparison of theCollins-Soper kernel extracted by the method pursuedhere with perturbative predictions, which are accurate inthe region b T (cid:28) Λ − , challenging. Precise comparisonsin this region will be an important test of systematics inthe lattice QCD approach. The infinite volume and con-tinuum limits must also ultimately be taken for a fullycontrolled result.Future studies would also gain significantly by exploit-ing the state-independence of the Collins-Soper kernel toobtain multiple constraints on the kernel from the samecalculation and thus test systematic effects. An alterna-tive, complementary, approach to extracting the Collins-Soper kernel from Eq. (2) was proposed in Ref. [37]. Thisstrategy uses the Mellin moments of the expressions sothat one only needs to calculate the quasi beam functionor its derivatives at b z = 0, which reduces the computa-tional cost and has the advantage that renormalizationfactors cancel in the ratio. This approach also, however,requires a nontrivial integration over the TMDPDF thatis extracted from experiments over a limited kinematicrange. Comparison of results of the two approaches willalso be valuable in future calculations.Despite the significant challenges described above,the results presented here suggest that controlled first-principles calculations of the Collins-Soper kernel at non-perturbative scales as large as b T ∼ b T ∈ (0 . , .
0) fm will be sufficient to differentiate differentmodels of the Collins-Soper kernel and will thereby pro-vide important input for fitting low-energy SIDIS data.Ultimately, larger values of b T , e.g. b T (cid:46) Acknowledgments
The authors thank Will Detmold, Markus Ebert,Michael Engelhardt, Andrew Pochinsky, and Iain Stew-art for helpful discussions, and Michael Endres for pro- viding the gauge field configurations used in this project.Calculations were performed using the Qlua [56] andChroma [57] software libraries. This work is supported inpart by the U.S. Department of Energy, Office of Science,Office of Nuclear Physics, under grant Contracts No. de-sc0011090, No. de-sc0012704 and within the frameworkof the TMD Topical Collaboration. PES is additionallysupported by the National Science Foundation under CA-REER Award No. 1841699. MLW was additionally sup-ported in part by an MIT Pappalardo fellowship. Thisresearch used resources of the National Energy ResearchScientific Computing Center (NERSC), a U.S. Depart-ment of Energy Office of Science User Facility operatedunder Contract No. DE-AC02-05CH11231, and the Ex-treme Science and Engineering Discovery Environment(XSEDE), which is supported by National Science Foun-dation grant No. ACI-1548562. This manuscript has beenauthored by Fermi Research Alliance, LLC under Con-tract No. DE-AC02-07CH11359 with the U.S. Depart-ment of Energy, Office of Science, Office of High EnergyPhysics.
Appendix A: Fits to three- and two-point functions
As shown in Eq. (16), ratios of three-point and two-point correlation functions asymptote to the bare quasibeam function in the double limit { τ, t − τ } → ∞ . Ratioscomputed at finite t and τ , however, will include contri-butions from matrix elements in excited states. Two-point correlation functions have the spectral representa-tion C ( t, (cid:126)P ) = (cid:88) n Z n(cid:126)P aE n(cid:126)P (cid:16) e − E n(cid:126)P t + e − E n(cid:126)P ( T − t ) (cid:17) , (A1)where n labels QCD energy eigenstates with the quan-tum numbers of a pion with momentum (cid:126)P , E n(cid:126)P is theenergy of state n , (cid:113) Z n(cid:126)P = 2 E n(cid:126)P (cid:104) | π (cid:126)P | n (cid:105) are overlapfactors for the interpolating operator π (cid:126)P onto state n ,and thermal effects arising from the finite Euclideantime extent T are included. Dependence of Z n(cid:126)P on thetype of source/sink interpolating operators used (wallor momentum-smeared) is suppressed in Eq. (A1) andthroughout this section. Fits of lattice QCD two-pointcorrelation function results to Eq. (A1) can be used toextract the ground-state energies E (cid:126)P shown in Fig. 2, aswell as excited-state energies and overlap factors.Three-point functions have an analogous spectral rep-resentation2 C Γ ,i ( t, τ, b µ , a, η, (cid:126)P ) = (cid:88) m,n (cid:113) ( Z m(cid:126)P ) ∗ Z n(cid:126)P aE m(cid:126)P E n(cid:126)P e − E m(cid:126)P ( t − τ ) e − E n(cid:126)P τ (cid:10) m (cid:12)(cid:12) O i Γ ( b µ , z µ , η ) (cid:12)(cid:12) n (cid:11) , (A2)where n, m index energy eigenstates. Combing Eq. (A1) and Eq. (A2) and isolating the ground-state contributionsyields a spectral representation for ratios of three- and two-point functions: C Γ ,i ( t, τ, b µ , a, η, (cid:126)P ) C ( t, a, (cid:126)P ) × (cid:88) n Z n(cid:126)P E n(cid:126)P (cid:32) Z (cid:126)P E (cid:126)P (cid:33) − (cid:16) e − ( E n(cid:126)P − E (cid:126)P ) t + e − E n(cid:126)P T e ( E n(cid:126)P + E (cid:126)P ) t (cid:17) = B bareΓ ( b z ,(cid:126)b T , a, η, P z ) (cid:88) ( n,m ) (cid:54) =(0 , (cid:113) ( Z m(cid:126)P ) ∗ Z n(cid:126)P ( E (cid:126)P ) E m(cid:126)P E n(cid:126)P Z (cid:126)P e − ( E m(cid:126)P − E (cid:126)P )( t − τ ) e − ( E n(cid:126)P − E (cid:126)P ) τ (cid:10) m (cid:12)(cid:12) O i Γ ( b µ , z µ , η ) (cid:12)(cid:12) n (cid:11)(cid:10) (cid:12)(cid:12) O i Γ ( b µ , z µ , η ) (cid:12)(cid:12) (cid:11) . (A3)After determining the spectral quantities appearing onthe left-hand-side of Eq. (A3) from fits to lattice QCDresults for C , where in practice the sum over statesis truncated at n = N states as discussed below, the barequasi beam functions and other parameters appearingon the right-hand-side of Eq. (A3) can be determinedfrom fits to lattice QCD results for three-point to two-point function ratios. Fitting directly to these ratioshas the advantages that ground-state overlap factors can-cel exactly between three- and two-point functions andthat correlated ratios are determined more precisely thanthree-point functions alone. Including the additional fac-tor on the left-hand-side of Eq. (A3), which depends onlyon energies and overlaps obtained in two-point functionfits, removes the need to model excited-state contami-nation in the two-point function during fits to the ratio(which would require fitting several additional parame-ters entering χ -minimization nonlinearly) without spoil-ing these correlations.Three-point correlation functions are computed for sixsource/sink separations t/a ∈ { , , , , , } andall operator insertion points 0 < τ < t . Signal-to-noiseratios of two-point and three-point correlation functionsare proportional to e − ( E (cid:126)P − m π ) t , where m π is the pionmass, and for n z ≥ t/a ∈ { , , , } are used in fits. Correlated χ -minimization fits of two-point functions to Eq. (A1), followed subsequently by fitsto Eq. (A3), are performed in an automated manner asfollows: • The minimum source/sink separation t min are var-ied over the range 2 ≤ t min ≤ t max − t plateau , where t max is chosen to be the largest t for which thesignal-to-noise ratio of C ( t, a, (cid:126)P ) is greater thana fixed value (a threshold of 2 is used in the re-sults presented here) and t plateau is a free parame-ter specifying the minimum number of points in afit (results presented here use t plateau = 3). The re- striction t min ≥ t min within this range, correlated χ -minimizationfits to Eq. (A1) are performed using two-point func-tion results with t min ≤ t ≤ t max . The two-pointfunction fitting procedure is identical to that de-scribed in Appendix B of Ref. [58]. To summa-rize, one-state fits are performed first, followed bytwo-state fits. If the Akaike information criterion(AIC) [59] is improved sufficiently by the additionof a second state to the fit (a threshold of ∆AIC < − N dof , where N dof is the number of degrees offreedom in the fit, is used in the final results), thena three-state fit is performed and the same crite-rion is used to judge whether the three-state fit ispreferred. This procedure is repeated until addingadditional states does not sufficiently improve thefit, in order to select the optimal number of statesto include in the fit for each t min . The best-fit pa-rameters are determined using nonlinear optimiza-tion for the energies E n(cid:126)P , with linear systems ofequations solved to determine Z n(cid:126)P at each iteration.Covariance matrices are determined using optimalshrinkage [60] as described in Refs. [58, 61] in orderto reduce finite-statistic effects leading to poorlyconditioned sample covariance matrices. Severalchecks on numerical χ optimization described inRef. [58] are then performed to verify the reliabil-ity of the fit. If these checks are passed, an accept-able two-point function model has been found forthis choice of t min , and three- to two-point func-tion ratios are subsequently analyzed using fits toEq. (A3). • The minimum source/operator and source/sinkseparations corresponding to τ ∈ [ τ min , τ max ] arevaried over the ranges 2 ≤ τ min ≤ ( t min − t plateau ) / ≤ t − τ max ≤ ( t min − t plateau ) /
2. Three-point to two-point ratios using all available t ∈ t min , t max ] and τ ∈ [ τ min , τ max ] are multipliedby the factor in brackets on the left-hand-side ofEq. (A3). A correlated χ -minimization fit is per-formed to extract the parameters on the right-hand-side of Eq. (A3) using the same methods ap-plied to two-point functions. The excited-state ma-trix elements appearing in Eq. (A3) enter the χ function linearly, and their optimal values are de-termined by solving a linear system of equationsat each step of iterative nonlinear optimization forthe energies appearing in Eq. (A3) as done in vari-able projection methods [62, 63]. Because the low-lying spectrum is imperfectly modeled by few-statefits, the energies extracted from fits to Eq. (A3)are not constrained to identically equal the ener-gies extracted from fits to Eq. (A1) (although thespectrum determined from Eq. (A1) is used to pro-vide initial conditions for nonlinear optimization).Optimal shrinkage is used to define the covariancematrix. The best-fit ground-state matrix elementdefines ( B bareΓ ) f for fit range choice f defined by t min , τ min and τ max . Fits where two solvers dis-agree on ( B bareΓ ) f by more than a specified toler-ance (10 − is used in final results) are discarded inanalogy to the reliability checks applied to fits totwo-point functions [58]. • Confidence intervals for ground-state matrix ele-ments and other fit parameters are determined us-ing bootstrap resampling; see Refs. [64, 65] for re-views. Fits to Eq. (A3) are repeated N boot times( N boot = 200 is used in final results) using ensem-bles constructed by randomly resampling with re-placement from the two- and three-point functionsin a correlated manner. Statistical uncertainties onfit parameters are obtained from empirical confi-dence intervals of bootstrap fit results as detailedin Ref. [58]. The 68% confidence interval defines( δB bareΓ ) f . Further reliability checks are applied:the median of the bootstrap distribution is veri-fied to be within a specified tolerance of ( B bareΓ ) f (2 σ is used in final results), and uncorrelated fit re-sults are verified to be within a specified toleranceof ( B bareΓ ) f (5 σ is used in the final results). Fitspassing all reliability checks define an ensemble of f = { , . . . , N success } successful fit range choices. • A weighted average of ground-state matrix elementresults from all successful fits is used to determinethe final results. Each fit result ( B bareΓ ) f providesan unbiased estimate of the bare quasi beam func-tion (in the infinite-statistics limit), and the rela-tive weights between successful fits are arbitrary inthe large statistics limit. Given a finite statisticalensemble, it is advantageous to define a weightedaverage that penalizes fits with worse goodness-of-fit and fits with larger uncertainties. The weightedaverage procedure used in Refs. [58, 61] is followed: averages are defined by B bareΓ ≡ N success (cid:88) f =1 w f ( B bareΓ ) f , (A4a)( δ stat B bareΓ ) ≡ N success (cid:88) f =1 w f (cid:0) ( δB bareΓ ) f (cid:1) , (A4b)( δ sys B bareΓ ) ≡ N success (cid:88) f =1 w f (cid:0) ( B bareΓ ) f − B bareΓ (cid:1) , (A4c)( δB bareΓ ) ≡ (cid:113) ( δ stat B bareΓ ) + ( δ sys B bareΓ ) , (A4d)where the weights w f are defined as w f ≡ ˜ w f (cid:80) N success f =1 ˜ w f , (A5a) (cid:101) w f ≡ p f (cid:0) ( δB bareΓ ) f (cid:1) − (cid:80) N success f (cid:48) =1 p f (cid:48) (cid:0) ( δB bareΓ ) f (cid:1) − . (A5b)Here, p f = Γ( N dof / , χ f / / Γ( N dof /
2) where Γis the gamma function (not to be confused withthe Dirac spinor index elsewhere in this work) inorder to penalize fits with large χ /N dof and large( δB bareΓ ) f . See Refs. [58, 61] for further discussionof this procedure.The resulting average values and uncertainties com-puted as in Eq. (A4d) define the bare quasi beam func-tions B bareΓ ( b z ,(cid:126)b T , a, η, P z ) and δB bareΓ ( b z ,(cid:126)b T , a, η, P z )used in this work. A representative sample of three- andtwo-point function ratio fit results are shown in Fig. 9 forthe smallest and largest momenta used in this work andfor both small and large Wilson-line extents among theset studied here. Appendix B: Bare quasi beam functions
Additional examples of the real and imaginaryparts of the extracted bare quasi beam functions B bareΓ ( b z ,(cid:126)b T , a, η, P z ), defined in Eq. (5), and determinedas described in Sec. III, are shown for various parameterchoices in Fig. 10. A general trend can be observed thatat increasing b T both the real and imaginary parts of thefunctions become more asymmetric in P z b z . This asym-metry arises primarily from linear divergences in the barequasi beam function that are canceled when the renor-malization factors Z MS O γ Γ are included, as can be seen bycomparing Fig. 10 and Fig. 11. As discussed in Sec. III B,residual b z asymmetries in B MS γ visible in Fig. 11 couldarise from finite-volume effects coupled with imperfectcancellation of exponential b z dependence between thebare quasi beam functions and the renormalization fac-tors.4 □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ □ ○ △▽ ◇ ▯ - - - □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ □ ○ △▽ ◇ ▯ - - - □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ □ ○ △▽ ◇ ▯ - - - □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ ▯ □ ○ △▽ ◇ ▯ - - - - - - □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ □ ○ △ ▽ - - - - - □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ □ ○ △ ▽ - - - - □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ □ ○ △ ▽ - - - - - □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ □ ○ △ ▽ - - - - - - - - FIG. 9: Examples of fits to the ratio of three- and two-point functions R Γ ( t, τ, b µ , a, η, P z ) (Eq. (16)), obtained as describedin the text. Shaded bands of the same colors as the points show 68% bootstrap confidence intervals of the τ and t -dependentfits from the fit range (specifically the choice of t min , τ min , and τ max ) that had the highest weight in the weighted average ofsuccessful 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 different fit rangechoices. □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - □□□□□□□ □□□□□□□□□□□□□□ □□□□□□○○○○○○○○○ ○○○○○○○○○○ ○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - □□□□□□□□□□□□□□ □□□□□□□□□□□□□ ○○○○○○○○○○○○○○ ○○○○○○○○○○ ○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - FIG. 10: Bare quasi beam functions B bareΓ ( b z ,(cid:126)b T = b T (cid:126)e x , a, η, P z = n z π/L ), defined in Eq. (5), for various parameter choices. Appendix C: Renormalized beam functions
The renormalized quasi beam function is computed bycombining the bare quasi beam function determined fromthree- to two-point function ratios as described in Ap-pendix A with the renormalization factors computed inRef. [39], as shown in Eq. (11). The uncertainty on therenormalized quasi beam function is obtained by combin-ing the total uncertainties of the bare quasi beam func-tion and the renormalization factors in quadrature. Re- sults are computed using two different staple orientationscorresponding to (cid:126)b T = b T (cid:126)e x and (cid:126)b T = b T (cid:126)e y . Interchang-ing these orientations (cid:126)e x ↔ (cid:126)e y is an exact symmetry of B MS γ but is not a symmetry of B bareΓ for some Γ, andthus results with both orientations can be averaged onlyafter renormalization. A weighted average [66] is used tocombine the B MS γ results with (cid:126)b T = b T (cid:126)e x and (cid:126)b T = b T (cid:126)e y by B (cid:48) MS γ ( µ, b z , b T , a, η, b RT , P z ) ≡ (cid:88) k =1 w k B MS γ ( µ, b z , b T (cid:126)e k , a, η, b RT , P z ) ,δ stat B (cid:48) MS γ ( µ, b z , b T , a, η, b RT , P z ) ≡ (cid:88) k =1 w k δB MS γ ( µ, b z , b T (cid:126)e k , a, η, b RT , P z ) ,δ sys B (cid:48) MS γ ( µ, b z , b T , a, η, b RT , P z ) ≡ (cid:88) k =1 w k (cid:104) B (cid:48) MS γ ( µ, b z , b T , a, η, b RT , P z ) − B MS γ ( µ, b z , b T (cid:126)e k , a, η, b RT , P z ) (cid:105) ,δB (cid:48) MS γ ( µ, b z , b T , a, η, b RT , P z ) ≡ δ stat B (cid:48) MS γ ( µ, b z , b T , a, η, b RT , P z ) + δ sys B (cid:48) MS γ ( µ, b z , b T , a, η,(cid:126)b RT , P z ) , (C1)where the weights are chosen to sum to unity and to be proportional to the inverse variance of each result: w k ≡ ˜ w k (cid:80) k =1 ˜ w k , ˜ w k ≡ δB MS γ ( µ, b z , b T (cid:126)e k , a, η,(cid:126)b RT , P z ) . (C2)As shown for one example in Fig. 4, the renormalized quasi beam functions do not depend on η or b RT withinuncertainties. The formal extrapolation to η → ∞ , and an average over possible choices of b RT in the window a (cid:28) b RT (cid:28) Λ − QCD , are thus implemented with an analogous weighted average: B MS γ ( µ, b z , b T , a, P z ) ≡ (cid:88) η/a ∈{ , , } (cid:88) ( b RT /a )=2 w η,b RT B (cid:48) MS γ ( µ, b z , b T , a, η, b RT , P z ) ,δ stat B MS γ ( µ, b z , b T , a, P z ) ≡ (cid:88) η/a ∈{ , , } (cid:88) ( b RT /a )=2 w η,b RT δB (cid:48) MS γ ( µ, b z , b T , a, η, b RT , P z ) ,δ sys B MS γ ( µ, b z , b T , a, P z ) ≡ (cid:88) η/a ∈{ , , } (cid:88) ( b RT /a )=2 w η,b RT (cid:20) B MS γ ( µ, b z , b T , a, η, b RT , P z ) − B (cid:48) MS γ ( µ, b z , b T , a, η, b RT , P z ) (cid:21) ,δB MS γ ( µ, b z , b T , a, P z ) ≡ δ stat B MS γ ( µ, b z , b T , a, P z ) + δ sys B MS γ ( µ, b z , b T , a, P z ) , (C3)where the weights are w η,b RT ≡ ˜ w η,b RT (cid:80) η/a ∈{ , , } (cid:80) b RT /a )=2 ˜ w η,b RT , ˜ w η,b RT ≡ δB (cid:48) MS γ ( µ, b z , b T , a, η, b RT , P z ) . (C4)The resulting renormalized quasi beam functions are shown in Fig. 11 along with fits to the Hermite and Bernsteinfunctional forms shown in Eqs. (17-18).Joint fits to Eqs. (17-18) for the real and imaginary parts of B MS γ ( b z , b T , a, P z ) for fixed b T , all b z ∈ [ − η +7 □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - □□□□□ □□□□□□□□□□□□□□□□□ □□□□□○○○○○○○○ ○○○○○○○○○○○ ○○○○○○○○ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - □□□□□□□□□□□□□ □□□□□□□□□□□□□□○ ○○○○○○○○○○○○○ ○○○○○○○○○○○ ○○ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - - □□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - - FIG. 11: Averaged renormalized quasi beam function B MS γ ( b z , b T , a, P z = n z π/L ), for various parameter choices. The bandsshow fits to a functional form based on Bernstein polynomials (Eq. (18)), as described in the text. a, η − a ], and all P z corresponding to n z ∈ { , , } ,are performed using uncorrelated χ -minimization Inorder to determine the appropriate polynomial to usein Eqs. (17-18), the AIC is employed. Fits to Eq. (17)with Hermite polynomials of degree N = { , , } areperformed, and minimum-AIC fits with b T /a ∈ [1 , χ /N dof as given in Ta-ble II. Fits to Eq. (18) with Bernstein polynomials ofdegree N = { , , , } are similarly performed, andminimum-AIC fits are obtained with minimum χ /N dof as given in Table III. Both fit forms describe the numer-ical quasi beam function results, although the Hermiteform achieves a lower χ /N dof in particular for small b T .Substituting either of the fit forms in Eqs. (17-18) intoEq. (10) and analytically performing the Fourier trans-forms gives the result γ q, MS ζ = γ where γ is the fit pa-rameter appearing in Eqs. (17-18) (neglecting one-loopmatching effects as discussed in Sec. III B). Best-fit val-ues for γ from these fits therefore provide determinationsof γ q, MS ζ . The statistical uncertainties of these determi-nations are calculated by bootstrap resampling renor-malized quasi beam function results from a Gaussiandistribution with mean B MS γ and width δB MS γ , refittingeach resampled ensemble, and taking 68% empirical con-fidence intervals of γ in the resulting fits. This procedureprovides the results for γ q, MS ζ for the Hermite and Bern-stein fit forms shown in Fig. 7. Appendix D: CS kernel from DFT
In this appendix, a strategy to extract the Collins-Soper kernel from the quasi beam function via theDFT method, as proposed in Refs. [32–34], is discussed.Naively taking a DFT of the quasi beam function ob-tained in this study at different momenta (Fig. 12),and extracting the Collins-Soper evolution kernel us-ing Eq. (10), yields the results shown in Fig. 12(b).Clearly, the results from ratios formed using three dif-ferent momentum pairs have significant x -dependence,and are different from each other by over 3 σ at the peakvalues; convergence is not apparent. For this reason, un-truncated Fourier transforms using models of the large P z b z behavior of the quasi beam function are used as de- Uncorrelated fits are performed because correlations betweenquasi beam functions with different staple geometries are notaccounted for in the total statistical plus systematic uncertain-ties of each { Γ , b z ,(cid:126)b T , η } , which are determined using weightedaverages over multiple fit range choices as described in Ap-pendix A. Although nonzero correlations exist for different ma-trix elements computed using the same gauge field configura-tions, the systematic uncertainties arising from using uncorre-lated χ -minimization are expected to be small compared tothose inherent in modeling the b z dependence of the quasi beamfunction. □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△ □ ○ △ (a) The points show the DFT of the averaged renormalizedquasi beam function B MS γ ( b z , b T , a, P z = n z π/L ). The bandsshow the results of an untruncated Fourier transform applied tothe Bernstein polynomials fit to the data (Eq. (18)). □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△ □ ○ △ - - - (b) Collins-Soper kernel extracted based on the DFT of thequasi beam function shown in (a). The solid grey line shows theresult obtained using the Bernstein polynomial model fit. FIG. 12: Fourier transformed quasi beam functions and aDFT calculation of the Collins-Soper kernel. □□□□□□□□ □□□□□□□□□□□□□□□□□□□□ □□□□□□□□○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△ □ ○ △ - - - FIG. 13: Extraction of the Collins-Soper kernel fromthe toy model of Eq. (D1) with the DFT method at { P z , P z } = { , }× (2 π/ γ q, toy ζ ( b T ) = − . b T /a with the µ -dependence suppressed. scribed in the main text, rather than the DFT approach.9 b T N χ /N dof γ -0.27(2) -0.29(2) -0.25(3) -0.24(7) -0.28(3) -0.33(3) -0.38(4) -0.42(3) -0.45(3) -0.44(4) -0.44(4) -0.46(4) -0.50(7)TABLE II: Order N and minimum χ /N dof obtained in fits of the quasi beam functions to the Hermite polynomial basedmodel form of Eq. (17), performed as described in the text. b T N χ /N dof γ -0.14(3) -0.19(3) -0.22(2) -0.22(2) -0.27(2) -0.32(3) -0.41(3) -0.43(4) -0.47(5) -0.45(8) -0.58(12) -0.43(11) -0.73(28)TABLE III: As in Table II, for the Bernstein polynomial based model form of Eq. (18).
20 25 30 35 40 45 50 - - - - - - FIG. 14: Collins-Soper kernel from the toy model of Eq. (D1)with the DFT method at { P z , P z } = { , } × (2 π/
32) anddifferent η values. The Collins-Soper kernel is determined bytaking the average of the central peak and trough values near x = 0 .
5, and the result converges to the original value withincreasing η . The instability in Fig. 12(b) can be understood as aconsequence of the limited range of b z and η consideredin this study. From Eq. (10), the ratio of quasi beamfunctions should stabilize in a certain x -region for a givenmomentum pair { P z , P z } . Given a limited range of b z inthe Fourier transform, this stabilization can be expectedto be robust for values of x ∼ . x ∼ . | b z | ≤ ( η − a )in this calculation, the Fourier transform will introducean oscillatory term to the quasi beam function with fre-quency P z ( η − a ). The same issue has been encounteredin lattice QCD calculations of collinear PDFs [67–70],and certain model assumptions have been suggested toavoid this challenge [71, 72]. For P z = 6 π/L , η = 12 a ,as in this calculation, the period of this oscillation is T = 2 π/ ( P z ( η − a )) ∼ .
0, and the oscillatory behavioris not apparent within the region 0 < x <
1, as shownin Fig. 12(a). The oscillatory behavior will persist in ra-tios of quasi beam functions at different momenta P z , P z , as an interference between oscillations with frequen- cies P z ( η − a ) and P z ( η − a ). As a result, the Collins-Soper kernel extracted from the peak can be shifted sig-nificantly, which adds an important systematic error tonumerical calculations via this approach.In future calculations with increased ranges of η or P z ,such that there are more rapid oscillations of the DFTs ofquasi beam functions within the range 0 < x <
1, this ap-proach 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 atleast two cycles, which would allow it to be determinedby taking an average of the central local maximum andminimum within the oscillating region. To illustrate thispoint, a toy model for the quasi beam function ˜ B ns in x -space is constructed:˜ B ns ( x, b T , µ, P z ) = 10 C ns ( µ, xP z ) x (1 − x ) × exp (cid:104) − . b T /a ) ln ( xP z a ) − . b T /a ) (cid:105) , (D1)where x ∈ [0 , P z and b T are in lattice units. TheMS scale is set to µ = 2 . a = 0 .
06 fm and L = 32 a .To study the oscillatory behavior in this toy model, theinverse FT of the quasi beam function ˜ B ns ( x, b T , µ, P z )is taken first. Then, a DFT of the truncated quasibeam functions back to x -space is performed for η = { a, a, a } , and the Collins-Soper kernel is computedusing Eq. (10). The results are shown in Fig. 13. It is ap-parent that the kernel suffers from oscillations due to thetruncated DFT, and for η = 12 a the shape of the curve isqualitatively similar to those in Fig. 12(b). Moreover, thepeaks or local maximums around x = 0 . η choices. Nevertheless, for η = 24 a and 36 a , takingthe average of the central peak and trough values pro-vides a close approximation to the Collins-Soper kernel,as shown in Fig. 14. With more rapid oscillations, thisaveraging method will lead to more accurate results. Fu-ture calculations with larger lattices sizes and higher pionmomenta will thus likely enable reliable determination ofthe Collins-Soper kernel with the DFT method, althoughthe toy model results shown in Fig. 14 suggest that very0 □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - - - - - - □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ □ ○ △ - - - - - - - - - FIG. 15: Collins-Soper kernel for the toy model of Eq. (D1)determined using the position space approach at { P z , P z } = { , } × (2 π/ large η values may be required to achieve percent-levelprecision. Appendix E: Alternate approach in position space
In this appendix an alternate approach to extract theCollins-Soper kernel in b z -space is investigated, as sug-gested in Ref. [34]. By taking the FT of the matchingkernel C ns ( µ, xP z ), the Collins-Soper kernel can be ex-pressed as γ q, FI ζ ( µ, b T ) = 1ln( P z /P z ) × ln (cid:82) db z ¯ C ns ( µ, y − b z P z , P z ) P z B MS γ ( b z ,(cid:126)b T , µ, P z ) (cid:82) db z ¯ C ns ( µ, y − b z P z , P z ) P z B MS γ ( b z ,(cid:126)b T , µ, P z ) , (E1)where¯ C ns ( µ, b z P z , P z ) ≡ (cid:90) dx e i x ( b z P z ) (cid:2) C ns ( µ, xP z ) (cid:3) − , (E2) △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○□ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ ○ △ - - - - △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ ○ △ - - - - FIG. 16: Collins-Soper kernel computed in the position spaceapproach from the lattice data. Upper panel: extraction withForm I (Eq. (E2)); lower panel: extraction with Form II(Eq. (E4)). and the inverse of the matching kernel C ns ( µ, xP z ) isobtained by expanding in α s . An alternative form is γ q, FII ζ ( µ, b T ) = 1ln( P z /P z ) × ln (cid:82) db z ¯ C (cid:48) ns ( µ, y − b z P z , P z ) P z B MS γ ( b z ,(cid:126)b T , µ, P z ) (cid:82) db z ¯ C (cid:48) ns ( µ, y − b z P z , P z ) P z B MS γ ( b z ,(cid:126)b T , µ, P z ) , (E3)where¯ C (cid:48) ns ( µ, b z P z , P z ) ≡ (cid:90) dx e i x ( b z P z ) C ns ( µ, xP z ) . (E4)Eqs. (E1) and (E3) are denoted as Form I and FormII, respectively; studying both forms enables a consis-tency check. Similar to the results obtained via the DFTmethod outlined in App. D, the Collins-Soper kernel ob-tained by either Form I or Form II should not depend onthe value of y or on the momentum pair { P z , P z } , whichprovides another handle on the relevant systematic un-certainties.Collins-Soper kernels extracted with the position-spaceapproach for the toy model of Eq. (D1) are shown inFig. 15. The two forms do not yield consistent answers,which indicates that they are not numerically equivalent1and that the perturbative convergence is lost in either orboth of the convolution integrals. Nevertheless, it is clearthat with Form I the extracted Collins-Soper kernel doesnot stabilize to the correct result with increasing η , whilewith Form II the ratio stabilizes around the true value forsufficiently large η . With sufficiently large η it is possiblethat this approach may provide a reliable determinationof the Collins-Soper kernel, although this will need to beinvestigated carefully in future work.The results of applying the position space approach tothe lattice QCD results in this study are shown in Fig. 16,which is compared to the Collins-Soper kernel extractedfrom the fits with the Bernstein polynomial model, dis- cussed in the main text. With the range of η values com-puted in the numerical study there is no plateau in the y -space Collins-Soper kernels, and the different choices ofmomentum pairs do not appear to converge. Althoughthe result extracted at { P z , P z } = { , } × π/L is con-sistent with the results extracted using the model fitsat the minimum value, this consistency is not found atdifferent b T values. Further investigation is needed toconfirm whether the position-space approach via Form Ior II can provide a valuable consistency check against theDFT approach with larger physical lattice volumes usedin calculations. [1] J. C. Collins and D. E. Soper, Nucl. Phys. B193 , 381(1981), [Erratum: Nucl. Phys.B213,545(1983)].[2] J. C. Collins and D. E. Soper, Nucl. Phys.
B197 , 446(1982).[3] J. C. Collins, D. E. Soper, and G. F. Sterman, Nucl.Phys.
B250 , 199 (1985).[4] F. Landry, R. Brock, G. Ladinsky, and C. P. Yuan, Phys.Rev.
D63 , 013004 (2001), hep-ph/9905391.[5] F. Landry, R. Brock, P. M. Nadolsky, and C. P. Yuan,Phys. Rev.
D67 , 073016 (2003), hep-ph/0212159.[6] A. V. Konychev and P. M. Nadolsky, Phys. Lett.
B633 ,710 (2006), hep-ph/0506225.[7] P. Sun, J. Isaacson, C. P. Yuan, and F. Yuan, Int. J.Mod. Phys.
A33 , 1841006 (2018), 1406.3073.[8] U. D’Alesio, M. G. Echevarria, S. Melis, and I. Scimemi,JHEP , 098 (2014), 1407.3311.[9] M. G. Echevarria, A. Idilbi, Z.-B. Kang, and I. Vitev,Phys. Rev. D89 , 074013 (2014), 1401.5078.[10] Z.-B. Kang, A. Prokudin, P. Sun, and F. Yuan, Phys.Rev.
D93 , 014009 (2016), 1505.05589.[11] A. Bacchetta, F. Delcarro, C. Pisano, M. Radici,and A. Signori, JHEP , 081 (2017), [Erratum:JHEP06,051(2019)], 1703.10157.[12] I. Scimemi and A. Vladimirov, Eur. Phys. J. C78 , 89(2018), 1706.01473.[13] V. Bertone, I. Scimemi, and A. Vladimirov, JHEP ,028 (2019), 1902.08474.[14] I. Scimemi and A. Vladimirov, JHEP , 137 (2020),1912.06532.[15] A. Bacchetta, V. Bertone, C. Bissolotti, G. Bozzi, F. Del-carro, F. Piacenza, and M. Radici (2019), 1912.07550.[16] F. Gautheron et al. (COMPASS) (2010).[17] J. Dudek et al., Eur. Phys. J. A48 , 187 (2012),1208.1244.[18] E.-C. Aschenauer et al. (2015), 1501.01220.[19] A. Accardi et al., Eur. Phys. J.
A52 , 268 (2016),1212.1701.[20] M. G. A. Buffing, Z.-B. Kang, K. Lee, and X. Liu (2018),1812.07549.[21] D. Gutierrez-Reyes, I. Scimemi, W. J. Waalewijn, andL. Zoppi, JHEP , 031 (2019), 1904.04259.[22] A. A. Vladimirov (2020), 2003.02288.[23] B. U. Musch, P. Hagler, J. W. Negele, and A. Schafer,Phys. Rev. D83 , 094507 (2011), 1011.1213.[24] B. U. Musch, P. Hagler, M. Engelhardt, J. W. Negele, and A. Schafer, Phys. Rev.
D85 , 094510 (2012),1111.4249.[25] M. Engelhardt, P. H¨agler, B. Musch, J. Negele, andA. Sch¨afer, Phys. Rev.
D93 , 054501 (2016), 1506.07826.[26] B. Yoon, T. Bhattacharya, M. Engelhardt, J. Green,R. Gupta, P. H¨agler, B. Musch, J. Negele, A. Pochin-sky, and S. Syritsyn, in
Proceedings, 33rd Interna-tional Symposium on Lattice Field Theory (Lattice 2015):Kobe, Japan, July 14-18, 2015 , SISSA (SISSA, 2015),1601.05717.[27] B. Yoon, M. Engelhardt, R. Gupta, T. Bhattacharya,J. R. Green, B. U. Musch, J. W. Negele, A. V. Pochinsky,A. Sch¨afer, and S. N. Syritsyn, Phys. Rev.
D96 , 094508(2017), 1706.03406.[28] X. Ji, Phys. Rev. Lett. , 262002 (2013), 1305.1539.[29] X. Ji, Sci. China Phys. Mech. Astron. , 1407 (2014),1404.6680.[30] X. Ji, P. Sun, X. Xiong, and F. Yuan, Phys. Rev. D91 ,074009 (2015), 1405.7640.[31] X. Ji, L.-C. Jin, F. Yuan, J.-H. Zhang, and Y. Zhao,Phys. Rev.
D99 , 114006 (2019), 1801.05930.[32] M. A. Ebert, I. W. Stewart, and Y. Zhao, Phys. Rev.
D99 , 034505 (2019), 1811.00026.[33] M. A. Ebert, I. W. Stewart, and Y. Zhao, JHEP , 037(2019), 1901.03685.[34] M. A. Ebert, I. W. Stewart, and Y. Zhao, JHEP , 099(2020), 1910.08569.[35] X. Ji, Y. Liu, and Y.-S. Liu, Nucl. Phys. B955 , 115054(2020), 1910.11415.[36] X. Ji, Y. Liu, and Y.-S. Liu (2019), 1911.03840.[37] A. A. Vladimirov and A. Sch¨afer, Phys. Rev.
D101 ,074517 (2020), 2002.07527.[38] M. Constantinou, H. Panagopoulos, and G. Spanoudes,Phys. Rev.
D99 , 074508 (2019), 1901.03862.[39] P. Shanahan, M. L. Wagman, and Y. Zhao, Phys. Rev.
D101 , 074505 (2020), 1911.00800.[40] J. R. Green, K. Jansen, and F. Steffens, Phys. Rev.
D101 , 074509 (2020), 2002.09408.[41] W. Detmold and M. G. Endres, Phys. Rev.
D97 , 074507(2018), 1801.06132.[42] M. G. Endres, R. C. Brower, W. Detmold, K. Orginos,and A. V. Pochinsky, Phys. Rev.
D92 , 114516 (2015),1510.04675.[43] M. Asakawa, T. Hatsuda, T. Iritani, E. Itou, M. Ki-tazawa, and H. Suzuki (2015), 1503.06516. [44] M. L¨uscher, JHEP , 071 (2010), [Erratum:JHEP03,092(2014)], 1006.4518.[45] B. Sheikholeslami and R. Wohlert, Nucl. Phys. B259 ,572 (1985).[46] G. S. Bali, B. Lang, B. U. Musch, and A. Sch¨afer, Phys.Rev.
D93 , 094515 (2016), 1602.05525.[47] Y. Li and H. X. Zhu, Phys. Rev. Lett. , 022004(2017), 1604.01404.[48] A. A. Vladimirov, Phys. Rev. Lett. , 062001 (2017),1610.05791.[49] R. A. Brice˜no, J. V. Guerrero, M. T. Hansen, and C. J.Monahan, Phys. Rev.
D98 , 014511 (2018), 1805.01034.[50] L. Piegl and W. Tiller,
The NURBS Book (Springer-Verlag, New York, NY, USA, 1996), 2nd ed.[51] S. Bethke, Eur. Phys. J.
C64 , 689 (2009), 0908.1135.[52] T. van Ritbergen, J. A. M. Vermaseren, and S. A. Larin,Phys. Lett.
B400 , 379 (1997), hep-ph/9701390.[53] P. Sun and F. Yuan, Phys. Rev.
D88 , 034016 (2013),1304.5037.[54] P. Sun and F. Yuan, Phys. Rev.
D88 , 114012 (2013),1308.5003.[55] J. Collins and T. Rogers, Phys. Rev.
D91 , 074020 (2015),1412.3820.[56] A. Pochinsky, Qlua. https://usqcd.lns.mit.edu/qlua.[57] R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD),Nucl. Phys. Proc. Suppl. , 832 (2005), hep-lat/0409003.[58] S. R. Beane et al. (2020), 2003.12130.[59] H. Akaike, IEEE Transactions on Automatic Control ,716 (1974), ISSN 2334-3303.[60] O. Ledoit and M. Wolf, Journal of MultivariateAnalysis , 365 (2004), ISSN 0047-259X, URL . [61] E. Rinaldi, S. Syritsyn, M. L. Wagman, M. I. Buchoff,C. Schroeder, and J. Wasem, Phys. Rev. D99 , 074510(2019), 1901.07519.[62] G. Golub and V. Pereyra, Inverse Problems , R1(2003), URL http://stacks.iop.org/0266-5611/19/i=2/a=201 .[63] D. P. O’Leary and B. W. Rust, Computational Optimiza-tion and Applications , 579 (2013).[64] A. C. Davison and D. V. Hinkley, The Basic Bootstraps (Cambridge University Press, 1997), p. 11–69, Cam-bridge Series in Statistical and Probabilistic Mathemat-ics.[65] P. Young,
Everything you wanted to know about dataanalysis and fitting but were afraid to ask (2012),1210.3781.[66] S. Aoki et al. (Flavour Lattice Averaging Group), Eur.Phys. J.
C80 , 113 (2020), 1902.08191.[67] M. Constantinou and H. Panagopoulos, Phys. Rev.
D96 ,054506 (2017), 1705.11193.[68] J.-W. Chen, T. Ishikawa, L. Jin, H.-W. Lin, Y.-B. Yang,J.-H. Zhang, and Y. Zhao, Phys. Rev.
D97 , 014505(2018), 1706.01295.[69] T. Izubuchi, L. Jin, C. Kallidonis, N. Karthik, S. Mukher-jee, P. Petreczky, C. Shugert, and S. Syritsyn, Phys. Rev.