Complete flavor decomposition of the spin and momentum fraction of the proton using lattice QCD simulations at physical pion mass
C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, H. Panagopoulos, G. Spanoudes
CComplete flavor decomposition of the spin and momentum fraction of the protonusing lattice QCD simulations at physical pion mass
C. Alexandrou , , S. Bacchio , M. Constantinou , J. Finkenrath ,K. Hadjiyiannakou , , K. Jansen , G. Koutsou , H. Panagopoulos , G. Spanoudes (Extended Twisted Mass Collaboration) Department of Physics, University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus Computation-based Science and Technology Research Center,The Cyprus Institute, 20 Kavafi Str., Nicosia 2121, Cyprus Department of Physics, Temple University, 1925 N. 12th Street, Philadelphia, PA 19122-1801, USA NIC, DESY, Platanenallee 6, D-15738 Zeuthen, Germany
We evaluate the gluon and quark contributions to the spin of the proton using an ensemble of gaugeconfiguration generated at physical pion mass. We compute all valence and sea quark contributionsto high accuracy. We perform a non-perturbative renormalization for both quark and gluon matrixelements. We find that the contribution of the up, down, strange and charm quarks to the protonintrinsic spin is (cid:80) q = u,d,s,c ∆Σ q + = 0 . (cid:80) q = u,d,s,c J q + = 0 . J g = 0 . J = J q + J g = 0 . . . . PACS numbers:Keywords:
I. INTRODUCTION
The spin decomposition of the proton reveals impor-tant information about its non-perturbative structure.Since the proton is composed of quarks and gluons, itis expected that its spin arises from the intrinsic spinand orbital angular momentum of its constituents. Thefirst attempts to measure the proton spin were performedat SLAC in the E80 [1, 2] and E130 [3, 4] series of exper-iments. The successful quark model that describes wellproperties of the low-lying hadrons predicted that all thespin is carried by the three valence quarks. The first ma-jor surprise came from the measurements of the EuropeanMuon Collaboration (EMC) [5, 6] that determine the pro-ton spin-dependent structure function down to x = 0 . s + . Phenomenological analyses point to a negativevalue but the error is large, giving values of ∆Σ s + rang-ing from − . − . q + = q + ¯ q to denotethe sum from quark and antiquark contributions to theintrinsic spin and momentum fraction. Results from in-clusive DIS experiments have, however, small sensitivityto the gluon helicity ∆ g . In contrast, polarized proton-proton collisions, in particular jet or hadron productionat high transverse momentum available from the Rela-tivistic Heavy Ion Collider (RHIC) [19–21] at BNL pro-vide tighter constraints on (cid:82) . . ∆ g ( x ) dx = 0 . +0 . − . .Despite the tremendous progress in the determination ofthe gluon helicity, large uncertainties remain mostly inthe small- x range. Thanks to its large kinematic reachin x and Q , the planned Electron-Ion Collider(EIC) [22]will provide significantly more input to constrain ∆ g .While experiments play a crucial role in the under-standing of the sources of the proton spin, they need tobe complemented by phenomenological analyses, whichinvolve model dependence and parameterizations. Lat-tice QCD (LQCD), on the other hand, provides the initio non-perturbative framework that is suitable to addressthe key questions of how the nucleon spin and momen- a r X i v : . [ h e p - l a t ] M a r tum is distributed among its constituents using directlythe QCD Lagrangian. Tremendous progress has beenmade in simulating lattice QCD in recent years. State-of-the-art simulations are being performed with dynamicalup, down and strange quarks with mass tuned to theirphysical values (referred to as the physical point). A sub-set of simulations also include a dynamical charm quarkwith mass fixed to approximately its physical value. Thisprogress was made possible using efficient algorithms andin particular multigrid solvers [23] that were developedfor twisted mass fermions [24].A number of recent lattice QCD studies were car-ried out to extract the intrinsic spin carried by eachquark flavor. They include previous works by the Ex-tended Twisted Mass Collaboration (ETMC) [25–27], byPNDME [28] and by χ QCD [29]. First attempts to com-pute the gluon average momentum fraction were car-ried out by the pioneering work of the QCDSF collab-oration [30, 31] in the quenched approximations. Re-sults on the gluon momentum fraction using dynamicalgauge field configurations appeared only recently. Theyused mostly simulations with larger than physical pionmass relying on chiral extrapolations to obtain final re-sults [32–34]. A first attempt to fully decompose the nu-cleon spin was carried out by χ QCD [35] in the quenchedapproximation, followed by a study of the gluon spin [36]using 2 + 1 dynamical fermions on four lattice spacingsand four volumes including an ensemble with physicalvalues for the quark masses. ETMC was the first to com-pute the gluon momentum fraction directly at the physi-cal point without the need of a chiral extrapolation [26].The latter is a significant progress, since chiral extrapo-lations in the nucleon sector introduce uncontrolled sys-tematic errors.In this study we will provide the complete quark fla-vor decomposition of the proton spin. This requires thecomputation of both valence and sea quark contributions.It also includes the computation of the gluon contribu-tions to the spin and momentum fraction of the pro-ton. In order to evaluate the quark loop contributionsthat are computationally very demanding, we apply im-proved techniques that are developed and implementedon graphics cards (GPUs) [37], as well as noise reduc-tion techniques [38, 39]. This work updates our pre-vious results on the proton spin presented in Ref. [26]in several respects: i) While Ref. [26] used an ensembleof twisted mass fermions generated with two degeneratelight quarks ( N f = 2) [40], we here use an ensemble oftwisted mass fermions [41, 42] that includes, besides thelight quarks, the strange and the charm quarks all withmasses fixed to their physical values ( N f = 2 + 1 + 1); ii)we perform a more elaborated analysis of excited statecontributions; iii) we use larger statistics; iv) we computethe gluon contribution to the proton spin taking into ac-count the generalized form factor B (0); and v) we usenon-perturbative renormalization not only for the quarkoperators but also for the gluon operator.The remainder of this paper is organized as follows: In Section II we provide the theoretical basis for the nu-cleon spin decomposition [43]. Sections III, IV describethe methodology to extract the nucleon bare matrix el-ements needed, while Section V explains the renormal-ization procedure and the conversion to the MS scheme.Our final results are discussed in Section VI and com-pared with other studies in Section VII. Finally, in Sec-tion VIII we summarize our findings and conclude. II. NUCLEON SPIN DECOMPOSITION
A key object for the study of the spin decompositionis the QCD energy-momentum tensor (EMT) T µν . Thesymmetric part of the EMT can be separated [44] intotwo terms, the traceless term, denoted by ¯ T µν , and thetrace term ˆ T µν as T µν = ¯ T µν + ˆ T µν . (1)Only the traceless part is relevant for this study. Keep-ing only the gauge-invariant parts of ¯ T µν , this can beexpressed in terms of the gluon part ¯ T µν ; g and the quarkpart ¯ T µν ; q as ¯ T µν = ¯ T µν ; g + ¯ T µν ; q . (2)where ¯ T µν ; g = F { µρ F ν } ρ (3)and ¯ T µν ; q = ¯ ψiγ { µ ←→ D ν } ψ , (4)where F µν is the gluon field-strength tensor and thenotation {· · · } means symmetrization over the indicesin the parenthesis and subtraction of the trace. Thesymmetrized covariant derivative ←→ D is defined as ←→ D =( ←− D + −→ D ) / M ij can be writtenin terms of the EMT as M αµν = ¯ T αν x µ − ¯ T αµ x ν (5)and the i-th component of the angular momentum oper-ator as J i = 12 (cid:15) ijk (cid:90) d xM jk ( x ) . (6)Substituting Eq.(3) into Eq.(6), as discussed in Refs. [43,45], the gauge invariant gluon angular momentum oper-ators is (cid:126)J g = (cid:90) d x ( (cid:126)x × ( (cid:126)E × (cid:126)B )) (7)where (cid:126)E and (cid:126)B are the chromo-electric and chromo-magnetic fields. Substituting Eq.(4) into Eq.(6), we ob-tain the gauge-invariant quark angular momentum oper-ator [43, 45], (cid:126)J q = (cid:90) d x (cid:20) ¯ ψ (cid:126)γγ ψ + ¯ ψ ( (cid:126)x × i −→ D ) ψ (cid:21) . (8)The first term in Eq.(8) is the quark intrinsic spin op-erator and the second term is the orbital angular mo-mentum. Putting gluon and quark operators togetherwe have that (cid:126)J = (cid:126)J g + (cid:126)J q = (cid:126)J g + (cid:32) (cid:126) Σ q (cid:126)L q (cid:33) . (9)This is the so called Ji decomposition [43] which doesnot allow to decompose J g any further in a gauge invari-ant manner. Jaffe and Manohar suggested a non gauge-invariant way to decompose further [46] the gluon angularmomentum, with the issue of the gauge-invariance beingaddressed in Ref. [47]. In this work we use Ji’s decom-position [43] and thus compute the total gluon angularmomentum J g .In order to compute the nucleon spin, we need to evalu-ate the nucleon matrix elements of the EMT. They can bedecomposed into three generalized form factors (GFFs)in Minkowski space as follows [45] (cid:104) N ( p (cid:48) , s (cid:48) ) | T µν ; q,g | N ( p, s ) (cid:105) = ¯ u N ( p (cid:48) , s (cid:48) ) (cid:20) A q,g ( q ) γ { µ P ν } + B q,g ( q ) iσ { µρ q ρ P ν } m N + C q,g ( q ) q { µ q ν } m N (cid:21) u N ( p, s ) (10)where u N the nucleon spinor with initial (final) momen-tum p ( p (cid:48) ) and spin s ( s (cid:48) ), P = ( p (cid:48) + p ) / q = p (cid:48) − p the momentum transfer. A ( q ), B ( q ) and C ( q ) are the three GFFs. In the forwardlimit, A q,g (0) gives the quark and gluon average momen-tum fraction (cid:104) x (cid:105) q,g . Summing over all quark and gluoncontributions gives the momentum sum (cid:104) x (cid:105) q + (cid:104) x (cid:105) g = 1.As shown in Ref. [43] the nucleon spin can be written interms of A and B in the forward limit as J = 12 [ A q (0) + A g (0) + B q (0) + B g (0)] , (11) where we consider a reference spin axis. The spin sum J = together with the momentum sum are satisfiedif B q (0) + B g (0) = 0. Although A q,g (0) and thus theaverage momentum fractions are extracted from the nu-cleon matrix element directly at zero momentum trans-fer, B q,g (0) can only be computed at non-zero momentumtransfer requiring an extrapolation to Q = 0.Since we have a direct way to compute the quark con-tribution J q and the intrinsic spin ∆Σ q we can deter-mine the quark orbital angular momentum by L q = J q −
12 ∆Σ q . (12) III. COMPUTATION OF THE BARE NUCLEONMATRIX ELEMENTSA. Ensembles of gauge configurations
In Table I we give the parameters of the N f =2 + 1 + 1 ensemble analyzed in this work denoted bycB211.072.64 [48]. For completeness we also list the pa-rameters of the N f = 2 ensemble analyzed in our previ-ous study [26], referred to as cA2.09.48. In both casesthe lattice spacing is determined using the mass of thenucleon [48–51].The ensembles are produced using the Iwasaki [52] im-proved gauge action and the twisted mass fermion formu-lation [41, 42]. A clover term [53] was added to stabilizethe simulations. The twisted mass fermion formulation isvery well suited for hadron structure providing an auto-matic O ( a ) improvement [42] with no need of improvingthe operators. TABLE I: Simulation parameters for the cB211.072.64 [48] and cA2.09.48 [40] ensembles. c SW is the value of the clovercoefficient, β = 6 /g where g is the coupling constant, N f is the number of dynamical quark flavors in the simulation, a isthe lattice spacing, V the lattice volume in lattice units, m π the pion mass, m N the nucleon mass, and L the spatial latticelength in physical units. For the parameters of the cA2.09.48 ensemble, the second error arises from the systematic error on thedetermination of the lattice spacing due to the extrapolation to the physical value of m π [40]. For the cB211.072.64 ensemblethis systematic error is negligible.ensemble c SW β N f a [fm] V am π m π L am N m N /m π m π [GeV] L [fm]cB211.072.64 1.69 1.778 2+1+1 0.0801(4) 64 ×
128 0.05658(6) 3.62 0.3813(19) 6.74(3) 0.1393(7) 5.12(3)cA2.09.48 1.57551 2.1 2 0.0938(3)(1) 48 ×
96 0.06208(2) 2.98 0.4436(11) 7.15(2) 0.1306(4)(2) 4.50(1)
B. Construction of correlation functions
To compute the nucleon matrix elements one needsto evaluate two- and three-point functions in Euclidean space. To create states with the quantum numbers of thenucleon we use as interpolating field J N ( t, (cid:126)x ) = (cid:15) abc u a ( x ) (cid:2) u bT ( x ) C γ d c ( x ) (cid:3) (13)where u ( x ) , d ( x ) are the up, down quark fields and C isthe charge conjugation matrix. The interpolating field inEq. (13) does not only create the nucleon state but alsoexcited states with the quantum numbers of the nucleon,including multi-particle states. In order to increase theoverlap of the interpolating field with the ground statewe employ Gaussian smearing [54, 55] on the quark fieldsas well as APE smearing [56] on the gauge links enteringthe hopping matrix of the smearing function. For moredetails about how we tune these smearing parameters aregiven in Refs. [49, 51]. The nucleon two-point function isgiven by C (Γ , (cid:126)p ; t s , t ) = (cid:88) (cid:126)x s e − i ( (cid:126)x s − (cid:126)x ) · (cid:126)p × Tr (cid:2) Γ (cid:104)J N ( t s , (cid:126)x s ) ¯ J N ( t , (cid:126)x ) (cid:105) (cid:3) , (14)where x is the initial lattice site at which states with thequantum numbers of the nucleon are created, referredto as source position and x s is the site where they areannihilated, referred to as sink. An appropriate operator O µν probes the quarks and gluons within the nucleonat a lattice site x ins referred to as insertion point. Theresulting three-point function is given by C µν (Γ ρ , (cid:126)q, (cid:126)p (cid:48) ; t s , t ins , t )= (cid:88) (cid:126)x ins ,(cid:126)x s e i ( (cid:126)x ins − (cid:126)x ) · (cid:126)q e − i ( (cid:126)x s − (cid:126)x ) · (cid:126)p (cid:48) × Tr (cid:2) Γ ρ (cid:104)J N ( t s , (cid:126)x s ) O µν ( t ins , (cid:126)x ins ) ¯ J N ( t , (cid:126)x ) (cid:105) (cid:3) . (15)The operator O µν may represent the EMT with twoLorentz indices or the helicity operator with one Lorentzindex. The Euclidean momentum transfer squared isgiven Q = − ( p (cid:48) − p ) and Γ ρ is the projector actingon the spin indices. We consider Γ = (1 + γ ) andΓ k = i Γ γ γ k taking the non-relativistic representationof γ µ . C. Analysis of correlation functions to extract thenucleon matrix elements
The information about the desired nucleon matrix ele-ment is contained in the three-point correlation functionof Eq.(15). In order to extract it, we construct appropri-ate combinations of three- to two-point functions, whichin the large Euclidean time limit, cancel the time depen-dence arising from the time propagation and the overlapterms between the interpolating field and the nucleonstate. An optimal choice that benefits from correlationsis the ratio [57–60]. R µν (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ; t s , t ins ) = C µν (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ; t s , t ins ) C (Γ , (cid:126)p (cid:48) ; t s ) × (cid:115) C (Γ , (cid:126)p ; t s − t ins ) C (Γ , (cid:126)p (cid:48) ; t ins ) C (Γ , (cid:126)p (cid:48) ; t s ) C (Γ , (cid:126)p (cid:48) ; t s − t ins ) C (Γ , (cid:126)p ; t ins ) C (Γ , (cid:126)p ; t s ) . (16) The sink and insertion time separations t s and t ins aretaken relative to the source. In the ratio of Eq.(16), tak-ing the limits ( t s − t ins ) (cid:29) a and t ins (cid:29) a , with a thelattice spacing, the nucleon state dominates. When thishappens, the ratio becomes independent of time R µν (Γ ρ ; (cid:126)p (cid:48) , (cid:126)p ; t s ; t ins ) t s − t ins (cid:29) a −−−−−−−→ t ins (cid:29) a Π µν (Γ ρ ; (cid:126)p (cid:48) , (cid:126)p ) (17)and yields the desired nucleon matrix element. In prac-tice, ( t s − t ins ) and t ins cannot be taken arbitrarily large,since the signal-to-noise ratio decays exponentially withthe sink-source time separation. Therefore, one needsto take ( t s − t ins ) and t ins large enough so that the nu-cleon state dominates in the ratio. To identify when thishappens is a delicate process. We use three methods tocheck for convergence to the nucleon state as summarizedbelow. Plateau method:
The ratio of Eq. (16) can be written asΠ µν (Γ ρ ; (cid:126)p (cid:48) , (cid:126)p ) + O ( e − ∆ E ( t s − t ins ) ) + O ( e − ∆ Et ins ) + · · · (18)where the first term is time-independent and contribu-tions from excited states are exponentially suppressed.In Eq. (18) ∆ E is the energy gap between the nucleonstate and the first excited state. In order to extract thenucleon matrix element of the operator of interest, weseek to identify nucleon state dominance by looking fora range of values of t ins for which the ratio of Eq. (16)is time-independent (plateau region). We fit the ratioto a constant within the plateau region and seek to seeconvergence in the extracted fit values as we increase t s .If such a convergence can be demonstrated, then the de-sired nucleon matrix element can be extracted. Summation method:
One can sum over t ins the ratio ofEq. (16) [61, 62] to obtain R µν summed (Γ ρ ; (cid:126)p (cid:48) , (cid:126)p ; t s ) = t s − a (cid:88) t ins =2 a R µν (Γ ρ ; (cid:126)p (cid:48) , (cid:126)p ; t s ; t ins ) = c + Π µν (Γ ρ ; (cid:126)p (cid:48) , (cid:126)p ) × t s + O ( e − ∆ Et s ) . (19)Assuming the nucleon state dominates, Π µν (Γ ρ ; (cid:126)p (cid:48) , (cid:126)p ) isextracted from the slope of a linear fit with respect to t s .As in the case of the plateau method, we probe conver-gence by increasing the lower value of t s , denoted by t low s used in the linear fit, until the resulting value converges.While both plateau and summation methods assume thatthe ground state dominates, the exponential suppressionof excited states in the summation is faster and approx-imately corresponds to using twice the sink-source timeseparation t s in the plateau method. Two-state fit method:
In this method we explicitly in-clude the contributions from the first excited state. Wethus expand the two- and three-point function correlatorsentering in the ratio of Eq.(16) to obtain C ( (cid:126)p, t s ) = c ( (cid:126)p ) e − E ( (cid:126)p ) t s + c ( (cid:126)p ) e − E ( (cid:126)p ) t s + · · · , (20)and C µν (Γ ρ , (cid:126)p (cid:48) , (cid:126)p, t s , t ins ) = A µν , (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ) e − E ( (cid:126)p (cid:48) )( t s − t ins ) − E ( (cid:126)p ) t ins + A µν , (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ) e − E ( (cid:126)p (cid:48) )( t s − t ins ) − E ( (cid:126)p ) t ins + A µν , (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ) e − E ( (cid:126)p (cid:48) )( t s − t ins ) − E ( (cid:126)p ) t ins + A µν , (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ) e − E ( (cid:126)p (cid:48) )( t s − t ins ) − E ( (cid:126)p ) t ins + · · · , (21)where c ( (cid:126)p ) and c ( (cid:126)p ) are the overlaps of the groundand first excited state with the interpolating field and E ( (cid:126)p ) and E ( (cid:126)p ) the corresponding energies. The pa-rameters A µνi,j , are the matrix elements of the i, j statesmultiplied by the corresponding overlap terms. Note that A µν , (cid:54) = A µν , for non-zero momentum transfer. Our pro-cedure to determine these parameters is as follows: wefirst fit the effective mass using the two-point functionof Eq. (20) at (cid:126)p = (cid:126) (cid:126)p to ex-tract the nucleon mass m N , c ( (cid:126)p ) /c ( (cid:126)p ) and ∆ E ( (cid:126)p ) = E ( (cid:126)p ) − E ( (cid:126)p ), where for E ( (cid:126)p ) we use the dispersion re-lation E ( (cid:126)p ) = (cid:112) m N + (cid:126)p . In Fig. 1 we compare theenergy E ( p ) extracted directly from the finite momen-tum two-point function and the dispersion relation. Ascan be seen, the dispersion relation is well satisfied andholds for all the momenta that are used in this study. In- [GeV ]0.91.01.11.21.31.4 E N ( p ) [ G e V ] FIG. 1: Red points show the energy of the nucleon E N ( (cid:126)p ) inGeV as extracted from finite momentum two-point functionsand the grey band shows the dispersion relation E N ( (cid:126)p ) = (cid:112) m N + (cid:126)p as a function of (cid:126)p in GeV . serting the expressions of Eqs. (20) and (21) in Eq. (16)and using E , ∆ E and c /c extracted from the two-point correlators we fit the resulting ratio to extract theremaining four parameters, M µν , ≡ A µν , / (cid:112) c ( (cid:126)p (cid:48) ) c ( (cid:126)p ), A µν , / A µν , , A µν , / A µν , and A µν , / A µν , . The first parame-ter M µν , is the desired ground state matrix element andthe rest are excited states contributions. In the case of zero momentum transfer, A µν , = A µν , and we only havethree parameters to determine.For all the three methods we minimize χ defined as χ = r T C r , (22)where C is the covariance matrix, r = y − f ( p , x ) theresidual vector between our data y and the model func-tion f ( p , x ) and p a vector holding the parameters of thefit.After determining the nucleon matrix element we canextract the charges, moments and GFFs. For zero mo-mentum transfer we haveΠ i (Γ k ; (cid:126) ,(cid:126)
0) = g A δ ik (23)in the case of the helicity operator. The average momen-tum fraction can be obtained from the matrix elementof the one-derivative vector operator at zero momentumtransfer, Π (Γ ; (cid:126) ,(cid:126)
0) = − m N (cid:104) x (cid:105) , (24)Π i (Γ ; p i , p i ) = ip i (cid:104) x (cid:105) , (25)where in the second expression the nucleon is boostedin the i th direction with momentum p i . As already dis-cussed in connection to Eq. (10), B (0) cannot be ex-tracted directly. We thus compute B ( Q ) as a functionof Q and extrapolate to zero Q . More details about theprocedure to extract B ( Q ) will be given in Sec. IV B. D. Connected and disconnected three-pointfunctions and statistics
The three-point function, defined in Eq. (15), receivescontribution from two types of diagrams: one when theoperator couples directly to a valence quark, known asthe connected contribution and one when the operatorcouples to a sea quark resulting into a quark loop, knownas the disconnected contribution. The gluon operator asdefined in Eq. (3), produces a closed gluon loop and thusa disconnected contribution. To evaluate the connectedcontributions we use standard techniques that involvethe computation of the sequential propagator throughthe sink. In this approach the sink-source time separa-tion, the projector and the momentum at the sink (cid:126)p (cid:48) arekept fixed. We perform the computation of the connectedthree-point functions fixing the sink momentum to zero,i.e we set (cid:126)p (cid:48) = (cid:126)
0. We then compute the sequential prop-agator for both the unpolarized and polarized projectorsΓ and Γ k , k = 1 , ,
3, respectively.In total we analyze 750 configurations separated by 4trajectories. We use seven values of the sink-source timeseparation t s ranging from 0.64 fm to 1.60 fm. In order tokeep the signal-to-noise ratio approximately constant weincrease the number of source positions as we increase t s .The statistics used for each value of t s are given in Ta-ble II. A range of t s and the increasingly larger statisticsfor larger t s allow us to better check excited state effectsand to thus reliably extract the nucleon matrix elementsof interest. TABLE II: Parameters used for the evaluation of the con-nected three-point functions. In the first column we give thevalue of t s in lattice units and in the second column in physicalunits. In all cases 750 gauge configurations are analyzed. Inthe third column we give the number of source positions, andin the fourth column the total number of measurements. Thelast column gives N inv , which is the total number of inversionsper configuration. t s /a t s [fm] N srcs N meas N inv The disconnected contribution involves the discon-nected quark loop correlated with the nucleon two-pointcorrelator. The disconnected quark loop is given by L ( t ins , (cid:126)q ) = (cid:88) (cid:126)x ins Tr (cid:2) D − ( x ins ; x ins ) G (cid:3) e i(cid:126)q · (cid:126)x ins , (26)where D − ( x ins ; x ins ) is the quark propagator that startsand ends at the same point x ins and G is an appropriatelychosen γ -structure. For the helicity operator we use (cid:126)γγ and for the quark part of EMT γ µ ←→ D ν . A direct compu-tation of the quark loops would need inversions from allspatial points on the lattice, making the evaluation unfea-sible for our lattice size. We, therefore, employ stochas-tic techniques combined with dilution schemes [63] thattake into account the sparsity of the Dirac operator andits decay properties. Namely, we employ the hierarchi-cal probing technique [38], which provides a partitioningscheme that eliminates contributions from neighboringpoints in the trace of Eq. (26) up to a certain coloringdistance 2 k . Using Hadamard vectors as the basis vec-tors for the partitioning, one needs 2 d ∗ ( k − vectors,where d =4 for a 4-dimensional partitioning. Note thatthe computational resources required are proportional tothe number of Hadamard vectors, and therefore, in d =4dimensions, increase 16-fold each time the probing dis-tance 2 k doubles. Contributions entering from points be-yond the probing distance are expected to be suppresseddue to the exponential decay of the quark propagator andare treated with standard noise vectors that suppress alloff-diagonal contributions by 1 / √ N r , i.e.1 N r (cid:88) r | ξ r (cid:105)(cid:104) ξ r | = 1 + O (cid:18) √ N r (cid:19) , (27)where N r is the size of the stochastic ensemble. Hier-archical probing has been employed with great success in previous studies for an ensemble with a pion massof 317 MeV [64]. For simulations at the physical point,it is expected and confirmed [49] that a larger probingdistance is required since the light quark propagator de-cays more slowly because of the smaller quark mass. Weavoid the need of increasing the distance by combininghierarchical probing with deflation of the low modes [65].Namely, for the light quarks we construct the low modecontribution to the quark loops by computing exactlythe smallest eigenvalues and corresponding eigenvectorsof the squared Dirac operator and combine them withthe contribution from the remaining modes, which areestimated using hierarchical probing. Additionally, weemploy the one-end trick [66], also used in our previousstudies [26, 27, 37] and fully dilute in spin and color.For the calculation of the nucleon matrix element of thegluonic part of the EMT, we use the gluon field strengthtensor F µν ( x ) = i g (cid:20) U µ ( x ) U ν ( x + a ˆ µ ) U † µ ( x + a ˆ ν ) U † ν ( x )+ U ν ( x ) U † µ ( x + a ˆ ν − a ˆ µ ) U † ν ( x − a ˆ µ ) U µ ( x − a ˆ µ )+ U † µ ( x − ˆ µ ) U † ν ( x − a ˆ ν − a ˆ µ ) U µ ( x − a ˆ ν − a ˆ µ ) U ν ( x − a ˆ ν )+ U † ν ( x − a ˆ ν ) U µ ( x − a ˆ ν ) U ν ( x − a ˆ ν + a ˆ µ ) U † µ ( x ) − h.c (cid:21) , (28)with g the bare coupling constant. For the gauge linksentering the field strength tensor we apply stout smear-ing [67] with parameter ρ = 0 .
129 [33]. As will be dis-cussed in Sec. IV A, we investigate the signal-to-noise ra-tio as we increase the number of stout steps.While the evaluation of the gluon loop is computa-tionally cheap because no inversions are needed, the cal-culation of the quark loops is very expensive. We usethe same combination of methods for the calculation ofthe flavors of the quark loops except for deflation, whichis only used for the light quark loops. The parametersused for the evaluation of the quark loops are collectedin Table III. Two hundred low modes of the square Diracoperator are computed in order to reduce the stochasticnoise in the computation of the light quark loops. For thecharm quark we use a coloring distance 2 in hierarchicalprobing instead of 2 used for the light and the strangequark loop since charm quarks are relatively heavy. In-stead we compute 12 stochastic vectors. We evaluatethe nucleon two-point functions for two hundred ran-domly chosen source positions which sufficiently reducesthe gauge noise for large enough sink-source time sepa-rations of the disconnected three-point functions. Sincethey are available we use the same number of two pointfunctions for all source-sink time separations.In summary, we perform in total 12,690,000 inversionsfor the connected and 16,272,000 for the disconnectedcontributions by employing the DD- α AMG solver and itsQUDA version [24, 68, 69] to accelerate the inversions.Using GPUs and the DD- α AMG solver are essential toobtain the required statistics.
TABLE III: Parameters and statistics used for the evaluationof the disconnected three-point functions. The number of con-figurations analyzed is N cnfs = 750 and the number of sourcepositions used for the evaluation of the two-point functions is N srcs = 200 per gauge configuration. In the case of the lightquarks, we compute the lowest 200 modes exactly and deflatebefore computing the higher modes stochastically. N r thenumber of noise vectors, and N Had the number of Hadamardvectors. N sc = 12 corresponds to spin-color dilution and N inv is the total number of inversions per configuration.Flavor N def N r N Had N sc N inv light 200 1 512 12 6144strange 0 1 512 12 6144charm 0 12 32 12 4608 IV. BARE NUCLEON MATRIX ELEMENTS
As already discussed, for the decomposition of the nu-cleon spin we need the axial charges for each quark flavor,which give the quark helicities, the average momentumfractions and B (0). The extraction of the axial chargesor ∆Σ q for the cB211.072.64 ensemble is presented inRef. [25], while the evaluation of the isovector A and B in Ref. [51]. In this section we focus in the extractionof the remaining quantities needed for the full decompo-sition of the nucleon spin. A. Average momentum fraction (cid:104) x (cid:105) The average momentum fraction (cid:104) x (cid:105) is extracted di-rectly from the nucleon matrix element at zero momen-tum transfer from the Eqs. (24), (25). In the case ofthe connected contribution, since we only have access tothree-point functions with (cid:126)p (cid:48) = (cid:126)
0, we are restricted touse Eq. (24). In Fig. 2 we show the bare ratio whichleads to the extraction of the connected contribution tothe isoscalar (cid:104) x (cid:105) u + + d + B . The ratio for each t s has beenconstructed between three- and two-point functions withthe same source positions to benefit from the correlationsbetween numerator and denominator. One can easily ob-serve a clear contamination from excited states at smalltime separations. In fact for t s < t s increases.We note that we exclude t ins /a = 0 , , t s /a − , t s /a sincethey do not carry physical information. For t s (cid:38) .
12 fmwe fit within the plateau region discarding five pointsfrom left and right, thus t ins ∈ [2 + τ, t s − τ −
2] with τ /a = 5 for all t s values. This range is found to yield agood χ / d . o . f . .In Fig. 3 we show the summed ratio of Eq. (19). Ascan be seen, a linear fit describes well the results. ins t s /2 [fm]0.450.500.550.600.650.700.750.80 R ( t s , t i n s ) → › x fi u + + d + B ( C o nn . ) t s =0.64 fmt s =0.80 fmt s =0.96 fm t s =1.12 fmt s =1.28 fm t s =1.44 fmt s =1.60 fm FIG. 2: The ratio of Eq. (16) for zero momentum from wherethe connected contribution to (cid:104) x (cid:105) u + + d + B using Eq. (24) is ex-tracted as a function of t ins for source-sink time separations t s /a = 8 , , , , , ,
20 using blue circles, orange downtriangles, up green triangles, left red triangles, right purpletriangles, brown rhombus, magenta crosses, respectively. Thebands show a constant fit to the points within the range ofthe band. s [fm]45678 X t i n s R ( t s , t i n s ) → › x fi u + + d + B ( C o nn . ) FIG. 3: The summed ratio of Eq. (19) is shown as a functionof t s (red circles). The slope yields the connected (cid:104) x (cid:105) u + + d + B .The linear fits are shown by the blue bands as we increase,from top to bottom, the smallest value of t s used in the fit, t low s . Our approach for the two-state fits has been discussedin Sec. III C. We extract m N , ∆ E and the overlap ra-tio c /c using the full statistics of two-point functionproduced with 264 source positions. However, as dis-cussed above, for the construction of the ratio we use thesame source positions for three- and two-point functions.We fit the resulting ratios simultaneously for all values t s ≥ t low s . We vary t low s , to check convergence of the ex-tracted nucleon matrix element. We show the resultingfits to the ratios in Fig. 4 for t low s /a = 12. Additionally,we plot in the middle panel the predicted time depen-dence of the ratio when fixing t ins = t s / t low s /a = 12.In the same panel we include values extracted using theplateau method. For t s /a <
14 where no plateau could beidentified we plot the midpoint t s /
2. As can be seen, thetwo-state fit predicts well the residual time-dependence ofvalues extracted using the plateau method. It also showsthat the plateau values, even for t s = 1 . t s > t s ∼ a = 2 .
08 fm requiringan order of magnitude more statistics as compared to the statistics used for t s = 20 a = 1 . t low s . They also agree with the value extractedfrom the summation method when t low s = 16 a = 1 .
28 fmis used in the fit, as shown in the right panel of Fig. 4.We thus take as our final determination of the connectedbare (cid:104) x (cid:105) u + + d + B the value extracted from the two-state fitfor t low s /a = 12 with χ / d . o . f = 1 .
2. The selected valueis shown by the horizontal grey band spanning the wholerange of Fig. 4. We find for the connected bare isoscalarmomentum fraction (cid:104) x (cid:105) u + + d + B = 0 . . (29) t ins t s /2 [fm]0.250.350.450.550.650.75 › x fi u + + d + B ( C o nn . ) t s [fm] 0.6 0.8 1.0 1.2 1.4 t lows [fm] . . . . . . . . . . FIG. 4: Excited state analysis for determining the connected isoscalar average momentum fraction (cid:104) x (cid:105) u + + d + B using Eq. (24).In the left panel, we show results for the ratio of Eq. (16) with symbol and color notation as in Fig. 2. The results are shownas a function of the insertion time t ins shifted by t s /
2. The dotted lines and associated error bands are the resulting two-statefits. In the middle panel, we show the plateau values or middle point when no plateau is identified, as a function of source-sinkseparation using the same symbol used for the ratio in the left panel for the same t s . The grey band is the predicted time-dependence of the ratio using the parameters extracted from the two-state fit when t low s = 12 a = 0 .
96 fm. In the right panel,we show values of the connected (cid:104) x (cid:105) u + + d + B extracted using the two-state fit (black squares) and the summation method (greenfilled triangles) as a function of t low s together with the χ / d . o . f for each fit. The open symbol shows the selected value for theconnected (cid:104) x (cid:105) u + + d + B with the grey band spanning the whole range of the figure being its statistical error. t ins t s /2 [fm]0.100.150.200.250.30 › x fi u + d + B ( C o nn . ) t s [fm] 0.6 0.8 1.0 1.2 1.4 t lows [fm] . . . . . . . . . . FIG. 5: Excited state analysis for determining the isovector average momentum fraction (cid:104) x (cid:105) u + − d + B using Eq. (24). The notationfollows that in Fig. 4. The isovector average momentum fraction (cid:104) x (cid:105) u + − d + B for the cB211.072.64 ensemble is reported in Ref. [51]. t ins t s /2 [fm]0.000.050.100.150.20 › x fi u + + d + B ( D i s c . ) t s [fm] 0.5 0.6 0.7 t lows [fm] FIG. 6: Excited state analysis for determining the disconnected contribution to the isoscalar average momentum fraction usingEq. (25). The notation is the same as that in Fig. 4. The sink-source time separations shown are t s /a = 6 , , , , , , , t ins t s /2 [fm]0.000.020.040.060.08 › x fi s + B t s [fm] 0.5 0.6 0.7 t lows [fm] FIG. 7: Excited state analysis for determining the strange average momentum fraction as extracted from Eq.(25). The notationis the same as that in Fig. 6.
For completeness and easy reference we repeat the anal-ysis following the same procedure as for the connected (cid:104) x (cid:105) u + + d + B . The results are shown in Fig. 5. As can beenseen, the effect of excited states is similar to what isobserved for the connected isoscalar case. The value de-termined from the two-state fit for t low s /a = 8 varies onlyvery mildly as we increase t low s and it is in agreementwith the value extracted from the summation methodfor t low s /a = 14. We thus select as our final value theone extracted from the two-state fit using t low s /a = 8,obtaining (cid:104) x (cid:105) u + − d + B = 0 . . (30)In Figs. 6 and 7 we present our results for the quarkdisconnected contributions to the isoscalar and strangeaverage momentum fractions. Since for the disconnectedcontribution one can use a boosted nucleon without theneed of additional inversions, one can extract it both fromthe diagonal part of the EMT as in Eq. (24) and from thenon-diagonal as in Eq. (25). If we use the diagonal part ofEMT, there is a large non-zero vacuum expectation value,which, although it cancels after the trace subtraction,it leads to large statistical fluctuations. In the case of Eq. (25), where the off-diagonal components enter, thisproblem does not arise. We thus boost the nucleon usingthe first non-zero momentum, namely (cid:126)p = ˆ n π/L withˆ n = (1 , ,
0) and all other permutations and we averageover the three directions and two orientations to obtaina good signal-to-noise ratio as presented in Figs. 6 and 7.Unlike the connected contributions, for both light dis-connected and strange, the ratios show fast convergence,indicating that excited states are suppressed, within ourstatistical uncertainties. We thus perform a fit within theplateau range that includes t ins ∈ [3 a, t s − a ]. The valuesextracted from the plateau fits converge to a constant andare consistent with the results extracted from the sum-mation method. We take the weighted average of theconverged plateau values, namely for the plateau valuesextracted for t s > . t low s in the range[6 a, a ] corroborating the fact that for these quantitiesexcited states contamination is suppressed compared tothe statistical error. We find for the disconnected contri-0bution to the isoscalar average momentum fraction (cid:104) x (cid:105) u + + d + B = 0 . (cid:104) x (cid:105) s + B = 0 . (cid:104) x (cid:105) c + B = 0 . , (33)which is compatible with zero. t ins t s /2 [fm]0.40.20.00.20.40.60.81.01.2 › x fi g B ( ¯ T ; g ) t s [fm] 0.5 0.6 0.7 0.8 0.9 t lows [fm]0.6 0.4 0.2 0.0 0.2 0.4 0.6 t ins t s /2 [fm]0.40.20.00.20.40.60.81.01.2 › x fi g B ( ¯ T i ; g ) t s [fm] 0.5 0.6 0.7 0.8 0.9 t lows [fm] FIG. 8: Excited state analysis for determining the gluon average momentum fraction (cid:104) x (cid:105) gB . In the upper panel we show (cid:104) x (cid:105) gB extracted using Eq.(24) and in the lower panel using Eq.(25). In both cases we use stout smearing with n S = 10 steps. Thenotation is the same as that in Fig. 6. We show results for t s /a = 6 , , , , , ,
18 with blue circles, orange down triangles,up green triangles, left red triangles, right purple triangles, brown rhombus, and magenta crosses, respectively.
In Fig. 8 we present our analysis for the gluon aver-age momentum fraction. For this case we employ stoutsmearing on the gauge links entering in the field strengthtensor of Eq. (28) to improve the signal of the gluonicpart of the EMT. We show the case where the numberof stout steps n S is 10. We analyze both the diagonaland off-diagonal components of EMT given by Eqs. (24)and (25). When using the diagonal components, due tothe subtraction of a large trace, large gauge fluctuationsare observed. This is analogous to the quark discon-nected contributions discussed above. Although the vac-uum expectation value for the traceless part of the EMT, (cid:104) | ¯ T g | (cid:105) , is compatible with zero, as expected by thesubtraction of the trace, we find that subtracting it fromthe corresponding nucleon matrix element significantlyimproves the signal due to the correlation between thetwo terms. For the vacuum expectation value (cid:104) | ¯ T ig | (cid:105) we find that it is also compatible with zero but subtract-ing it from the nucleon matrix element does not improvethe signal. Therefore, in this case, subtraction of the vacuum expectation value is not performed. The gluonicratios using the diagonal and non-diagonal elements ofEMT are shown in Fig. 8. For both cases the plateausvalues, obtained by fitting in the range t ins /a ∈ [3 , t s − t s , show convergence and agreement with theresults extracted using the summation method. Thus,we use the plateau values for t s (cid:38) (cid:104) x (cid:105) gB isdetermined from the matrix element of the off-diagonalelements of EMT. We find (cid:104) x (cid:105) g B = 0 . n S = 10 stout smearing steps. We note that, fordisconnected quantities, where effects from excited satesare significantly milder, in combination with the largerstatistical errors, two-state fits are not reliable and are1therefore not presented. n S › x fi g B › x fi gB ( ¯ T ) › x fi gB ( ¯ T ) FIG. 9: Bare results for (cid:104) x (cid:105) gB as a function of the number ofstout smearing steps. With red squares are results extractedusing Eq.(25) and with blue circles results extracted usingEq.(24). The open symbol shows the selected value given inEq.(34) with its associated error band. By performing the same analysis for different steps ofstout smearing we can investigate the dependence on thesmearing steps n S . In Fig. 9 we show the dependenceof the extracted value of (cid:104) x (cid:105) gB on the number of stoutsmearing steps when using both the diagonal and off-diagonal elements of EMT. As can be seen, the errors decrease as n S increases and the values converge when n S (cid:38)
8. This means that the renormalization functionsshould also converge for n S (cid:38)
8, since the renormalizedmatrix element should be independent of the stout smear-ing. Details about the renormalization will be providedin Sec. V. B. B at zero momentum transfer As already discussed in connection to Eq.(10), directaccess to B (0) is not possible due to the vanishing ofthe kinematical factor in front of B ( Q ). Therefore,one needs to compute B ( Q ) for finite Q and extrapo-late to Q = 0 using a fit Ansatz. In order to accomplishthis, one has to isolate from the other two GFFs appear-ing in the decomposition of Eq.(10). We thus need tocompute the three-point function of the one-derivativevector operator for finite momentum transfer using bothunpolarized and polarized projectors and for both diag-onal and off-diagonal elements of the traceless EMT.To isolate B ( Q ) from A ( Q ) and C ( Q ) one hasto first minimize χ = (cid:88) ρ,µ,ν (cid:88) (cid:126)p (cid:48) ,(cid:126)p ∈ Q (cid:20) G µν (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ) F ( Q ; t s , t ins ) − R µν (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ; t s , t ins ) w µν (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ; t s , t ins ) (cid:21) , (35)where R is the ratio of Eq. (16) and w its statistical error.The kinematical coefficients G are defined in Appendix A.The three form factors are the components of F ( Q ; t s , t ins ) = A ( Q ; t s , t ins ) B ( Q ; t s , t ins ) C ( Q ; t s , t ins ) . (36)The time dependence t s , t ins appears due to contributionsfrom excited states that will be analyzed using the meth-ods discussed in Sec. III C. In the following discussionwe suppress the time dependence for simplicity. As dis-cussed, one can extract the form factors by minimizingthe χ in Eq. (35) or alternatively show that it is equiv-alent to F = V † Σ − U † ˜ R (37)where˜ R µν (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ) ≡ [ w µν (Γ ρ , (cid:126)p (cid:48) , (cid:126)p )] − R µν (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ) , ˜ G µν (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ) ≡ [ w µν (Γ ρ , (cid:126)p (cid:48) , (cid:126)p )] − G µν (Γ ρ , (cid:126)p (cid:48) , (cid:126)p ) , and˜ G = U Σ V, (38) where we compute the Singular Value Decomposition(SVD) of ˜ G in the last line. U is a hermitian N × N matrix with N being the number of combinations of µ , ν , ρ and components of (cid:126)p (cid:48) , (cid:126)p that contribute to the same Q . V is a hermitian 3 × N (cid:29) N × G . In the following we use the latter approach since noexplicit minimization is needed.After extracting B B ( Q ; t s , t ins ) using the SVDmethod we investigate its dependence on t s and t ins fol-lowing the same procedure as for the average momen-tum fraction. In Fig. 10 we present the analysis forthe connected contribution to the bare isoscalar B u + + d + B for two representative values of the momentum trans-fer. As for the case of the connected contributions to theisoscalar average momentum fraction, we observe largeeffects from excited states. As t s increases the valuechanges from negative to positive. Two-state fits yieldconsistent results as we increase t low s . The value extractedfrom the two-state fit for t low s /a = 8 is in agreement withthe value determined using the summation method for t low s ≥ .
12 fm and we thus select it as our final value.2 t ins t s /2 [fm]0.100.050.000.050.10 B u + + d + ; B ( C o nn . ) t s [fm] 0.6 0.8 1.0 1.2 1.4 t lows [fm]0.8 0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.8 t ins t s /2 [fm]0.100.050.000.050.100.15 B u + + d + ; B ( C o nn . ) t s [fm] 0.6 0.8 1.0 1.2 1.4 t lows [fm] FIG. 10: Excited state analysis for determining the connected contribution to bare B u + + d + B ( Q ) for Q = 0 .
32 GeV (top) and Q = 0 .
60 GeV (bottom). The notation is the same as that in Fig. 4. Q [GeV ] B u + + d + ; B ( Q ) ( C o nn . ) FIG. 11: The connected bare B u + + d + B ( Q ) as a function of Q . The red band is the result of a constant fit while the blueof a linear fit with fit ranges designated by the range of theband. The open symbol is the extrapolated value at Q = 0. In Fig. 11 we show results for the connected bare GFF B u + + d + B ( Q ) as a function of Q up to 1 GeV . As canbe seen, the Q behavior is relatively flat for small valuesof Q . In order to extrapolate to zero momentum we usea dipole form B B ( Q ) = B B (0)(1 + Q M ) (39)supported by the quark-soliton model in the large N c -limit for Q ≤ [70], where M the mass ofthe dipole. Since we are interested in fitting the small Q -dependence where the GFF is flat, one can expandEq. (39) as B B ( Q ) = B B (0) (cid:18) − Q M (cid:19) (40)for Q M (cid:28)
1. In the zeroth approximation, Eq. (40) yieldsa constant and in the first approximation a linear func-tion of Q M . In Fig. 11 we show the fits to both a constantand linear forms up to Q = 0 . . As can be seen,both constant and linear fits are in agreement, confirm-ing that the Q /M is negligible. We take as our valuethe one extracted from the constant fit, to obtain for theconnected isoscalar B u + + d + B (0) = 0 . . (41)Following the same procedure we extract the dis-connected light and strange quark contributions to B B ( Q ). In Fig. 12 we show the results as a functionof Q . Although the results are rather noisy for bothquark flavors we can fit the Q -dependence to the formof Eq. (39). We find for disconnected isoscalar B u + + d + B (0) = − . B s + B (0) = − . . (43)The gluon GFF B g B ( Q ) is shown in Fig. 13 as afunction of Q and the value extracted using a constant3 Q [GeV ] B u + + d + ; B ( Q )( D i s c . ) Q [GeV ] B s + ; B ( Q ) FIG. 12: Disconnected contributions to B u + + d + B ( Q ) in thetop panel and strange contributions to B s + B ( Q ) in the bot-tom panel with notation as in Fig. 11. Results are extractedusing the plateau method for a small separation t s /a = 6. Q [GeV ] B g ; B ( Q ) FIG. 13: Gluon contribution to B g B with notation as inFig. 11. Results are extracted using the plateau method fora separation t s /a = 10. fit is B g B (0) = − . . (44)In Table IV we tabulate our values on the bare resultsfor (cid:104) x (cid:105) and B (0). TABLE IV: Bare results for the momentum fraction (cid:104) x (cid:105) = A (0) and B (0) for the isovector, isoscalar connected, isoscalardisconnected, strange, charm, and gluon contributions. For details on the extraction of the isovector B u + − d + B see Ref. [51]. u + − d + u + + d + (Connected) u + + d + (Disconnected) s + c + g (cid:104) x (cid:105) B (0) 0.130(36) 0.018(25) -0.038(38) -0.017(18) - -0.049(40) V. RENORMALIZATION
One of the main ingredients of this work is the renor-malization of the quark and gluon parts of the EMTof Eqs. (4) and (3). The renormalization of each partrequires a dedicated calculation, and in this section weclassify them in multiplicative renormalization functions( Z qq , Z gg ) and mixing coefficients ( Z qg , Z gq ). The lat-ter are needed to disentangle the quark and gluon mo-mentum fractions from the bare matrix element of theoperators of Eqs. (3) and (4). Therefore, a 2 × (cid:104) x (cid:105) q + R = Z qq (cid:104) x (cid:105) q + B + Z qg (cid:104) x (cid:105) gB , (45) (cid:104) x (cid:105) gR = Z gg (cid:104) x (cid:105) gB + Z gq (cid:104) x (cid:105) q + B . (46)In the above equations, (cid:104) x (cid:105) q + is understood to be the fla-vor singlet combination that sums the up, down, strangeand charm quark contributions. The subscript R ( B )represents the renormalized (bare) matrix elements. Wenote that a complete calculation of the 2 × Z gg and Z qq , since the mixing coefficients Z gq and Z qg are small as shown in Ref. [33] using the one-loop lattice perturbation theory. Furthermore, note thatgluon and quark EMT operators also mix with gauge-variant operators (BRS-variations and operators vanish-ing by the equations of motion), which have zero nucleonmatrix elements and are not considered.In the following subsections we present the non-perturbative renormalization of the quark EMT, the non-perturbative renormalization of the gluon EMT and ourestimates for the mixing coefficients as extracted froma calculation within one-loop lattice perturbation the-ory [33]. A. Quark EMT renormalization
The quark EMT is renormalized non-perturbativelyusing an analysis within the Rome-Southampton scheme(RI (cid:48) scheme) [71]. This is a very convenient prescriptionfor non-perturbative calculations, and is obtained by ap-plying the conditions Z q = 112 Tr (cid:2) ( S L ( p )) − S Born ( p ) (cid:3)(cid:12)(cid:12)(cid:12) p = µ , (47) Z − q Z qq
112 Tr (cid:104) Γ L qEMT ( p ) (cid:0) Γ BornqEMT (cid:1) − ( p ) (cid:105)(cid:12)(cid:12)(cid:12) p = µ = 1 . (48)The trace in the above conditions is taken over spinand color indices. The momentum p of the vertex func-tions is set to the RI (cid:48) scale, µ . Note that the vertexfunction Γ L qEMT is amputated in the above condition.Also, S Born (Γ BornqEMT ) is the tree-level value of the quarkpropagator (quark operator). We employ the momentumsource method introduced in Ref. [72], which offers highstatistical accuracy using a small number of gauge con-figurations as demonstrated for twisted mass fermions inRefs. [73–75]. Discretization effects and other systematicuncertainties in the renormalization functions (Z-factors)can be amplified based on the choice of the momentum.A way around this problem is the use of momenta withequal spatial components. The temporal component isthen chosen such that the ratio P ≡ (cid:80) i p i / ( (cid:80) i p i ) isless than 0.3 [76]. Such a ratio is relevant to Lorentznon-invariant contributions present in perturbative cal-culations of Green’s functions beyond leading order in a [77]. Therefore, a large value of P a effects in the non-perturbative estimates too. Themomenta employed in this work for the quark EMT areof the form ap ≡ π (cid:18) n t T /a , n x L/a , n x L/a , n x L/a (cid:19) , n t ∈ [2 , , n x ∈ [2 , , (49) taking all combinations of n t and n x that satisfy P < . ≤ ( a p ) ≤ T and L are the temporal andspatial extent of the lattice, and correspond to T /a = 48,
L/a = 24 for the N f = 4 ensembles that are generatedspecifically for the renormalization program at the samecoupling constant as the cB211.072.64 ensemble.An important aspect of our renormalization programis the improvement of the non-perturbative estimates bysubtracting finite- a effects [75, 78], calculated to one-loopin lattice perturbation theory and to all orders in thelattice spacing, O ( g a ∞ ). Note that the dimensionlessquantity appearing in the perturbative expressions is ap (for massless fermions). Z qq has two components depending on the indices ofthe operator defined in Eq. (4). Z qq corresponds to thequark EMT operator with µ = ν , while Z qq to µ (cid:54) = ν . Z qq and Z qq renormalize the bare matrix elements ofthe quark EMT operator obtained with the same con-straints on the external indices. Thus, in our work weuse Z qq for the connected contribution and Z qq for thedisconnected ones, as described in Sec. IV A.For the proper extraction of Z qq we use five N f = 4 en-sembles at different pion masses reproducing a β value of1.778 to match the N f = 2 + 1 + 1 ensemble on which thebare matrix elements have been calculated. The N f = 4ensembles correspond to a pion mass that ranges between350 and 520 MeV, allowing one to take the chiral limit.More details on the N f = 4 ensembles can be found inRef. [51]. The chiral extrapolation is performed using aquadratic fit with respect to the pion mass of the form Z RI (cid:48) ( µ ) + ¯ z RI (cid:48) ( µ ) · m π , (50)where Z RI (cid:48) and ¯ z RI (cid:48) depend on the scheme and the scale.The pion mass dependence of Z RI (cid:48) qq is found to be verymild, as demonstrated in Fig. 14 where we show the datafrom the five ensembles for a representative renormaliza-tion scale ( a µ ) = 2. Same conclusions hold for Z qq . Z qq FIG. 14: Z RI (cid:48) qq at a scale ( a µ ) = 2, as a function of the pionmass squared in lattice units. The dashed lines correspondto the chiral extrapolation using Eq. (50) leading to a valueshown with open circle in the chiral limit ( Z RI (cid:48) qq ). Once the chiral extrapolation is performed, we applythe subtraction of artifacts calculated in one-loop lattice5perturbation theory. This subtraction procedure leads toimproved estimates, as it significantly reduces discretiza-tion effects. The next step is the conversion to the MS-scheme, which is commonly used to compare to exper-imental and phenomenological values. The conversionprocedure is applied on the Z-factors obtained on eachinitial RI (cid:48) scale ( a µ ), with a simultaneous evolution toa MS scale, chosen to be µ =2 GeV. Assuming absenseof mixing, the conversion and evolution uses the inter-mediate Renormalization Group Invariant (RGI) scheme,which is scale independent and relates the Z-factors be-tween the two schemes: Z RGI qq = Z RI (cid:48) qq ( µ ) ∆ Z RI (cid:48) qq ( µ )= Z MS qq (2 GeV) ∆ Z MS qq (2 GeV) . (51)Therefore, the appropriate factor to multiply Z RI (cid:48) qq is [96] C RI (cid:48) , MS qq ( µ , ≡ Z MS qq (2 GeV) Z RI (cid:48) qq ( µ ) = ∆ Z RI (cid:48) qq ( µ )∆ Z MS qq (2 GeV) . (52)The quantity ∆ Z S qq ( µ ) is expressed in terms of the β -function and the anomalous dimension γ Sqq ≡ γ S of theoperator∆ Z S qq ( µ ) = (cid:32) β g S ( µ ) π (cid:33) − γ β × exp (cid:40)(cid:90) g S ( µ )0 d g (cid:48) (cid:18) γ S ( g (cid:48) ) β S ( g (cid:48) ) + γ β g (cid:48) (cid:19)(cid:41) , (53)and may be expanded to all orders of the coupling con-stant. The expression for the quark EMT operator isknown to three-loops in perturbation theory and can befound in Ref. [75] and references therein.The conversion and evolution is followed by a fitto eliminate the residual dependence on aµ using theAnsatz Z qq ( a µ ) = Z qq + z qq · ( a µ ) . (54)For both Z qq and Z qq we distinguish between the sin-glet and the non-singlet case, which is necessary for theproper renormalization and the flavor decomposition pre-sented in Sec. VI. Their difference is known to be verysmall as it first appears in two-loop perturbation the-ory [79]. Here, we calculate also the singlet renormal-ization function non-perturbatively, which requires bothconnected and disconnected contributions to the vertexfunctions.In the upper panel of Fig. 15 we plot the nonsingletvalues of Z MS qq (at 2 GeV) with and without the subtrac-tion of the corresponding O ( g a ∞ ) contributions. Wealso include the singlet Z ( s ) , MS qq , after subtraction of the O ( g a ∞ ) terms. We find that the singlet and nonsingletrenormalization functions are compatible within uncer-tainties, with the singlet being more noisy, due to the inclusion of the disconnected contributions. In the lowerpanel of Fig. 15 we show the corresponding quantities for Z MS qq . While the singlet one has large uncertainties inthis case too, it is smaller than the statistical errors of Z ( s ) , MS qq . In both cases we have subtracted the vacuumexpectation value.We find that the subtraction procedure improves sig-nificantly the data, leading to smaller dependence on theinitial scale ( a µ ) . As can be seen from the plot, the O ( g a ∞ )-terms capture a large part of the discretiza-tion effects. The subtraction of finite- a terms from thenon-perturbative estimates of Z MS qq (as well as Z MS qq ) re-duces the slope with respect to ( a µ ) , between mo-menta with the same n x value and different n t . Asan example, let us consider the class of momenta with n x = 3 (( a µ ) ∈ [2 − . z unsub qq = 0 . z sub qq = 0 . a µ ) considered in this work. Z MSqq1 (Unsub. ) Z
MSqq1 (Sub. ) Z (s)MSqq1 (Sub. ) ) Z MSqq2 (Unsub. ) Z
MSqq2 (Sub. ) Z (s)MSqq2 (Sub. )
FIG. 15: Z MS qq (upper) and Z MS qq (bottom) as a function of theinitial RI (cid:48) scale ( a µ ) . The purely non-perturbative dataare shown with green crosses, and the improved estimatesafter the subtraction of O ( g a ∞ )-terms are shown with redcircles. The blue squares show results of the singlet case aftersubstraction of lattice artifacts. The dashed lines show thefit using Eq. (54), and the extrapolated values with an opensymbol. The final estimates for Z qq and Z qq are determinedusing the fit interval ( a µ ) (cid:15) [2 − Z MS qq = 1 . , (55) Z MS qq = 1 . . (56)6The numbers in the first and second parenthesis corre-spond to the statistical and systematic uncertainties, re-spectively. The source of systematic error is related tothe ( a µ ) → aµ ) = 2) and higher (( aµ ) = 7) fit rangesand taking the largest deviation as the systematic error.We emphasize that the procedure of improving the Z-factors utilizing lattice perturbation theory, has impor-tant implication on the spin and momentum decomposi-tion: use of the unimproved Z-factors would underesti-mate both the intrinsic spin and the quark momentumfraction by 5%. B. Gluon EMT renormalization
Similarly to the case of the quark EMT, we renormal-ize the gluon EMT non-perturbatively. This is a cru-cial improvement compared to our previous work [26, 33]in which we used Z gg from one-loop perturbation the-ory. The renormalization condition for Z gg involves thegluon-field renormalization function Z g , which in the RIscheme reads Z g = N c −
12 3 / ˆ p (cid:80) ρ (cid:104) Tr[ A ρ ( p ) A ρ ( − p )] (cid:105) (cid:12)(cid:12)(cid:12) p = µ , (57) Z gg = N c − Z g (2ˆ p ˆ p i ) / (ˆ p ) (cid:104) Tr[ A ρ ( p ) ¯ T i ; g A ρ ( − p )] (cid:105) (cid:12)(cid:12)(cid:12)(cid:12) ρ (cid:54) = i (cid:54) =4 ,pρ =0 ,pi = µi . (58)In the above equations N c is the number of colors and a ˆ p j = 2 sin ( a p j / N c − / A ρ ( p ) = a (cid:88) x e + ip · ( x +ˆ ρ/ (cid:34)(cid:32) U ρ ( x ) − U † ρ ( x )2 iag (cid:33) − (Trace) (cid:35) (59)and the gluon propagator in momentum space is given as (cid:104) A ρ ( p ) A ρ ( − p ) (cid:105) in the ρ direction. The numerator 3 / ˆ p ofEq. (57) is the tree-level expressions for the gluon propa-gator in the Landau gauge, in which the Lorentz indiceshave been set equal to each other and are summed over.Similarly, (2ˆ p ˆ p i ) / (ˆ p ) in Eq. (58) is the non-amputatedtree-level value corresponding to the gluon EMT in theLandau gauge. In this study we focus on the ¯ T ig case,since as shown in Fig. 9 it is significantly more precise.This setup justifies the presence of Z − g in Eq. (58),instead of Z +1 g . For simplicity in the notation, the de-pendence of Z g and Z gg on the RI scale µ is implied.Unlike the case of Z qq , here we use non-amputated vertexfunctions for the gluon EMT. Such a choice is desirable,as the 4 × Z g given in Eq. (57) is convenient, asthere is a sum over the Lorentz indices of the gluon fields.While a similar condition could be imposed on Z gg , we do not sum over ρ in Eq. (58). Instead we choose the index ρ to be different from the Lorentz indices of the operator(4 and i ). This has the advantage that any mixing withother gluon operators [80] vanishes automatically, at leastto one-loop level.We also explore an alternative definition of Z gg as pro-posed in Ref. [34]. In such condition for Z g there is nosummation over ρ , which is set equal to the ρ index of theexternal gluon fields in Eq. (58), and has the constraint p ρ = 0. Therefore, it is convenient to eliminate Z g fromEq. (58), obtaining Z gg = 2ˆ p ˆ p i (cid:104) A ρ ( p ) A ρ ( − p ) (cid:105) ˆ p (cid:104) A ρ ( p ) T ig A ρ ( − p ) (cid:105) (cid:12)(cid:12)(cid:12)(cid:12) ρ (cid:54) = i (cid:54) =4 ,p ρ =0 ,p i = µ i . (60)One major difference in the calculation of Z gg ascompared to Z qq is the need to reduce the high noise-to-signal ratio appearing in the calculation of gluonicquantities. To this end, some equivalent renormalizationprescriptions have been proposed to reduce the statis-tical uncertainties. In the discussion that follows wewill investigate three methods for the extraction of Z gg .Note that for zero stout steps, n S = 0, all these methodsreduce to the same equation. Method 1:
Application of stout smearing only onthe operator ¯ T ig in Eq. (60), while the external gluonfields remain unsmeared. Method 2:
Application of stout smearing in boththe operator and the external gluon fields of Eq. (60) assuggested in Ref. [32]. Since our action is not smearedone would need to apply reweighting in the calculationof both the Z-factors and matrix elements. We followRef. [32] and assume that its effect is negligible on therenormalization function.
Method 3:
A generalization of Method 1 as suggestedin Ref. [81], in which we multiply Eq. (60) by the ratio R (( aµ ) ) ≡ f (( aµ ) ) f (( aµ ) → , (61)where f (( aµ ) ) = (cid:104) Tr[ A sρ ( p ) A sρ ( − p )] (cid:105)(cid:104) Tr[ A ρ ( p ) A ρ ( − p )] (cid:105) (cid:12)(cid:12)(cid:12) p = µ . (62)Presence of an index s implies stout smearing with n S steps. The multiplication of Eq. (60) by Eq. (61) leads tothe same Z gg in the ( a µ ) → R (( aµ ) ) → T ig and the gluon fields.The vertex functions entering Eq. (60) are calculatedon one N f = 4 ensemble with β = 1 .
778 and volume12 ×
24 with a pion mass of 350 MeV. While the Z-factors7are defined in the massless limit, the use of a single en-semble is sufficient given the negligible pion mass depen-dence observed for the quark case shown in Fig. 14, wherethe results are very precise. Focusing on a single ensem-ble allows one to reach higher statistics. As a purely glu-onic quantity it is susceptible to large gauge fluctuationsand therefore about 31,000 configurations are analyzedto reduce the noise.In contrast to the quark case, the momenta for the ver-tex functions cannot have the same spatial components,due to the constraint p ρ = 0 in the renormalization con-dition. Consequently, the value of P Z gg for momenta satis-fying P < .
4, and within the range 1 ≤ ( a µ ) ≤ C RI , MS gg ( µ , ¯ µ ) = 1+ g π (cid:18) N c − N f (cid:18)
53 + log( ¯ µ µ ) (cid:19)(cid:19) , (63)which must multiply Z RI gg . In the expression above, ¯ µ isthe renormalization scale in the MS scheme and µ inthe RI scheme, as defined previously. This conversionfactor is consistent with the expression of Ref. [81], withthe only difference a global minus sign in the one-loopexpression, due to a different definition.We obtain Z gg in the MS scheme evolving the valuesfrom every scale ( a µ ), and therefore, we must eliminateany residual dependence on ( a µ ) by taking the limit( a µ ) →
0. While such a fit is linear in Z qq wheredemocratic momenta can be used, here the residual de-pendence on ( a µ ) may be polynomial. We identifytwo sources for this behavior: i. Truncation effects in theconversion factor C gg , which is only known to one-looporder; ii. finite- a effects due to p ρ = 0, making it unreli-able to go to high ( a µ ) values. The latter effect can bereduced by a similar subtraction of finite- a effects as inthe case of Z qq . Since these have not yet been calculatedwe cannot make the subtraction in the current data.
1. Method 1
The most straightforward and conceptually concretemethod to extract Z gg is to use Eq. (60) with stoutsmearing applied only on the gauge links of the oper-ator ¯ T ig . In Fig. 16 we show ( Z MS gg ) − as a function ofthe initial scale squared for various smearing steps. Aswe increase the number of stout steps, the ( a µ ) → a µ ) dependence since more smearing alters thediscretization effects between the numerator and denom-inator in Eq. (60). For the examples shown in Fig. 16 we Z ( M S ) gg ( n S =3 ) Z ( M S ) gg ( n S =5 ) ) Z ( M S ) gg ( n S =10 ) FIG. 16: The inverse Z gg in the MS as a function of theinitial RI scale ( aµ ) , using method 1. From top to bottomthe plots show three cases for the number of stout smearingsteps, namely n S = (3 , ,
10) with (filled green triangles, filledblue squares, filled red circles). The extrapolated values in the( aµ ) → use a polynomial of second degree with respect to ( a µ ) for three and five steps, and third degree for 10 steps ofstout smearing. We observe that the value of ( Z MS gg ) − increases ( Z MS gg decreases) with the stout steps, which isexpected, as the value of the bare matrix elements in-creases with the stout steps. As will be discussed later,we find that the renormalized matrix element is indepen-dent of the number of stout steps (see Fig. 21).
2. Method 2
The difference between method 1 and method 2 is theuse of stout smearing on the gauge links used to constructthe gluon fields entering Eq. (60). As already mentioned,this would need reweighting. This was assumed to benegligible in Ref. [32] and we also neglect it here. InFig. 17 we show Z MS gg from method 2 for selected num-ber of stout steps including zero steps, the same for allthree methods. Without smearing there is no noticeabledependence on the scale ( a µ ) allowing us to fit to aconstant while as increasing the number of steps the de-pendence becomes linear. We note that smearing alsothe gluon field, provides a better correlation with the op-erator for higher momenta allowing us to investigate upto ( a µ ) = 7. It is worth mentioning that while thereis a big jump on the extrapolated value between n S = 0and n S = 5, between n S = 5 and n S = 10 the differenceis relatively small.8 Z M S gg ( n S =0 ) Z M S gg ( n S =5 ) ) Z M S gg ( n S =10 ) FIG. 17: Z MS gg using method 2 with the same notation asin Fig. 16. From top to bottom we show results for n S =(0 , ,
3. Method 3
This method is more involved since one has to com-pute first Eq. (62) and extrapolate to ( aµ ) →
0. InFig. 18 we show the ratio of smeared to unsmeared gluonpropagators. Note that this ratio does not alter the con-version factor for Z gg . As we increase the number ofsmearing steps the discretization effects between the nu-merator and the denominator change, not canceling inthe ratio leading to a more curved behavior. We fit theresults up to a forth order polynomial since the gluonpropagators alone are very precise.In Fig. 19 we show the Z MS gg when the Eq. (61) multipliesEq. (60). The resulting behavior is fitted to a constantsince most of the systematics are cancelled when multi-plying with Eq. (61).It is interesting to compare the final extrapolated val-ues of Z MS gg among the three methods. The results fromeach method should agree since they renormalize thesame bare matrix elements. Such a comparison will givean indication of additional discretization effects, whichmight remain after the ( a µ ) → n S . We find thatall three methods are overall compatible as a function ofnumber of stout smearing steps. It is worth mentioningthat the Z MS gg has a strong dependence on n S going fromzero steps up to five steps, whereas increasing further thesteps the Z-factor does not change significantly. In Fig. 9 ) f (( a ) ) n S = 1n S = 3n S = 5 FIG. 18: The ratio of Eq. (62) as a function of the initialRI scale ( aµ ) for stout smearing steps, n S = (1 , ,
5) with(green triangles, blue squares, red circles). The extrapolatedvalues in the ( aµ ) → Z M S gg ( n S =1 ) Z M S gg ( n S =3 ) ) Z M S gg ( n S =5 ) FIG. 19: Z MS gg using method 3 with notation as in Fig. 16with number of stout smearing steps n S = (1 , ,
5) from topto bottom. we demonstrated that the bare matrix element shows anincrease from zero steps up to 12 steps, albeit large errorsfor small number of steps, tending to converge after 10steps.One important consistency check for the calculationof Z gg is the comparison of the renormalized matrix ele-ments between different methods, which is demonstratedin Fig. 21. For simplicity, we neglect the mixing for thisdiscussion. Such mixing is found to be very small (seeSubsec. V C), and thus, does not alter the main conclu-sions of Fig. 21. The multiplication of the bare matrix9 S Z M S gg Method 1Method 2Method 3
FIG. 20: Z MS gg as a function of the number of stout smearingsteps, for the three methods described in the text. For n S = 0all three methods reduce to the same method given with theblack cross. element (cid:104) x (cid:105) g by Z MS gg is shown in Fig. 21 for the threemethods investigated. As can be seen, the three meth-ods yield compatible results for all stout steps. While thestout smearing does not alter the values of the renormal-ized matrix element, it has the advantage that it reducesthe gauge fluctuations. n S › x fi g R Method 1Method 2Method 3
FIG. 21: The renormalized gluon average momentum fractioncomputed from ¯ T i ; g using the three methods described in themain text. The selected method and value is given with theopen symbol and the corresponding green horizontal band. The chosen value for Z gg in the MS scheme at 2 GeV isobtained at 10 stout steps using method 1 as shown inFig. 21. Z MS gg = 1 . , (64)which is a conservative choice as it has the largest sta-tistical uncertainty compared to the other methods (seeFig. 20). A systematic has been added by varying thehighest point in the polynomial fit from 4 to 3.5 of theinitial RI scale ( aµ ) . The value above will be used inthe 2 × C. Mixing between fermion and gluon operators
The renormalization of the quark and gluon EMT ismore complicated as compared to other operators stud- ied within hadron structure (e.g., intrinsic spin). This isdue to their mixing, resulting into a 2 × Z qq , Z qg , Z gg , Z gq of the mixing matrixcan be obtained within lattice perturbation theory, fol-lowing the procedure of our previous work on the gluonEMT [33]. In particular, to one loop level, one needs tocalculate the diagrams of Fig. 22. (a)(b)(c)(d)FIG. 22: One-loop diagrams contributing to Z qq , Z qg , Z gg , Z gq . Diagrams (a) and (b) have an in-sertion of the quark operator of Eq. (4) (filled square) andexternal quark (solid lines) and gluon (wavy lines) fields,respectively. Diagrams (c) and (d) have an insertion of thegluon operator of Eq. (3) (filled circle) and external gluonand quark fields, respectively. Here we are interested in the extraction of Z gq and Z qg from our perturbative calculation, as Z gg and Z qq arecomputed non-perturbatively. In the calculation withinperturbation theory we use up to two steps of stoutsmearing for the gluon EMT. This limitation is posedby the fast increase of algebraic expressions (millions ofterms), for higher number of stout steps. We find that thepolynomial nature of the perturbative renormalization0functions with respect to the stout parameter, leads to aconvergence at a small number of stout steps. This hasbeen confirmed in other calculations with stout smear-ing [33, 83]. Using the lattice spacing and coupling con-stant of the ensemble under study we extract the mixingcoefficients: Z MS qg = 0 . , (65) Z MS qg = 0 . , (66) Z MS gq = − . . (67) VI. RESULTS
In this section we give the renormalized matrix ele-ments, by combining the bare matrix elements extractedin Sec. IV and the renormalization factors in Sec. V yield-ing our physical results. The renormalized results areobtained from the expressions X q + R = Z qq X q + B + δZ qq N f (cid:88) q = u,d,s,c X q + B + Z qg N f X gB (68) X gR = Z gg X gB + Z gq (cid:88) q = u,d,s,c X q + B (69)where X = (cid:104) x (cid:105) , J , and δZ qq the difference between sin-glet and non-singlet Z qq and N f = 4 since we have fourflavors in the sea. In order to fully decompose the quarkflavors we use the corresponding isovector results fromRefs. [25, 51], which are also given in Table VI for com-pleteness.In Fig. 23 we show our results for the proton aver-age momentum fraction for the up, down, strange andcharm quarks, for the gluons as well as their sum. Theup quark is the largest quark contribution, namely about35%, twice bigger than the down quark. The strangequark contributes significantly smaller, namely about 5%and the charm is restricted to about 2%. The gluonhas a significant contribution of about 45%. Summingall the contributions results to (cid:80) q = u,d,s,c (cid:104) x (cid:105) q + R + (cid:104) x (cid:105) gR =104 . . u + d + s + c + X q + = u, d, s, c g Total0.00.20.40.60.81.01.2 › x fi q + , g p . ( . ) % . ( . ) % . ( . ) % . ( . ) % . ( . ) % . ( . ) % . ( . ) % FIG. 23: The decomposition of the proton average momen-tum fraction (cid:104) x (cid:105) . We show the contribution of the up (redbar), down (green bar), strange (blue bar), charm (orangebar), quarks and their sum (purple bar), the gluon (cyan bar)and the total sum (grey bar). Note that what is shown is thecontribution of both the quarks and antiquarks ( q + = q + ¯ q ).Whenever two overlapping bars appear the inner bar denotesthe purely connected contribution while the outer one is thetotal contribution which includes disconnected taking into ac-count also the mixing. The error bars for the former areomitted while for the latter are shown explicitly on the bars.The percentages written in the figure are for the total con-tribution. The dashed horizontal line is the momentum sum.Results are given in MS scheme at 2 GeV. The (cid:80) q = u,d,s B q + (0) + B g (0) is expected to vanish torespect the momentum and spin sums, as pointed out byEq. (11). We find for the renormalized values that (cid:88) q = u,d,s B q + ,R (0) + B g ,R (0) = − . ∆Σ q + = g q + A . These are taken from Ref. [25] and in-cluded in Table V, for easy reference. The up quark has alarge contribution, up to about 85% of the proton intrin-sic spin. The down quark contributes about half com-pared to the up and with opposite sign. The strange andcharm quarks also contribute negatively with the latterbeing about five times smaller than the former giving a1% contribution. The total ∆Σ q + is in agreement withthe upper bound from COMPASS [84].Having both the quark angular momentum and thequark intrinsic spin allows us to extract the orbital an-gular momentum using Eq. (12). For a direct calcula-tion using TMDs see Ref. [85]. Our results are shownin Fig. 26. The orbital angular momentum of the up1 u + d + s + c + X q + = u, d, s, c g Total0.00.10.20.30.40.50.6 J q + , g p . ( . ) % . ( . ) % . ( . ) % . ( . ) % . ( . ) % . ( . ) % . ( . ) % FIG. 24: The decomposition of the proton spin J . The colornotation of the bars is as in Fig. 23. The dashed horizontal lineindicates the observed proton spin value and the percentageis given relative to the total proton spin. Results are given inMS scheme at 2 GeV. u + d + s + c + Total0.30.20.10.00.10.20.30.40.50.6 q + p . ( . ) % - . ( . ) % - . ( . ) % - . ( . ) % . ( . ) % FIG. 25: Results for the intrinsic quark spin ∆Σ contri-butions to the proton spin decomposed into up (red bar),down (green bar), strange (blue bar), and charm (orange bar).The total contribution of the four flavor is also shown (greybar) [25]. The dashed horizontal line is the observed protonspin and the percent numbers are given relative to it. Resultsin MS scheme at 2 GeV. quark is negative reducing the total angular momentumcontribution of the up quark to the proton spin. Thecontribution of the down quark to the orbital angularmomentum is positive almost canceling the negative in-trinsic spin contribution resulting to a relatively smallpositive contribution to the spin of the proton.Our final values for each quark flavor and gluon contri-bution to the intrinsic spin, angular momentum, orbitalangular momentum and momentum fraction of the pro-ton are summarized in Table V. u + d + s + c + Total0.30.20.10.00.10.20.30.40.50.6 L q + p - . ( . ) % . ( . ) % . ( . ) % . ( . ) % . ( . ) % FIG. 26: Orbital angular momentum L contributions to theproton spin. The color notation is as in Fig. 25. The dashedhorizontal line denotes the observed proton spin and the per-centage is given relative to the it. Results are given in MS at2 GeV.TABLE V: Results for the proton for the average momentumfraction (cid:104) x (cid:105) , the intrinsic quark spin ∆Σ [25], the total an-gular momentum J and the orbital angular momentum L inthe MS scheme at 2 GeV. Results are given separately for theup ( u + ), down ( d + ), strange ( s + ), charm ( c + ) and for gluons( g ) where for the quarks, results include the antiquarks con-tribution. The sum over quarks and gluons is also given asTot. (cid:104) x (cid:105) J ∆Σ L u + d + s + c + g (cid:104) x (cid:105) , J and ∆Σ. (cid:104) x (cid:105) J ∆Σ u + − d + The results given in Table V and VI are obtained us-ing one ensemble of twisted mass fermions. Therefore, itis not possible to quantitatively determine finite latticespacing and volume systematics. However, in Ref. [86]several N f = 2 twisted mass fermion ensembles wereanalyzed with pion masses in the range of 260 MeVto 470 MeV and lattice spacings a =0.089, 0.070 and0.056 fm as well as for two different volumes. A contin-uum extrapolation at a given value of the pion mass wasperformed. We found negligible O ( a )-terms yielding aflat continuum extrapolation. Therefore, we expect thatcut-off effects will be small also for our current ensemble.2 VII. COMPARISON WITH OTHER STUDIES
In order to evaluate the contribution of each quark fla-vor to the proton spin and momentum one needs to com-pute the quark-disconnected diagrams as done here. Theevaluation of these contributions is much more challeng-ing as compared to the connected ones, in particular atthe physical point. This is the main reason that most lat-tice QCD studies to date have mostly computed isovectorquantities for which the aforementioned diagrams cancel.For instance, in the case of the axial charge, which isan isovector quantity, there are numerous studies [87],whereas for the individual quark flavor axial charges g q + A ≡ ∆Σ q + results computed directly at the physicalpoint are still scarce. In order to make a comparisonwith other lattice QCD studies, we include results ob-tained using a chiral extrapolation. We limit ourselvesto comparing results that were obtained having at leastone ensemble with close to physical pion mass, meaningbelow 180 MeV. Although such a chiral extrapolationmay introduce uncontrolled systematics that are absentfrom the results reported here, it allows for a comparisonwith published lattice QCD results on these quantities.We begin with ∆Σ q + and consider the following lat-tice QCD studies:i) The χ QCD collaboration analyzed three N f = 2+1gauge ensembles of domain-wall fermion (DWF)generated by the RBC/UKQCD collaboration withpion masses 171 ,
302 and 337 MeV and lattice spacings of 0 . , .
111 and 0 .
083 fm. They usedoverlap fermions in the valence sector. They per-formed a combined fit in order to extrapolate tothe physical pion mass, the continuum and infinitevolume limits [29].ii) The PNDME collaboration analyzed several N f =2 + 1 + 1 gauge ensembles of highly Improved Stag-gered Quarks (HISQ) generated by the MILC Col-laboration. They used Wilson clover fermions inthe valence sector. For the connected contribu-tions they analyzed eleven ensembles with m π (cid:39) , ,
135 MeV and lattice spacings a (cid:39) . , . , . , .
06 fm. The disconnected contri-butions were computed on a subset of these ensem-bles. The strange quark contributions were com-puted on seven ensembles using all lattice spacingsbut only one physical point ensemble was analyzed;the light disconnected were computed on six ensem-bles for two values of m π = (315 , N f = 2 ensem-ble of twisted mass fermion with m π = 130 MeV, a = 0 .
094 fm and Lm π = 3 [27].
12 u + This WorkETMC '17 QCD '18PNDME '18 NNPDFpol1.1 DSSV08 JAM15 JAM17
12 d +
12 s +
12 c + FIG. 27: Results for ∆Σ q + . From left to right: for u + , d + , s + and c + quarks. Red, green, blue, and orange denote latticeQCD results, with filled symbols being results that are computed directly at the physical point and open symbols resultsthat were obtained after a chiral extrapolation. The inner error bar is the statistical error and the outer one the total thatincludes systematic errors. In particular, red circles show the results using the cB211.072.64 ensemble of this work and reportedin Ref. [25] with the associated error band. Green squares show ETM Collaboration results from Ref. [27]; blue upwardspointing triangles show results from χ QCD [29]; and orange left-pointing triangles from PNDME [28]. Results from globalfits of polarized PDFs are shown with black symbols and right triangles, pentagons, diamonds, and rhombus being fromNNPDFpol.1 [17], DSSV08 [11], JAM15 [88], JAM17 [89], respectively.
In Fig. 27 we show a comparison of our results on the intrinsic spin ∆Σ q + to the aforementioned lattice QCD3studies. As can be seen, there is an agreement amongdifferent lattice QCD analyses. In addition, we compareto the results extracted from global-fit analysis of polar-ized parton distribution. The values from the analysis ofthe cB211.072.64 ensemble of this work for the up anddown quarks agree very well with the phenomenologi-cal extractions. We note that the precision reached iscomparable to that of the phenomenological values. Forthe strange quark contribution ∆Σ s + lattice QCD re-sults achieve a better accuracy than the results extractedfrom global fits and point to a smaller value as comparedto those from DSSV08 [11] and JAM15 [88]. Our resultsfor ∆Σ c + predict a non-zero value, showing small butnon-zero charm quark effects in the nucleon.For the comparison of the average momentum fractionof each quark flavor and the gluon we consider latticeQCD results from the following groups:i) The χ QCD collaboration using the same gauge en- sembles as described for the case of the intrin-sic quark spin. In addition, they included in theanalysis an ensemble with m π = 139 MeV and a = 0 .
114 fm [81]. Despite the fact that a phys-ical point ensemble is included, a chiral extrapola-tion is still performed in order to extract their valueat the physical point. In using more precise resultsfor heavier than physical pion mass ensembles theirresult at the physical point is weighted less in thefit. This procedure may yield better precision atthe physical point but it can also potentially intro-duce an unknown systematic error due to the chiralextrapolation.ii) The ETM collaboration using the same setup as forthe intrinsic quark spin. › x fi u + This WorkETMC '17 QCD '18QCD( -extr.) '18 NNPDF3.1 CT14 MMHT2014 ABMP2016 CJ15 HERAPDF2.0 › x fi d + › x fi s + › x fi g FIG. 28: From left to right we show results for the nucleon average momentum fraction (cid:104) x (cid:105) for the up, down, and strangequark flavors as well as for the gluon. Red circles are the results of this work with the associated error band, green squaresshow results from the ETM Collaboration [26] and upwards-pointing triangles are from the χ QCD Collaboration [81] withfilled symbols being the results obtained directly at the physical point and open symbols using a chiral extrapolation. Resultsfrom the global fit analyses are show in black left- right- pointing triangles, pentagons, diamonds, rhombus, and down-pointingtriangles from NNPDF3.1 [90], CT14 [91], MMHT2014 [92], ABMP2016 [93], CJ15 [94], HERAPDF2.0 [95], respectively.
In Fig. 28 we compare the results for the average mo-mentum fraction for each quark flavor. The results high-light the improvement achieved in the current analysis ascompared to the two previous direct determinations usingphysical point ensembles by the ETM [26] and χ QCD [81]Collaborations. This is mostly due to the precise deter-mination of the quark loop (disconnected) contributionsusing our improved techniques. Additionally, our currentdetermination is in remarkable agreement with the phe-nomenological extractions resolving a long standing dis- crepancy between lattice QCD results and experimentaldeterminations. In Fig. 28 we also include a comparisonof the gluon momentum fraction where we only show lat-tice results with non-perturbative renormalization, thusexcluding the previous result from the ETM Collabora-tion [26]. There is agreement between the result of thisstudy and the one from the χ QCD Collaboration as wellas with the phenomenological determinations, which arevery precise compared to the current lattice QCD values.For the angular momentum and orbital angular momen- tum, the quark decomposition at the physical point has4 J u + J d + J s + J g FIG. 29: Results for the angular momentum J for each quark flavor and the gluon. With red circles are results from this workand with green squares results from our previous study [26]. The notation is as in Fig. 28. only been computed by the ETM Collaboration [26].We show the comparison between these two studies inFigs. 29 and 30, respectively. The results of this workhave improved accuracy for all quark flavors for this classof observables. Both J u + and L u + have smaller valueswhile the rest are in agreement with our previous study. This is due the fact that more sink-source time separa-tions are used reaching larger separations with more ac-curacy. This leads to a better extraction of the nucleonmatrix element. For J g the result of this study is the onlyone available with a non-perturbative renormalization atthe physical point. L u + L d + L s + FIG. 30: Results for the orbital angular momentum L for each quark flavor. The notation is as in Fig. 29. VIII. CONCLUSIONS
This work updates the ETM Collaboration results ofRef. [26] by making six major improvements: i) an anal-ysis of an ensemble of N f = 2 + 1 + 1 of twisted massfermions adding dynamical strange and charm quarks ascompared to an N f = 2 ensemble; ii) a more accurateevaluation of the disconnected contributions yielding themost accurate lattice QCD determination of these quan-tities directly at the physical point to date; iii) the anal-ysis of larger sink-source time separations at higher ac-curacy eliminating more reliably excited states contribu-tions in dominant connected contributions; iv) the com-putation of the GFF B ( Q ) needed for J q,g ; v) the non-perturbative evaluation of the difference between the sin-glet and the non-singlet renormalization functions for allthe relevant operators; and vi) non-perturbative renor-malization of the gluon momentum fraction and angular momentum J g .The major outcomes of this work are:i) The contribution of quarks to the intrinsic protonspin is found to be: (cid:80) q = u,d,s,c ∆Σ q + = 0 . . ≤ ∆Σ ≤ .
18 [84].ii) The verification of the momentum sum for theproton computing all the contributions: (cid:104) x (cid:105) u + + (cid:104) x (cid:105) d + + (cid:104) x (cid:105) s + + (cid:104) x (cid:105) c + + (cid:104) x (cid:105) g = 0 . . . . . . J u + + J d + + J s + + J c + + J g = 0 . . . . . . (cid:80) q = u,d,s,c L q + = 0 . N f =2+1+1 physical ensembles with finer latticespacings and bigger volumes to perform the continuumand infinite volume extrapolations. Acknowledgments [1] M. J. Alguard et al., Phys. Rev. Lett. , 1261 (1976),[,289(1976)].[2] M. J. Alguard et al., Phys. Rev. Lett. , 70 (1978),[,294(1978)].[3] G. Baum et al., Phys. Rev. Lett. , 2000 (1980),[,298(1980)].[4] G. Baum et al., Phys. Rev. Lett. , 1135 (1983),[,302(1983)].[5] J. Ashman et al. (European Muon), Phys. Lett. B206 ,364 (1988), [,340(1987)].[6] J. Ashman et al. (European Muon), Nucl. Phys.
B328 ,1 (1989), [,351(1989)].[7] B. Adeva et al. (Spin Muon (SMC)), Nucl. Instrum.Meth.
A343 , 363 (1994).[8] J. Ball et al., Nucl. Instrum. Meth.
A498 , 101 (2003).[9] C. A. Aidala, S. D. Bass, D. Hasch, and G. K. Mallot,Rev. Mod. Phys. , 655 (2013), 1209.2803.[10] D. de Florian, R. Sassot, M. Stratmann, and W. Vogel-sang, Phys. Rev. Lett. , 072001 (2008), 0804.0422.[11] D. de Florian, R. Sassot, M. Stratmann, and W. Vogel-sang, Phys. Rev. D80 , 034030 (2009), 0904.3821.[12] J. Blumlein and H. Bottcher, Nucl. Phys.
B841 , 205(2010), 1005.3113.[13] E. Leader, A. V. Sidorov, and D. B. Stamenov, Phys.Rev.
D82 , 114018 (2010), 1010.0574.[14] R. D. Ball, S. Forte, A. Guffanti, E. R. Nocera, G. Ridolfi,and J. Rojo (NNPDF), Nucl. Phys.
B874 , 36 (2013),1303.7236.[15] A. Deur, S. J. Brodsky, and G. F. De Tramond, Rept.Prog. Phys. (2019), 1807.05250.[16] H.-W. Lin et al., Prog. Part. Nucl. Phys. , 107 (2018), 1711.07916.[17] E. R. Nocera, R. D. Ball, S. Forte, G. Ridolfi, and J. Rojo(NNPDF), Nucl. Phys. B887 , 276 (2014), 1406.5539.[18] X. Liu and B.-Q. Ma, Eur. Phys. J.
C79 , 409 (2019),1905.02360.[19] E. C. Aschenauer et al. (2013), 1304.0079.[20] P. Djawotho (STAR), Nuovo Cim.
C036 , 35 (2013),1303.0543.[21] A. Adare et al. (PHENIX), Phys. Rev.
D90 , 012007(2014), 1402.6296.[22] A. Accardi et al., Eur. Phys. J.
A52 , 268 (2016),1212.1701.[23] M. A. Clark, B. Jo, A. Strelchenko, M. Cheng, A. Gamb-hir, and R. Brower (2016), 1612.07873.[24] C. Alexandrou, S. Bacchio, and J. Finkenrath, Comput.Phys. Commun. , 51 (2019), 1805.09584.[25] C. Alexandrou, S. Bacchio, M. Constantinou, K. Had-jiyiannakou, K. Jansen, G. Koutsou, and A. VaqueroAviles-Casco (2019), 1909.00485.[26] C. Alexandrou, M. Constantinou, K. Hadjiyiannakou,K. Jansen, C. Kallidonis, G. Koutsou, A. Vaquero Avils-Casco, and C. Wiese, Phys. Rev. Lett. , 142002(2017), 1706.02973.[27] C. Alexandrou, M. Constantinou, K. Hadjiyiannakou,K. Jansen, C. Kallidonis, G. Koutsou, and A. Va-quero Aviles-Casco, Phys. Rev.
D96 , 054507 (2017),1705.03399.[28] H.-W. Lin, R. Gupta, B. Yoon, Y.-C. Jang, and T. Bhat-tacharya, Phys. Rev.
D98 , 094512 (2018), 1806.10604.[29] J. Liang, Y.-B. Yang, T. Draper, M. Gong, and K.-F.Liu, Phys. Rev.
D98 , 074505 (2018), 1806.08366. [30] M. Gockeler, R. Horsley, E.-M. Ilgenfritz, H. Oelrich,H. Perlt, P. E. L. Rakow, G. Schierholz, A. Schiller, andP. Stephenson, Nucl. Phys. Proc. Suppl. , 324 (1997),hep-lat/9608017.[31] R. Horsley, R. Millo, Y. Nakamura, H. Perlt, D. Pleiter,P. E. L. Rakow, G. Schierholz, A. Schiller, F. Winter, andJ. M. Zanotti (QCDSF, UKQCD), Phys. Lett. B714 , 312(2012), 1205.6410.[32] P. E. Shanahan and W. Detmold, Phys. Rev.
D99 ,014511 (2019), 1810.04626.[33] C. Alexandrou, M. Constantinou, K. Hadjiyiannakou,K. Jansen, H. Panagopoulos, and C. Wiese, Phys. Rev.
D96 , 054503 (2017), 1611.06901.[34] Y.-B. Yang, M. Gong, J. Liang, H.-W. Lin, K.-F. Liu,D. Pefkou, and P. Shanahan, Phys. Rev.
D98 , 074506(2018), 1805.00531.[35] M. Deka et al., Phys. Rev.
D91 , 014505 (2015),1312.4816.[36] Y.-B. Yang, R. S. Sufian, A. Alexandru, T. Draper, M. J.Glatzmaier, K.-F. Liu, and Y. Zhao, Phys. Rev. Lett. , 102001 (2017), 1609.05937.[37] C. Alexandrou, M. Constantinou, V. Drach, K. Had-jiyiannakou, K. Jansen, G. Koutsou, A. Strelchenko, andA. Vaquero, Comput. Phys. Commun. , 1370 (2014),1309.2256.[38] A. Stathopoulos, J. Laeuchli, and K. Orginos (2013),1302.4018.[39] C. Michael and C. Urbach (ETM), PoS
LATTICE2007 ,122 (2007), 0709.4564.[40] A. Abdel-Rehim et al. (ETM), Phys. Rev.
D95 , 094515(2017), 1507.05068.[41] R. Frezzotti, P. A. Grassi, S. Sint, and P. Weisz (Alpha),JHEP , 058 (2001), hep-lat/0101001.[42] R. Frezzotti and G. C. Rossi, JHEP , 007 (2004), hep-lat/0306014.[43] X.-D. Ji, Phys. Rev. Lett. , 610 (1997), hep-ph/9603249.[44] X.-D. Ji, Phys. Rev. D52 , 271 (1995), hep-ph/9502213.[45] X.-D. Ji, J. Phys.
G24 , 1181 (1998), hep-ph/9807358.[46] R. Jaffe and A. Manohar, Nuclear PhysicsB , 509 (1990), ISSN 0550-3213, URL .[47] X.-S. Chen, X.-F. Lu, W.-M. Sun, F. Wang, and T. Gold-man, Phys. Rev. Lett. , 232002 (2008), 0806.3166.[48] C. Alexandrou et al., Phys. Rev.
D98 , 054518 (2018),1807.00495.[49] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finken-rath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, andA. Vaquero Aviles-Casco, Phys. Rev.
D100 , 014509(2019), 1812.10311.[50] C. Alexandrou and C. Kallidonis, Phys. Rev.
D96 ,034511 (2017), 1704.02647.[51] C. Alexandrou et al. (2019), 1908.10706.[52] Y. Iwasaki, Nucl. Phys.
B258 , 141 (1985).[53] B. Sheikholeslami and R. Wohlert, Nucl. Phys.
B259 ,572 (1985).[54] C. Alexandrou, S. Gusken, F. Jegerlehner, K. Schilling,and R. Sommer, Nucl. Phys.
B414 , 815 (1994), hep-lat/9211042.[55] S. Gusken, Nucl. Phys. Proc. Suppl. , 361 (1990).[56] M. Albanese et al. (APE), Phys. Lett. B192 , 163 (1987).[57] C. Alexandrou, M. Constantinou, S. Dinter, V. Drach,K. Jansen, C. Kallidonis, and G. Koutsou, Phys. Rev.
D88 , 014509 (2013), 1303.5979.[58] C. Alexandrou, M. Brinet, J. Carbonell, M. Constanti-nou, P. A. Harraud, P. Guichon, K. Jansen, T. Ko-rzec, and M. Papinutto, Phys. Rev.
D83 , 094502 (2011),1102.2208.[59] C. Alexandrou, G. Koutsou, J. W. Negele, andA. Tsapalis, Phys. Rev.
D74 , 034508 (2006), hep-lat/0605017.[60] P. Hagler, J. W. Negele, D. B. Renner, W. Schroers,T. Lippert, and K. Schilling (LHPC, SESAM), Phys.Rev.
D68 , 034505 (2003), hep-lat/0304018.[61] L. Maiani, G. Martinelli, M. L. Paciello, and B. Taglienti,Nucl. Phys.
B293 , 420 (1987).[62] S. Capitani, M. Della Morte, G. von Hippel, B. Jager,A. Juttner, B. Knippschild, H. B. Meyer, and H. Wittig,Phys. Rev.
D86 , 074502 (2012), 1205.0180.[63] W. Wilcox, in
Numerical challenges in lattice quantumchromodynamics. Proceedings, Joint InterdisciplinaryWorkshop, Wuppertal, Germany, August 22-24, 1999 (1999), pp. 127–141, hep-lat/9911013.[64] J. Green, S. Meinel, M. Engelhardt, S. Krieg, J. Laeuchli,J. Negele, K. Orginos, A. Pochinsky, and S. Syritsyn,Phys. Rev.
D92 , 031501 (2015), 1505.01803.[65] A. S. Gambhir, A. Stathopoulos, and K. Orginos, SIAMJ. Sci. Comput. , A532 (2017), 1603.05988.[66] C. McNeile and C. Michael (UKQCD), Phys. Rev. D73 ,074506 (2006), hep-lat/0603007.[67] C. Morningstar and M. J. Peardon, Phys. Rev.
D69 ,054501 (2004), hep-lat/0311018.[68] C. Alexandrou, S. Bacchio, J. Finkenrath, A. Frommer,K. Kahl, and M. Rottmann, Phys. Rev.
D94 , 114509(2016), 1610.02370.[69] S. Bacchio, C. Alexandrou, and J. Finkerath, EPJ WebConf. , 02002 (2018), 1710.06198.[70] K. Goeke, J. Grabis, J. Ossmann, M. V. Polyakov,P. Schweitzer, A. Silva, and D. Urbano, Phys. Rev.
D75 ,094021 (2007), hep-ph/0702030.[71] G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa,and A. Vladikas, Nucl. Phys.
B445 , 81 (1995), hep-lat/9411010.[72] M. G¨ockeler, R. Horsley, H. Oelrich, H. Perlt, D. Petters,P. E. L. Rakow, A. Sch¨afer, G. Schierholz, and A. Schiller,Nucl. Phys.
B544 , 699 (1999), hep-lat/9807044.[73] C. Alexandrou, M. Constantinou, T. Korzec,H. Panagopoulos, and F. Stylianou (ETM Collabora-tion), Phys. Rev.
D83 , 014503 (2011), arXiv:1006.1920.[74] C. Alexandrou, M. Constantinou, T. Korzec,H. Panagopoulos, and F. Stylianou, Phys.Rev.
D86 ,014505 (2012), [arXiv:1201.5025].[75] C. Alexandrou, M. Constantinou, and H. Panagopoulos(ETM), Phys. Rev.
D95 , 034505 (2017), 1509.00213.[76] M. Constantinou et al. (ETM), JHEP , 068 (2010),1004.1115.[77] M. Constantinou, V. Lubicz, H. Panagopoulos, andF. Stylianou, JHEP , 064 (2009), 0907.0381.[78] M. Constantinou, R. Horsley, H. Panagopoulos, H. Perlt,P. E. L. Rakow, G. Schierholz, A. Schiller, and J. M.Zanotti, Phys. Rev. D91 , 014502 (2015), 1408.6047.[79] M. Constantinou, M. Hadjiantonis, H. Panagopoulos,and G. Spanoudes, Phys. Rev.
D94 , 114513 (2016),1610.06744.[80] S. Caracciolo, P. Menotti, and A. Pelissetto, Nucl. Phys.
B375 , 195 (1992).[81] Y.-B. Yang, J. Liang, Y.-J. Bi, Y. Chen, T. Draper, K.- F. Liu, and Z. Liu, Phys. Rev. Lett. , 212001 (2018),1808.08677.[82] S. D. Joglekar and B. W. Lee, Annals Phys. , 160(1976).[83] M. Constantinou, M. Costa, and H. Panagopoulos, Phys.Rev. D88 , 034504 (2013), 1305.1870.[84] C. Adolph et al. (COMPASS), Phys. Lett.
B753 , 18(2016), 1503.08935.[85] M. Engelhardt, Phys. Rev.
D95 , 094505 (2017),1701.01536.[86] C. Alexandrou, J. Carbonell, M. Constantinou, P. A.Harraud, P. Guichon, K. Jansen, C. Kallidonis, T. Ko-rzec, and M. Papinutto, Phys. Rev.
D83 , 114513 (2011),1104.1600.[87] S. Aoki et al. (Flavour Lattice Averaging Group) (2019),1902.08191.[88] N. Sato, W. Melnitchouk, S. E. Kuhn, J. J. Ethier, andA. Accardi (Jefferson Lab Angular Momentum), Phys.Rev.
D93 , 074005 (2016), 1601.07782.[89] J. J. Ethier, N. Sato, and W. Melnitchouk, Phys. Rev. Lett. , 132001 (2017), 1705.05889.[90] R. D. Ball et al. (NNPDF), Eur. Phys. J.
C77 , 663(2017), 1706.00428.[91] S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston,P. Nadolsky, J. Pumplin, C. Schmidt, D. Stump, andC. P. Yuan, Phys. Rev.
D93 , 033006 (2016), 1506.07443.[92] L. A. Harland-Lang, A. D. Martin, P. Motylinski, andR. S. Thorne, Eur. Phys. J.
C75 , 204 (2015), 1412.3989.[93] S. Alekhin, J. Blmlein, S. Moch, and R. Placakyte, Phys.Rev.
D96 , 014011 (2017), 1701.05838.[94] A. Accardi, L. T. Brady, W. Melnitchouk, J. F. Owens,and N. Sato, Phys. Rev.
D93 , 114017 (2016), 1602.03154.[95] H. Abramowicz et al. (H1, ZEUS), Eur. Phys. J.
C75 ,580 (2015), 1506.06042.[96] In the case of the quark singlet EMT operator, the con-version to the MS scheme is actually more involved dueto its mixing with gluon EMT and other gluon opera-tors [81]. However, to one loop it is simplified to Eq. (52).
Appendix A: Expressions for GFFs
The following expressions are provided in Euclidean space. We suppress the Q = − q argument of the generalizedform factors, E N is the nucleon energy for three-momentum (cid:126)q , for the case (cid:126)p (cid:48) = (cid:126)
0, the kinematic factor K = (cid:112) m N / [ E N ( E N + m N )] and Latin indices ( k , n , and j ) take values 1, 2, and 3 with k (cid:54) = j while ρ takes values 1, 2,3, and 4. Π (Γ , (cid:126)q ) = A K (cid:18) − E N − E N m N − m N (cid:19) + B K (cid:18) − E N E N m N + E N m N − m N (cid:19) + C K (cid:18) E N − E N m N + E N m N − m N (cid:19) , (A1)Π (Γ n , (cid:126)q ) =0 , (A2)Π kk (Γ , (cid:126)q ) = A K (cid:18) E N m N q k m N (cid:19) + B K (cid:18) − E N m N + m N − q k E N m N + q k m N (cid:19) + C K (cid:18) − E N m N + m N q k E N m N + q k m N (cid:19) , (A3)Π kk (Γ n , (cid:126)q ) = A K (cid:18) − i (cid:15) k n ρ q k q ρ m N (cid:19) + B K (cid:18) − i (cid:15) k n ρ q k q ρ m N (cid:19) , (A4)Π k (Γ , (cid:126)q ) = A K (cid:18) − i q k − i q k E N m N (cid:19) + B K (cid:18) − i q k i q k E N m N (cid:19) + C K (cid:18) i q k − i q k E N m N (cid:19) , (A5)Π k (Γ n , (cid:126)q ) = A K (cid:18) − (cid:15) k n ρ (cid:18) q ρ q ρ E N m N (cid:19)(cid:19) + B K (cid:18) − (cid:15) k n ρ (cid:18) q ρ q ρ E N m N (cid:19)(cid:19) (A6)Π kj (Γ , (cid:126)q ) = A K q k q j m N + B K (cid:18) − q k q j E N m N + q k q j m N (cid:19) + C K (cid:18) q k q j E N m N + q k q j m N (cid:19) , (A7)Π kj (Γ n , (cid:126)q ) = A K (cid:18) − i (cid:15) k n ρ q j q ρ m N − i (cid:15) j n ρ q k q ρ m N (cid:19) + B K (cid:18) − i (cid:15) k n ρ q j q ρ m N − i (cid:15) j n ρ q k q ρ m N (cid:19) ..