From Ji to Jaffe-Manohar orbital angular momentum in Lattice QCD using a direct derivative method
M. Engelhardt, J. R. Green, N. Hasan, S. Krieg, S. Meinel, J. Negele, A. Pochinsky, S. Syritsyn
CCERN-TH-2020-133
From Ji to Jaffe-Manohar orbital angular momentum in Lattice QCD using a directderivative method
M. Engelhardt, ∗ J. R. Green, N. Hasan, S. Krieg,
3, 4
S. Meinel, J. Negele, A. Pochinsky, and S. Syritsyn
7, 8 Department of Physics, New Mexico State University, Las Cruces, NM 88003, USA Theoretical Physics Department, CERN, 1211 Geneva 23, Switzerland Bergische Universit¨at Wuppertal, 42119 Wuppertal, Germany IAS, J¨ulich Supercomputing Centre, Forschungszentrum J¨ulich, 52425 J¨ulich, Germany Department of Physics, University of Arizona, Tucson, AZ 85721, USA Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA
A Lattice QCD approach to quark orbital angular momentum in the proton based on generalizedtransverse momentum-dependent parton distributions (GTMDs) is enhanced methodologically byincorporating a direct derivative technique. This improvement removes a significant numerical biasthat had been seen to afflict results of a previous study. In particular, the value obtained for Jiquark orbital angular momentum is reconciled with the one obtained independently via Ji’s sumrule, validating the GMTD approach. Since GTMDs simultaneously contain information about thequark impact parameter and transverse momentum, they permit a direct evaluation of the crossproduct of the latter. They are defined through proton matrix elements of a quark bilocal operatorcontaining a Wilson line; the choice in Wilson line path allows one to continuously interpolate fromJi to Jaffe-Manohar quark orbital angular momentum. The latter is seen to be significantly enhancedin magnitude compared to Ji quark orbital angular momentum, confirming previous results.
I. INTRODUCTION
The manner in which the spin of the proton arises from the spins and orbital angular momenta of its quark andgluon constituents has been the object of sustained study. Efforts to resolve this so-called proton spin puzzle weresparked by the finding, in EMC experiments [1, 2], that the quark spins alone fail to provide a satisfactory account ofthe proton’s overall spin. Naturally, also the methods of Lattice QCD have been brought to bear on the problem, withthe standard calculational scheme relying on Ji’s sum rule [3]. The sum rule relates the total quark angular momentum J to specific moments of generalized parton distributions (GPDs), and by combining this with a calculation of thequark spin S [4, 5], one can then also isolate the quark orbital angular momentum L = J − S [6–12]. The morerecent studies have furthermore begun to gain control over the gluon angular momentum [10–12] and gluon spin [13]contributions.The aforementioned indirect approach to quark orbital angular momentum is limited specifically to the Ji decom-position of proton spin associated with Ji’s sum rule. However, the definition of quark orbital angular momentum inQCD is not unique, since the matter degrees of freedom in a gauge theory cannot be unambiguously separated fromthe gauge degrees of freedom. Quarks necessarily carry gauge fields along with them, and it is a matter of definitionto what extent these are included in the evaluation of quark orbital angular momentum. In addition to the Ji decom-position of proton spin, another widely studied decomposition scheme is the one due to Jaffe and Manohar [14]. Itpossesses the conceptual advantage of allowing for a partonic interpretation of the angular momentum distributions.A formulation that offers a direct path to evaluating quark orbital angular momentum and that encompasses bothof the aforementioned decompositions is the one in terms of generalized transverse momentum-dependent partondistributions (GTMDs) [15–17]. GTMDs, as functions of quark transverse momentum k T as well as momentumtransfer ∆ T , are related, through the Fourier conjugate pair (∆ T , b T ) (with b T denoting the quark impact parameter),to Wigner functions that allow one to sample directly orbital angular momentum b T × k T . GTMDs are defined througha quark bilocal operator containing a Wilson line, and it is the choice of path of the Wilson line that allows one to ∗ [email protected] a r X i v : . [ h e p - l a t ] A ug access different definitions of quark orbital angular momentum. In particular, it was realized in [18] that Jaffe-Manoharorbital angular momentum is associated with a staple-shaped path in the limit of infinite staple length, an operatortype extensively studied in the context of standard transverse momentum-dependent parton distributions (TMDs).On the other hand, Ji orbital angular momentum results from using a straight path [19–22].An initial Lattice QCD exploration of the GTMD approach to quark orbital angular momentum was undertakenin [23]. By varying the staple length of a staple-shaped Wilson line path in small steps, a quasi-continuous, gauge-invariant interpolation between the Ji and Jaffe-Manohar limits was realized. In performing this study, it was possibleto take recourse to concepts and methods from standard Lattice TMD studies [24–28], since, as noted above, thesame operator structure enters; GTMD matrix elements only differ by their additional dependence on the momentumtransfer ∆ T . Jaffe-Manohar quark orbital angular momentum was seen to be significantly enhanced in magnituderelative to its Ji counterpart. The results obtained in [23] were, however, affected by one substantial shortcoming:Although the relative comparison between Ji and Jaffe-Manohar quark orbital angular momentum could be expectedto be trustworthy, in absolute terms, the orbital angular momenta were significantly underestimated for a technicalreason. Namely, the weighting by b T in orbital angular momentum b T × k T corresponds to computing a derivativewith respect to ∆ T of the relevant GTMD matrix element. This derivative was realized via a finite difference over amomentum interval that was much too large to yield an accurate estimate. This became directly apparent in comparingthe result for Ji quark orbital angular momentum with the corresponding result obtained using Ji’s sum rule. Thediscrepancy amounted to approximately a factor 2. The present work resolves this discrepancy by incorporating adirect derivative method to evaluate the ∆ T -derivative. This methodological improvement removes the described biasby construction, and will be seen to reconcile the results obtained through the GTMD approach and through Ji’s sumrule. This validates the GTMD approach as implemented here.The present study also uses pion mass m π = 317 MeV, which is significantly lower than the mass m π = 518 MeVused in the initial exploration [23]. II. GENERALIZED TRANSVERSE MOMENTUM-DEPENDENT (GTMD) APPROACH TO QUARKORBITAL ANGULAR MOMENTUM
As laid out in detail in [23], cf. also [16, 17], the longitudinal quark orbital angular momentum component L ina longitudinally polarized proton propagating in the 3-direction with momentum P can be evaluated within LatticeQCD in units of the number of valence quarks n via L n = 1 a (cid:15) ij ∂∂ ∆ T,j (Φ( a(cid:126)e i ) − Φ( − a(cid:126)e i ))Φ( a(cid:126)e i ) + Φ( − a(cid:126)e i ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∆ T =0 (1)(summation over i, j implied), with the proton matrix elementΦ( z T ) = (cid:104) P + ∆ T / , S = (cid:126)e | ψ ( − z T / γ + U [ − z T / , z T / ψ ( z T / | P − ∆ T / , S = (cid:126)e (cid:105) . (2)Here, (cid:126)e denotes the unit vector in longitudinal direction, whereas (cid:126)e i is a unit vector in a transverse direction; a denotesthe lattice spacing. The momentum transfer ∆ T and the operator separation z T are purely transverse and orthogonalto each other. U is a Wilson line connecting the quark operators ψ , ψ ; its path remains to be specified. This structurecan be understood as follows: In the limit z T →
0, the operator in Φ( z T ) reduces to the light-cone +-componentof the vector current, and therefore, at ∆ T = 0, Φ( z T ) simply counts valence quarks (up to a normalization factor2 P + ). This motivates the denominator in (1). Note the use of the nonlocal current with separation a , matching thenumerator; this will be revisited below. Consider now taking the ∆ T -derivative of Φ( z T ), and evaluating at ∆ T = 0.The momentum transfer ∆ T is Fourier conjugate to the quark impact parameter b T , and therefore, this operationamounts to weighting the counting of quarks by their impact parameter b T . Likewise, the operator separation z T is Fourier conjugate to the quark transverse momentum k T . Thus, taking the derivative with respect to z T andevaluating at z T = 0 amounts to weighting the counting of quarks by their transverse momentum k T . The numeratorin (1), together with the division by a , is the appropriate finite-difference realization of such a z T -derivative; atfinite lattice cutoff a , distances smaller than a cannot be resolved. In aggregate, therefore, (1) evaluates the total L = b T × k T of the quarks in the proton, normalized to the number of valence quarks n .An important role in the evaluation of quark orbital angular momentum falls to the Wilson line U in (2). In agauge theory, the matter degrees of freedom cannot be treated in complete isolation; they necessarily carry gaugefields with them, and the evaluation of quark orbital angular momentum depends on a definition of what part of theoverall gauge field to apportion to quarks as one decomposes orbital angular momentum into a quark and a gluonpart. The path of the Wilson line U carries this information. In the following, staple-shaped paths U ≡ U [ − z/ , ηv − z/ , ηv + z/ , z/
2] (3)will be considered, in which the points listed in the argument of U are connected by straight Wilson lines, as illustratedin Fig. 1. The direction of the staple legs is defined by the vector v , with their length scaled by the parameter η . For η = 0, the path reduces to a straight line between − z/ z/
2. For ease of notation, in the following, η will alsobe allowed to be negative to reverse the direction of the staple (keeping v fixed). This class of gauge links containstwo important limits: η = 0 corresponds to the Ji decomposition of proton spin [19–22], whereas η → ±∞ , with v pointing in a light-like direction, corresponds to the Jaffe-Manohar decomposition of proton spin [18, 20, 22]. Byvarying η continuously, a gauge-invariant interpolation between these two decompositions can be obtained. (cid:1) η → ∞ z ηv + z − z ηv − z v FIG. 1. Path of the gauge connection U , cf. (3), in the correlator (2). The operator in (2), extracting information about quark momentum from the proton state, is of the standard formused in the definition of transverse momentum-dependent parton distributions (TMDs) [29–31]. The matrix element(2) only differs from the standard TMD correlator by the introduction of the nonvanishing momentum transfer ∆ T inthe external states, defining a generalized TMD (GTMD) correlator [15] in which the quark momentum informationis supplemented by quark impact parameter information. Consequently, considerations from the standard TMDframework can be applied [32] in treating the TMD operator in (2). Physically, the staple-shaped gauge link pathincorporates the effect of final state interactions on the struck quark in a semi-inclusive deep inelastic scattering(SIDIS) process. The staple legs represent semi-classical quark paths along which gluon exchanges with the protonremnant are summed . Generalized to the impact parameter-dependent case, the staple-shaped gauge link path thusincorporates into the quark orbital angular momentum the torque accumulated by the struck quark as it is leavingthe proton [20]. For η = 0, this torque vanishes.From the formal point of view, the TMD operator contains divergences that are absorbed into renormalization andsoft factors in the standard TMD framework; these factors appear multiplicatively in the continuum theory [29, 31–33].They are identical for all four instances of Φ in (1) and the ratio is therefore designed to cancel them (of course, thesefactors do not depend on ∆ T , which only enters through the external states). This is the chief motivation for formingthe ratio (1) and employing the nonlocal operator with separation a in the denominator; in this way, operators inthe numerator and denominator match already at finite lattice spacing. Setting the number of valence quarks n tothe appropriate integer serves as the renormalization condition. Nonetheless, since the operator separations in (1)are proportional to the lattice spacing, additional ultraviolet divergences arise as the lattice spacing goes to zero.This is equivalent to the observation in momentum space that, even if one has constructed a renormalized TMD, k T -moments thereof may still diverge as one lets z T , which acts as a regulator on k T -integrals, go to zero. The specificscheme to control that divergence adopted in (1) amounts to identifying the transverse momentum cutoff z T with thelattice resolution a . To connect (1) with its counterpart in other renormalization schemes such as the standard M S scheme, an additional matching would be required that is not available at present. The numerical results presentedbelow suggest that this unquantified systematic uncertainty is minor. The scale evolution of quark orbital angularmomentum has been discussed in detail recently in [34], giving an estimate of the variation expected in the regime inwhich the lattice calculation is performed. The numerical calculation to follow is carried out at a single lattice spacing a ; it would be interesting to extend it to several lattice spacings to directly observe the scale evolution of the results.It should be noted that the multiplicative nature of the renormalization factors in the continuum theory does notstraightforwardly extend to the lattice theory. In general, the breaking of chiral symmetry engendered by the Wilsonfermion discretization used in this work generates operator mixing within the class of TMD operators [28, 35–37] thatprecludes a simple factoring out of renormalization factors and cancellation in the ratio (1). Also these effects willnot be studied quantitatively in the present work, and they constitute a further systematic uncertainty. A study ofthe Sivers shift ratio [28], in which the same TMD operator appears as in (2), revealed no significant operator mixingeffects at the level of statistical accuracy achieved in the calculation; that study included the gauge ensemble used To be specific, this interpretation pertains to the forward in time, η ≥ v points in thedirection opposite to the proton momentum. Note that the quark orbital angular momentum is an even function of η . also in the present work. This suggests that operator mixing effects also do not introduce a dominant systematic biasin the results obtained here.A standard way to regulate the rapidity divergences [38] contained in the TMD operator for light-like staple direction v is to take v off the light cone into the spacelike region [29, 31], while maintaining a zero transverse component, v T = 0. A convenient Lorentz-invariant way to characterize the direction of v is the Collins-Soper type parameterˆ ζ = v · P (cid:112) | v |√ P (4)on which the numerical results obtained below will also depend. Ultimately, one is interested in their large-ˆ ζ behavior;this corresponds to v approaching the light cone. Choosing v to be spacelike simultaneously facilitates a straightforwardconnection to the standard Lattice QCD methodology for evaluating hadronic matrix elements: Given that thetemporal dimension in Lattice QCD is Euclidean, serving to project out the hadronic ground state, the operatorsof which one evaluates matrix elements cannot be extended in physical time. However, once a spacelike vector v isadopted, the problem at hand can be boosted to a Lorentz frame in which v is purely spatial, and thus the entire TMDoperator in (2) exists at a single time. The lattice calculation can be performed in that frame. Maintaining v T = 0in the lattice frame, v will point purely in the (negative) 3-direction, v ≡ − (cid:126)e . In that case, ˆ ζ = P /m (where m isthe proton mass). Consequently, achieving large ˆ ζ requires large proton momentum component P ; this representsa significant challenge, cf. [27] for a study in the context of the Boer-Mulders effect. On the other hand, in thespecial case η = 0, the dependence on the staple direction v disappears; Ji quark orbital angular momentum is boost-invariant. The corresponding results obtained below will indeed be seen to be independent of P (or, equivalently, ˆ ζ ,if one formally maintains v ≡ − (cid:126)e in the lattice frame).On a lattice of finite extent with periodic boundary conditions in the spatial directions, momenta are quantized.Therefore, performing direct calculations only of the matrix element Φ( z T ), cf. (2), limits the accuracy in determiningits derivative with respect to ∆ T , by forcing one to evaluate finite differences over, in practice, substantial momentumincrements. In the initial study [23] of quark orbital momentum in the proton employing the GTMD approach laidout here, this was the dominant source of systematic uncertainty. It introduced a bias in the overall magnitude of thenumerical data approaching a factor 2. Although relative comparisons performed in [23], such as the one between Jiand Jaffe-Manohar orbital angular momentum, can be expected to be robust with respect to this bias, in absoluteterms, the Ji orbital angular momentum extracted in [23] displayed a significant discrepancy compared to the valueobtained independently via Ji’s sum rule. To remedy this dominant systematic bias, and thereby achieve agreementwith the result obtained via Ji’s sum rule within the statistical fluctuations, is the principal objective and advance ofthe present work. Other systematic uncertainties, such as the ones associated with renormalization discussed furtherabove, are, in comparison, minor, and are accordingly deferred to future work.To completely remove the systematic bias originating from finite difference evaluations of the ∆ T -derivative ofΦ( z T ), a direct derivative method is adopted in the present work. The detailed implementation of the method will bedescribed in the next section, following [39], where the method was first laid out in detail. Heuristically, the method isbased on the observation that arbitrarily small increments in overall momenta can be achieved by twisting the spatialboundary conditions of the quarks, or, equivalently, coupling the quarks to a constant U (1) background gauge field.For the purpose of computing a momentum derivative, this gauge field can be infinitesimally small, and can thereforebe treated perturbatively. In effect, this generates an additional vector current insertion in the diagram one wouldcalculate to obtain Φ( z T ). Thus, by instead directly performing a lattice calculation of the diagram containing theadditional vector current insertion, one directly accesses the ∆ T -derivative of Φ( z T ), excluding any systematic bias.A moderate price one pays is that the additional operator insertion will tend to somewhat increase the statisticalfluctuations in the calculation. III. LATTICE METHODOLOGY
To obtain the proton matrix element Φ( z T ), cf. (2), one calculates three-point functions C [ ˆ O ] together withtwo-point functions C , which are projected onto a definite proton momentum p (cid:48) = P + ∆ T / ∆ T at the operator insertion in C [ ˆ O ], C [ ˆ O ]( t, t f , p (cid:48) , p ) = (cid:88) x f , y e − i x f · p (cid:48) + i y · ( p (cid:48) − p ) tr[Γ pol (cid:104) n ( t f , x f ) ˆ O ( t, y )¯ n (0 , (cid:105) ] (5) C ( t f , p (cid:48) ) = (cid:88) x f e − i x f · p (cid:48) tr[Γ pol (cid:104) n ( t f , x f )¯ n (0 , (cid:105) ] . (6)By momentum conservation, a definite source momentum p = P − ∆ T / C [ ˆ O ]. The protoninterpolating fields n ( t, x ) are constructed in practice using Wuppertal-smeared quarks, to be discussed in more detailbelow, such as to optimize overlap with the true proton state. The projector Γ pol = (1+ γ ) (1 − iγ γ ) selects statespolarized in the 3-direction. As already discussed further above, the lattice calculation is performed in a Lorentz framein which the TMD operator ˆ O , specified in (2), exists at the single time t ; in particular, v = − (cid:126)e , corresponding toˆ ζ = P /m . In general, the three-point function C [ ˆ O ] contains both connected contributions, in which ˆ O is insertedinto a valence quark propagator, as well as disconnected contributions, in which ˆ O is inserted into a sea quark loop.In the present investigation, only the former are taken into account. The latter contributions, which are associatedwith significantly higher computational cost, are expected to be minor, and are excluded. In the isovector u − d quark channel, the disconnected contributions cancel exactly; in that case, no systematic uncertainty is associatedwith neglecting the disconnected diagrams.Having calculated the three-point and two-point functions, the matrix element (2) is obtained from the ratio2 E ( p (cid:48) ) C [ ˆ O ]( t, t f , p (cid:48) , p ) C ( t f , p (cid:48) ) −→ Φ( z T ) (7)which exhibits plateaus in t for 0 (cid:28) t (cid:28) t f , yielding Φ( z T ). Here, E ( p ) = E ( p (cid:48) ) denotes the energy of the initial andfinal proton states. Note that, for general choices of initial and final momenta, the ratio (7) has to be replaced bya more general expression [6, 7]; it is the specific symmetric choice of the initial and final momenta p = P − ∆ T / p (cid:48) = P + ∆ T / t f . In the present study, data for only one source-sink separation t f = 1 .
14 fm weregathered, and therefore it will not be possible to estimate excited state effects quantitatively. Previous form factorstudies on the same gauge ensemble [40, 41] showed that the importance of excited state effects depends substantiallyon the specific observable considered. In some cases, the systematic bias at the separation t f = 1 .
14 fm was seen tobe smaller than the statistical fluctuations, in other cases, it exceeded the latter by factors up to 2 to 3. Controllingfor excited state effects will certainly be desirable in future investigations.Consider now evaluating the ∆ T -derivative of Φ( z T ) at ∆ T = 0, as called for by (1). First, note that the two-pointfunction C is an even function of ∆ T ; therefore, only the derivative of the three-point function C [ ˆ O ] is needed ,while one can directly use C ( t f , P ) in (7). To proceed, it is useful to write C [ ˆ O ] explicitly in terms of theappropriate propagators, combining the coordinates for ease of notation into four-vectors, e.g., ( t, y ) ≡ y , C [ ˆ O ]( t, t f , p (cid:48) , p ) = (cid:88) x f , y e − i x f · P + i ( y − x f / · ∆ T · (8) · (cid:68) tr (cid:104) S n ¯ n Γ pol (0; x f ) G sm-pt ( x f , y − z T / γ + U ( y − z T / , y + z T / G pt-sm ( y + z T / , (cid:105)(cid:69) = (cid:88) x f , y e − i x f · P e − i ( x f − ( y − z T / · ∆ T / e i ( y + z T / · ∆ T / · (9) · (cid:28) tr (cid:20)(cid:16) γ G pt-sm ( y − z T / , x f ) γ S n ¯ n † Γ pol (0; x f ) (cid:17) † γ + U ( y − z T / , y + z T / G pt-sm ( y + z T / , (cid:21)(cid:29) = (cid:88) x f , y e − i x f · P · (10) · (cid:28) tr (cid:20)(cid:16) γ G pt-sm ( y − z T / , x f ; ∆ T / γ S n ¯ n † Γ pol (0; x f ) (cid:17) † γ + U ( y − z T / , y + z T / G pt-sm ( y + z T / , − ∆ T / (cid:21)(cid:29) . Here, S n ¯ n Γ pol denotes the standard sequential source formed at the sink time t f , as described in detail, e.g., in theAppendix of [43]; G pt-sm ( r, s ) is the standard quark propagator from a smeared source to a point sink. In the secondstep, the γ -hermiticity of the propagators was used, and in the third step, the twisted propagator G pt-sm ( r, s ; q ) = e − i ( r − s ) · q G pt-sm ( r, s ) (11)was introduced. Thus, the ∆ T -dependence has been entirely absorbed into the twisted propagators, and it is theirderivatives one needs in order to obtain the ∆ T -derivative of C [ ˆ O ]. Following [39, 42], in the absence of smearing, In applications requiring higher derivatives with respect to momentum transfer, such as calculations of charge radii [42], also derivativesof the two-point function enter. the derivative of the twisted point-to-point propagator can be cast in the form ∂∂q j G pt-pt ( r, s ; q ) (cid:12)(cid:12)(cid:12)(cid:12) q =0 = − i (cid:88) z G pt-pt ( r, z )Γ jV G pt-pt ( z, s ) (12)where the sum extends over the four-dimensional coordinate z , and the vector current insertion acts asΓ jV G pt-pt ( z, s ) = U † j ( z − (cid:126)e j ) 1 + γ j G pt-pt ( z − (cid:126)e j , s ) − U j ( z ) 1 − γ j G pt-pt ( z + (cid:126)e j , s ) . (13)Supplementing the point-to-point propagator with a smearing kernel, G pt-sm ( r, s ; q ) = e − i ( r − s ) · q (cid:88) u G pt-pt ( r, u ) K ( u, s ) = (cid:88) u G pt-pt ( r, u ; q ) K ( u, s ; q ) (14)where also the twisted smearing kernel K ( u, s ; q ) = e − i ( u − s ) · q K ( u, s ) (15)has been introduced, one has the derivative ∂∂q j G pt-sm ( r, s ; q ) (cid:12)(cid:12)(cid:12)(cid:12) q =0 = (cid:88) z G pt-pt ( r, z ) (cid:34) − i (cid:88) u Γ jV G pt-pt ( z, u ) K ( u, s ) + ∂∂q j K ( z, s ; q ) (cid:12)(cid:12)(cid:12)(cid:12) q =0 (cid:35) (16)where the derivative of the twisted smearing kernel will be treated below. Thus, in order to calculate the ∆ T -derivativeof C [ ˆ O ], one has to evaluate two additional propagators compared to a standard calculation of C [ ˆ O ] itself. Inthe latter case, one needs to evaluate the forward propagator from a smeared source K , and the backward propagatorfrom the smeared sequential source Kγ S n ¯ n † Γ pol ; now, one additionally needs propagators from the sources − i Γ jV G pt-pt K + ∂∂q j K (cid:12)(cid:12)(cid:12)(cid:12) q =0 , − i Γ jV G pt-pt Kγ S n ¯ n † Γ pol + (cid:32) ∂∂q j K (cid:12)(cid:12)(cid:12)(cid:12) q =0 (cid:33) γ S n ¯ n † Γ pol . (17)It remains to construct the derivative of the twisted smearing kernel [42]. A single step of (twisted) Wuppertalsmearing is defined by K ( u, s ; q ) = e − i ( u − s ) · q
11 + 6 α δ u,s + α (cid:88) j =1 (cid:104) U j ( u ) δ u + (cid:126)e j ,s + U † j ( u − (cid:126)e j ) δ u − (cid:126)e j ,s (cid:105) (18)= 11 + 6 α δ u,s + α (cid:88) j =1 (cid:104) e iq j U j ( u ) δ u + (cid:126)e j ,s + e − iq j U † j ( u − (cid:126)e j ) δ u − (cid:126)e j ,s (cid:105) (19)so that its derivative at zero momentum is K (cid:48) ( u, s ) ≡ ∂∂q j K ( u, s ; q ) (cid:12)(cid:12)(cid:12)(cid:12) q =0 = α α (cid:104) iU j ( u ) δ u + (cid:126)e j ,s − iU † j ( u − (cid:126)e j ) δ u − (cid:126)e j ,s (cid:105) . (20)If the smearing kernel K is given by N steps of Wuppertal smearing, K ( u, s ; q ) = (cid:88) w ,w ,...,w N − K ( u, w ; q ) K ( w , w ; q ) . . . K ( w N − , s ; q ) (21)then its derivative at zero momentum can be computed iteratively as K (cid:48) ≡ ( K N ) (cid:48) = K (cid:48) K N − + K ( K N − ) (cid:48) . (22)The numerical data for the present study were generated using a 2+1-flavor isotropic clover fermion ensemble on32 ×
96 lattices generated by R. Edwards, B. Jo´o and K. Orginos with lattice spacing a = 0 .
114 fm and pion mass m π = 317 MeV. The present investigation is therefore also significantly closer to the physical limit than the initialexploration [23]. A total of 23224 data samples was gathered on 968 gauge configurations. The Euclidean temporal (cid:224)(cid:224) (cid:224)(cid:224) (cid:224)(cid:224) (cid:225)(cid:225) (cid:236)(cid:236) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) Ζ(cid:96) L (cid:72) Η (cid:61) (cid:76) (cid:144) n (cid:72) Η (cid:61) (cid:76) u (cid:45) d quarks m Π (cid:61)
317 MeV (cid:165)
FIG. 2. Ji quark orbital angular momentum, i.e., the η = 0 limit, for the three values of ˆ ζ probed, with the average plottedat ˆ ζ = ∞ (open square). The filled diamond represents the value extracted at the same pion mass in the MS scheme at thescale µ = 4 GeV via Ji’s sum rule [7]. The isovector u − d quark combination was evaluated. The shown uncertainties arestatistical jackknife errors. separation between proton sources and sinks was t f = 10 a . HYP-smearing was applied to the lattice links used inconstructing the Wilson line U in (2). This leads to the renormalization and soft factors associated with the operatorin (2) corresponding more closely to their tree-level values even before their cancellation in the ratio (1). Three spatialproton momenta were used, P · L/ (2 π ) = (0 , , n P ) with n P = 0 , ,
2, where L = 32 a denotes the spatial lattice extent.Although the corresponding range of the Collins-Soper parameter ˆ ζ = 0 , . , .
63 is limited, and one cannot a prioriexpect to obtain a good indication of the large-ˆ ζ behavior with data in this range, it will be seen below that the resultsfor Jaffe-Manohar orbital angular momentum already appear to stabilize in the region of the two nonzero values ofˆ ζ . Corroboration concerning this suggested early onset of asymptotic behavior by future studies including larger ˆ ζ would certainly be desirable. IV. NUMERICAL RESULTS
Considering initially the special case of a straight Wilson line, η = 0, corresponding to Ji quark orbital angularmomentum, Fig. 2 displays the results obtained in the isovector case at the three available values of ˆ ζ . Recall thatˆ ζ is defined here through v = − (cid:126)e in the lattice frame even when η = 0; with that definition, on the other hand, Jiquark orbital angular momentum should then be independent of ˆ ζ , since v does not in fact enter its construction.This is borne out by the data in Fig. 2. The residual apparent trend in the data may be due to the deviation betweenthe lattice dispersion relation and the continuum dispersion relation used in the data analysis, but is also consistentwith statistical fluctuation.Performing a χ fit of a constant in ˆ ζ to the data yields the average value plotted at ˆ ζ = ∞ , taken here tolabel the physical limit. This result is confronted with an independent lattice determination of Ji quark orbitalangular momentum via Ji’s sum rule at the same pion mass, in the M S scheme at the scale µ = 4 GeV [7]. The twodeterminations are in good agreement; the discrepancy observed in the initial exploration [23] is entirely resolved. Thisvalidates the use of the GTMD approach, properly implemented in particular with respect to taking the ∆ T -derivativein (1), to calculating quark orbital angular momentum in the proton. The result corroborates the assumption thatthe various further systematic uncertainties noted above, i.e., stemming from renormalization and matching, operatormixing, or excited state contaminations are minor and do not bias the result beyond the statistical uncertainty.Departing from the η = 0 limit, one probes the torque [20] due to final state interactions accumulated by a quarkstruck in a deep inelastic scattering process as it exits the proton. The η = ±∞ limit corresponds to Jaffe-Manoharorbital angular momentum. By varying η gradually, a gauge-invariant, continuous interpolation between the Ji andJaffe-Manohar limits can be exhibited. This is shown in Fig. 3, again for the isovector u − d quark channel, and forthe three values of ˆ ζ probed. Note that the plots are normalized to the magnitude of the Ji quark orbital angularmomentum, i.e., the result obtained at η = 0.Evidently, the torque supplied by the final state interactions is appreciable, as already observed in [23]. Comparedto the initial Ji value, quark orbital angular momentum is enhanced in magnitude as one proceeds towards theasymptotic Jaffe-Manohar limit. Contrasting the three panels in Fig. 3, the effect strengthens as the Collins-Soper (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:225) (cid:225) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) Η (cid:200) v (cid:200) (cid:144) a (cid:72) L (cid:144) n (cid:76) (cid:144) (cid:200) L (cid:72) Η (cid:61) (cid:76) (cid:144) n (cid:72) Η (cid:61) (cid:76) (cid:200) u (cid:45) d quarks m Π (cid:61)
317 MeV
Ζ(cid:96) (cid:61) (cid:45)(cid:165) (cid:165) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:225) (cid:225) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) Η (cid:200) v (cid:200) (cid:144) a (cid:72) L (cid:144) n (cid:76) (cid:144) (cid:200) L (cid:72) Η (cid:61) (cid:76) (cid:144) n (cid:72) Η (cid:61) (cid:76) (cid:200) u (cid:45) d quarks m Π (cid:61)
317 MeV
Ζ(cid:96) (cid:61) (cid:45)(cid:165) (cid:165) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:225) (cid:225) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) Η (cid:200) v (cid:200) (cid:144) a (cid:72) L (cid:144) n (cid:76) (cid:144) (cid:200) L (cid:72) Η (cid:61) (cid:76) (cid:144) n (cid:72) Η (cid:61) (cid:76) (cid:200) u (cid:45) d quarks m Π (cid:61)
317 MeV
Ζ(cid:96) (cid:61) (cid:45)(cid:165) (cid:165)
FIG. 3. Quark orbital angular momentum as a function of staple length parameter η , normalized to the magnitude of Jiquark orbital angular momentum, i.e., the result obtained at η = 0. This quantity is even under η → − η , corresponding totime reversal. Accordingly, the | η | → ∞ extrapolated values are obtained by averaging the η > η < | η || v | /a = 7 , . . . , u − d quark combination is shown, with the three panelscorresponding to the three available values of ˆ ζ . The shown uncertainties are statistical jackknife errors. parameter ˆ ζ is increased, with Jaffe-Manohar quark orbital angular momentum enhanced by about 30% relative tothe Ji case for the two nonzero values of ˆ ζ . This is somewhat less strong than in the exploration [23]; whether this is agenuine physical trend associated with the change in pion mass from 518 MeV to 317 MeV, or whether it is an artefactof the systematic bias in the calculation in [23] cannot be decided at this point. The fact that the effect strengthenswith rising ˆ ζ suggests that it can be expected to persist into the ˆ ζ → ∞ limit. Fig. 4 displays the integrated torque,i.e., the difference between the Jaffe-Manohar and Ji quark orbital momenta, τ = L ( η = ∞ )3 n ( η = ∞ ) − L ( η =0)3 n ( η =0) , (23)as a function of ˆ ζ , normalized to the magnitude of Ji quark orbital momentum. An extrapolation to the ˆ ζ → ∞ (cid:224)(cid:224) (cid:224)(cid:224) (cid:224)(cid:224) (cid:225)(cid:225) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) Ζ(cid:96) Τ (cid:144) (cid:200) L (cid:72) Η (cid:61) (cid:76) (cid:144) n (cid:72) Η (cid:61) (cid:76) (cid:200) u (cid:45) d quarks m Π (cid:61)
317 MeV (cid:165)
FIG. 4. Integrated torque accumulated by a quark struck in a deep inelastic scattering process along its trajectory exiting theproton, as a function of Collins-Soper parameter ˆ ζ . The data pertain to the isovector u − d quark channel, and are normalizedto the magnitude of the η = 0 Ji orbital angular momentum. An ad hoc extrapolation to the ˆ ζ → ∞ limit is also exhibited.The shown uncertainties are statistical jackknife errors. (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224) (cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:225) (cid:225) (cid:237) (cid:237) (cid:231) (cid:231) (cid:45) (cid:45) (cid:45) (cid:45) Η (cid:200) v (cid:200) (cid:144) a (cid:72) L (cid:144) n (cid:76) (cid:144) (cid:200) L , u (cid:45) d (cid:72) Η (cid:61) (cid:76) (cid:144) n u (cid:45) d (cid:72) Η (cid:61) (cid:76) (cid:200) isoscalar2u quarksd quarks m Π (cid:61)
317 MeV
Ζ(cid:96) (cid:61) (cid:45)(cid:165) (cid:165) (cid:224) (cid:236) (cid:230)
FIG. 5. Flavor-separated quark orbital angular momentum as a function of staple length η , analogous to Fig. 3, at fixedˆ ζ = 0 . d quarks and for two u quarks (i.e., the u -quark data for L /n have been multiplied by 2to compensate for n = 2), as well as for the isoscalar total quark orbital angular momentum, i.e., the sum of the “ d ” and “2 u ”data. All results are still normalized by the magnitude of isovector Ji orbital angular momentum (thus, at η = 0, the “2 u ” and“ d ” data differ by unity). The shown uncertainties are statistical jackknife errors. limit is also shown. The ad hoc fit ansatz, A + B/ ˆ ζ , is not underpinned by a theoretical argument at this point, butis motivated by the good description it provides of the considerably more detailed data as a function of ˆ ζ availablefor the pion Boer-Mulders TMD ratio [27]. Auxiliary information concerning the expected large-ˆ ζ behavior would bedesirable to aid in sharpening the analysis. The ad hoc extrapolation indeed yields a signal in the ˆ ζ → ∞ limit.Generalizing to the flavor-separated case, it should be kept in mind that the additional disconnected contributionsthat arise compared to the isovector case have not been evaluated. These are, however, expected to be minor atthe pion mass m π = 317 MeV used in this calculation. Fig. 5 shows data analogous to Fig. 3 for one value of ˆ ζ ,exhibiting the behavior of d -quark and u -quark orbital angular momentum separately, as well as the total (isoscalar)quark orbital angular momentum. Here, the u -quark data have been normalized to two quarks, i.e., L /n in the u quark case has been multiplied by 2 to compensate for n = 2 for u quarks; hence the “2 u ” label. The isoscalar resultwas then obtained by simple addition of the “ d ” and “2 u ” data .As observed previously in [23], the strong cancellation of the u - and d -quark orbital angular momenta in the protonthat has long been known for the η = 0 Ji case [6, 7] extends to nonzero η and the Jaffe-Manohar limit. Only asmall negative contribution to the spin of the proton from quark orbital angular momentum remains. The data for This may differ slightly from calculating 3 L ,u + d /n u + d at finite statistics. (cid:224)(cid:224) (cid:224)(cid:224) (cid:224)(cid:224) (cid:236)(cid:236) (cid:236)(cid:236) (cid:236)(cid:236) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:225)(cid:225) (cid:237)(cid:237) (cid:231)(cid:231) (cid:45) (cid:45) Ζ(cid:96) Τ (cid:144) (cid:200) L , u (cid:45) d (cid:72) Η (cid:61) (cid:76) (cid:144) n u (cid:45) d (cid:72) Η (cid:61) (cid:76) (cid:200) m Π (cid:61)
317 MeV (cid:165) (cid:224) (cid:236) (cid:230) isoscalar2u quarksd quarks
FIG. 6. Flavor-separated integrated torque accumulated by a quark struck in a deep inelastic scattering process, as a functionof ˆ ζ , analogous to Fig. 4, together with ad hoc extrapolations to infinite ˆ ζ . Data for d quarks, for two u quarks (i.e., the u -quark data for τ have been multiplied by 2 to compensate for n = 2 in the u -quark case), and their sum, the isoscalar u + d quark combination, are shown. Data are normalized by the magnitude of isovector Ji quark orbital angular momentum, as inprevious figures. For better visibility, isoscalar data are slightly displaced horizontally. The shown uncertainties are statisticaljackknife errors. the flavor-separated integrated torque, cf. (23), are collected in Fig. 6 and extrapolated to the ˆ ζ → ∞ limit. At thepresent level of statistics, scarcely a signal is obtained for the flavor-separated integrated torque in that limit; thedata do appear compatible with the observations made at fixed ˆ ζ = 0 .
315 from Fig. 5.
V. CONCLUSIONS AND OUTLOOK
The chief advance of the present study is the adoption of a direct derivative method [39] in the GTMD approach toevaluating quark orbital angular momentum in the proton in Lattice QCD. The introduction of this method has led to areliable quantitative computation of the needed derivative, with respect to momentum transfer, of the relevant GTMDmatrix element. This is validated by the result obtained specifically for the quark orbital angular momentum definedthrough the Ji decomposition of proton spin; it agrees well with the corresponding result obtained independently inLattice QCD calculations relying on Ji’s sum rule [7]. The discrepancy observed in the initial exploration [23] is thusresolved.The agreement between the quark orbital angular momentum calculated in this work using the GTMD approachand the result from the Ji sum rule suggests that other systematic uncertainties, such as ones associated with excitedstate effects, renormalization and matching, as well as operator mixing are minor and do not rise to the level of thestatistical uncertainties of the present calculation.In the GTMD approach, one directly computes quark orbital angular momentum by weighting the appropriateWigner function (related to GTMDs via Fourier transformation) by b T × k T , where b T is the quark impact parameterand k T the quark transverse momentum [16, 17]. The aforementioned derivative with respect to momentum transfersupplies the weighting by b T , its Fourier conjugate. The information on k T , on the other hand, is supplied through thenonlocal TMD operator used in constructing the relevant proton matrix element (2). The treatment of renormalizationissues thus follows closely the methods used in more widely explored Lattice TMD calculations [24–28]. Ratios ofproton matrix elements are constructed to cancel renormalization and soft factors associated with the TMD operator.In effect, one evaluates quark orbital angular momentum in units of the number of valence quarks. Operator mixingeffects [28, 35–37] can spoil these cancellations, and though they appear not to play a significant role in the presentcalculation, these effects will ultimately have to be brought under control.The advantage of this nonlocal operator-based GTMD approach is that one can extend Lattice QCD calculationsbeyond the Ji decomposition of proton spin and establish a continuous, gauge-invariant interpolation from Ji [3] toJaffe-Manohar [14] quark orbital angular momentum. The corresponding information is contained in the choice ofWilson line path in the TMD operator. A straight Wilson line path yields Ji quark orbital angular momentum [19–21], and a staple-shaped Wilson line path, in the limit of infinite staple length, yields its Jaffe-Manohar counterpart[18, 20]. The difference between the two can be interpreted as the integrated torque accumulated by a quark struckin a deep inelastic scattering process as it exits the proton, through final state interactions [20]. It correspondsto a Qiu-Sterman type correlator [20–22]. The data obtained for this term in the present investigation allow one to1observe the gradual accumulation of torque by the quark until it attains the asymptotic Jaffe-Manohar orbital angularmomentum. The latter is enhanced in magnitude relative to the initial Ji value, with the integrated torque addingabout one third of the magnitude of the Ji orbital angular momentum at the pion mass m π = 317 MeV used in thiscalculation.To further sharpen the analysis of quark orbital angular momentum within the GTMD approach, calculations fora sequence of lattice spacings would be desirable, to allow for a direct study of the scale evolution. Also, a morecomprehensive exploration of the dependence on the Collins-Soper parameter ˆ ζ is warranted, to clarify whether theresults indeed already stabilize at the fairly low values of ˆ ζ employed in this work, as suggested by Fig. 4. Ultimately,the large ˆ ζ behavior determines the physical limit. Finally, of course, also further progress towards the physical pionmass must be made.The improved calculation of quark orbital angular momentum based on GTMDs achieved using the present method-ology opens the way to reliably compute other related quantities, such as quark spin-orbit correlations in the proton[16, 22]. Corresponding calculations are in progress, employing domain wall fermions to curtail the operator mixingeffects that are induced when the fermion discretization breaks chiral symmetry. ACKNOWLEDGMENTS
This work benefited from fruitful discussions with M. Burkardt, W. Detmold, R. Gupta, S. Liuti and C. Lorc´e.The lattice calculations performed in this work relied on code developed by B. Musch, as well as the Chroma [44]and Qlua [45] software suites. R. Edwards, B. Jo´o and K. Orginos provided the clover fermion ensemble, which wasgenerated using resources provided by XSEDE (supported by National Science Foundation Grant No. ACI-1053575).Computations were performed using resources of the National Energy Research Scientific Computing Center (NERSC),a U.S. DOE Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. This work wasfurthermore supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under grantDE-FG02-96ER40965 (M.E.), grant DE-SC-0011090 (J.N.), grant DE-SC0018121 (A.P.), as well as through the TMDTopical Collaboration (M.E., J.N.); and it was also supported by the U.S. Department of Energy, Office of Science,Office of High Energy Physics under Award Number DE-SC0009913 (S.M.). N.H. and S.K. received support fromDeutsche Forschungsgemeinschaft through grant SFB-TRR 55, and S.S. is supported by the U.S. National ScienceFoundation under CAREER Award PHY-1847893 and through the RHIC Physics Fellow Program of the RIKENBNL Research Center. [1] J. Ashman et al. [European Muon Collaboration], Phys. Lett.
B206 , 364 (1988).[2] J. Ashman et al. [European Muon Collaboration], Nucl. Phys.
B328 , 1 (1989).[3] X. Ji, Phys. Rev. Lett. , 610 (1997).[4] H.-W. Lin, R. Gupta, B. Yoon, Y.-C. Jang and T. Bhattacharya, Phys. Rev. D 98 , 094512 (2018).[5] J. Liang, Y.-B. Yang, T. Draper, M. Gong and K.-F. Liu, Phys. Rev.
D 98 , 074505 (2018).[6] P. H¨agler et al. [LHP Collaboration], Phys. Rev.
D 77 , 094502 (2008).[7] J. D. Bratt et al. [LHP Collaboration], Phys. Rev.
D 82 , 094502 (2010).[8] M. G¨ockeler, R. Horsley, D. Pleiter, P. E. L. Rakow, A. Sch¨afer, G. Schierholz and W. Schroers [QCDSF Collaboration],Phys. Rev. Lett. , 042002 (2004).[9] G. Bali, S. Collins, M. G¨ockeler, R. R¨odl, A. Sch¨afer and A. Sternbeck [RQCD Collaboration], Phys. Rev. D 100 , 014507(2019).[10] M. Deka, T. Doi, Y.-B. Yang, B. Chakraborty, S. J. Dong, T. Draper, M. Glatzmaier, M. Gong, H.-W. Lin, K.-F. Liu,D. Mankame, N. Mathur and T. Streuer, Phys. Rev.
D 91 (2015), 014505 (2015).[11] C. Alexandrou, M. Constantinou, K. Hadjiyiannakou, K. Jansen, C. Kallidonis, G. Koutsou, A. Vaquero Avil´es-Casco andC. Wiese, Phys. Rev. Lett. , 142002 (2017).[12] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, H. Panagopoulosand G. Spanoudes, Phys. Rev.
D 101 , 094513 (2020).[13] Y.-B. Yang, R. S. Sufian, A. Alexandru, T. Draper, M. J. Glatzmaier, K.-F. Liu and Y. Zhao, Phys. Rev. Lett. , 102001(2017).[14] R. Jaffe and A. Manohar, Nucl. Phys.
B337 , 509 (1990).[15] S. Meißner, A. Metz and M. Schlegel, JHEP , 056 (2009).[16] C. Lorc´e and B. Pasquini, Phys. Rev.
D 84 , 014015 (2011).[17] E. Leader and C. Lorc´e, Phys. Rept. , 163 (2014).[18] Y. Hatta, Phys. Lett.
B708 , 186 (2012).[19] X. Ji, X. Xiong and F. Yuan, Phys. Rev. Lett. , 152005 (2012). [20] M. Burkardt, Phys. Rev. D 88 , 014014 (2013).[21] A. Rajan, A. Courtoy, M. Engelhardt and S. Liuti, Phys. Rev.
D 94 , 034041 (2016).[22] A. Rajan, M. Engelhardt and S. Liuti, Phys. Rev.
D 98 , 074022 (2018).[23] M. Engelhardt, Phys. Rev.
D 95 , 094505 (2017).[24] P. H¨agler, B. Musch, J. W. Negele and A. Sch¨afer, Europhys. Lett. , 61001 (2009).[25] B. Musch, P. H¨agler, J. W. Negele and A. Sch¨afer, Phys. Rev. D 83 , 094507 (2011).[26] B. Musch, P. H¨agler, M. Engelhardt, J. Negele and A. Sch¨afer, Phys. Rev.
D 85 , 094510 (2012).[27] M. Engelhardt, P. H¨agler, B. Musch, J. Negele and A. Sch¨afer, Phys. Rev.
D 93 , 054501 (2016).[28] B. Yoon, M. Engelhardt, R. Gupta, T. Bhattacharya, J. Green, B. Musch, J. Negele, A. Pochinsky, A. Sch¨afer, andS. Syritsyn, Phys. Rev.
D 96 , 094508 (2017).[29] J. C. Collins,
Foundations of Perturbative QCD (Cambridge University Press, 2011).[30] X. Ji, J.-P. Ma and F. Yuan, Phys. Rev.
D 71 , 034005 (2005).[31] S. M. Aybat and T. Rogers, Phys. Rev.
D 83 , 114042 (2011).[32] M. Echevarria, A. Idilbi, K. Kanazawa, C. Lorc´e, A. Metz, B. Pasquini and M. Schlegel, Phys. Lett.
B759 , 336 (2016).[33] S. M. Aybat, J. C. Collins, J.-W. Qiu and T. C. Rogers, Phys. Rev.
D 85 , 034043 (2012).[34] Y. Hatta and X. Yao, Phys. Lett.
B798 , 134941 (2019).[35] M. Constantinou, H. Panagopoulos and G. Spanoudes, Phys. Rev.
D 99 , 074508 (2019).[36] J. R. Green, K. Jansen and F. Steffens, Phys. Rev.
D 101 , 074509 (2020).[37] P. Shanahan, M. Wagman and Y. Zhao, Phys. Rev.
D 101 , 074505 (2020).[38] J. C. Collins, T. C. Rogers and A. M. Stasto, Phys. Rev.
D 77 , 085009 (2008).[39] G. M. de Divitiis, R. Petronzio and N. Tantalo, Phys. Lett.
B718 , 589 (2012).[40] J. R. Green, S. Meinel, M. Engelhardt, S. Krieg, J. Laeuchli, J. Negele, K. Orginos, A. Pochinsky and S. Syritsyn, Phys.Rev.
D 92 , 031501 (2015).[41] J. R. Green, N. Hasan, S. Meinel, M. Engelhardt, S. Krieg, J. Laeuchli, J. Negele, K. Orginos, A. Pochinsky and S. Syritsyn,Phys. Rev.
D 95 , 114502 (2017).[42] N. Hasan, J. Green, S. Meinel, M. Engelhardt, S. Krieg, J. Negele, A. Pochinsky and S. Syritsyn, Phys. Rev.
D 97 , 034504(2018).[43] D. Dolgov, R. Brower, S. Capitani, P. Dreher, J. Negele, A. Pochinsky, D. Renner, N. Eicker, T. Lippert, K. Schilling,R. Edwards and M. Heller [LHP and SESAM Collaborations], Phys. Rev.
D 66 , 034506 (2002).[44] R. G. Edwards and B. Jo´o [SciDAC Collaboration], Nucl. Phys. Proc. Suppl. , 832 (2005).[45] A. Pochinsky,
Qlua , https://usqcd.lns.mit.edu/qluahttps://usqcd.lns.mit.edu/qlua