Moments of Nucleon Unpolarized, Polarized, and Transversity Parton Distribution Functions from Lattice QCD at the Physical Point

The second Mellin moments $\langle x\rangle$ of the nucleon's unpolarized, polarized, and transversity parton distribution functions (PDFs) are computed. Two lattice QCD ensembles at the physical pion mass are used: these were generated using a tree-level Symanzik-improved gauge action and 2+1 flavour tree-level improved Wilson Clover fermions coupling via 2-level HEX-smearing. The moments are extracted from forward matrix elements of local leading twist operators. We determine renomalization factors in RI-(S)MOM and match to $\overline{\mathrm{MS}}$ at scale $2\,\mathrm{GeV}$. Our findings show that operators that exhibit vanishing kinematics at zero momentum can have significantly reduced excited-state contamination. The resulting polarized moment is used to quantify the longitudinal contribution to the quark spin-orbit correlation. All our results agree within two sigma with previous lattice results.


I. INTRODUCTION
The distribution of the momentum and spin within a hadron is encoded by parton distribution functions (PDFs).Determining the PDFs is thus an indispensable ingredient to our understanding of the structure of hadrons [1][2][3].There have been various efforts of extracting the PDFs from global fits, for a recent summary see [4].The Lattice QCD community has also achieved remarkable strides in the computation of PDFs over the recent years [5][6][7].
In this study 1 , our focus centers on the evaluation of the second Mellin moment, denoted as ⟨x⟩ [5,[9][10][11], of unpolarized, polarized, and transversity PDFs.We achieve this through the examination of matrix elements of local twist-two operators [10,[12][13][14][15][16][17].This method does not require high momenta to suppress higher-twist contributions, as is needed in calculations that use nonlocal operators, for example the widely used quasi-PDF method [3,6].One of our objectives is to gain insights into the contamination stemming from excited states for different matrix elements and constraining the resulting uncertainty.To attain this objective, a comprehensive investigation of matrix elements at finite but modest momenta becomes imperative, as certain operators have nonvanishing matrix elements exclusively at nonzero momentum.Although the exploration of forward matrix elements of local operators at non-zero momentum is some- 1 Preliminary results were reported in [8].
This paper is organized as follows.In section II we explain our analysis chain and discuss in detail which operators are considered.The different steps of the analysis to extract the matrix elements are shown in section III.We continue the computation of moments in section IV where they become renormalized and averaged over the different results.Further, our findings are put in relation to other Lattice QCD results and global fits.In section V we utilize the moment of the polarized PDF to compute the quark spin-orbit correlation.Last, in section VI we summarize our findings.There are three appendices: appendix A shows extraction of the matrix element for each operator, which are summarized in appendix B. Finally, the calculation of renormalization factors is discussed in appendix C.

II. METHOD
Moments of PDFs can be obtained by calculating forward matrix elements of local leading-twist operators [13,[21][22][23]] Here, the symbol X denotes either V , A, or T , corresponding to the vector, axial, or tensor channels, respectively, and in the tensor case Γ α = σ βγ so that α is a compound index.These channels are associated with unpolarized, polarized, or transversity PDFs.Symmet-rizing the indices and subtracting traces is indicated by braces, {α, µ}.We specifically focus on the isovector channel, which involves the difference between O X for up and down quarks, O X (q = u) − O X (q = d), to avoid calculating disconnected diagrams.The left-right acting covariant derivative ) is constructed on the Euclidean lattice using central finite differences between neighboring lattice points, connected by appropriate gauge links.
It is well understood that these forward matrix elements are proportional to the desired moment ⟨x⟩ [10,12,13].The matrix element is given by In this equation, p represents the 4-momentum of the nucleon.
In the continuum, the operators described in Equation (1) form irreducible representations of the Lorentz group.However, in the context of Euclidean space, the Lorentz group is replaced by the orthogonal group [23].When we transition to the lattice, the orthogonal group further reduces to the hypercubic group H(4).This reduced symmetry can lead to certain operators mixing with lower-dimensional ones.Fortunately, for the specific one-derivative operator studied, such mixing does not occur [24].
Nevertheless, it is important to note that the Euclidean irreducible representations to which our operators belong are divided into multiple hypercubic irreducible representations.In our work, we adopt the common notation, where τ (b) a represents the a th irreducible representation of dimension b.Each of these hypercubic irreducible representations necessitates a distinct renormalization factor.To keep renormalization diagonal, we construct operators with well-defined hypercubic irreducible representations, as suggested by Göckeler et al. [23].
In practical terms, this implies that for each τ (b) a , we must compute the corresponding renormalization factor Z τ (b) a .This factor is subsequently applied to the matrix elements of an operator that transforms irreducibly under the given representation.As a result, we denote the renormalization factor for the operator O The matrix element described in Equation ( 2) can be determined on the lattice by considering the ratios of three-point and two-point correlation functions, as previously discussed in the literature, e.g.[12,13].The twopoint correlation function, denoted as quantifies the correlation between a nucleon source and a nucleon sink separated by a time interval τ .Here we use2 Γ pol = P + [1 − iγ 1 γ 2 ] with P + = (1+γ4) /2 and a nucleon interpolating operator of the form with smeared quark fields q.
The three-point correlation function, denoted as separates the source and sink nucleons by a time interval T while incorporating the operator of interest, O X , at time τ .From here on we let ⃗ p ′ = ⃗ p as indicated in (2).A visual representation is given by Figure 1.The matrix element is extracted in the limit where ≡ lim Once the matrix element is obtained, we can compute the moment by simply dividing the kinematic factor, ⟨x⟩ K = M, with where the a α,µ ∈ R are appropriate factors to express the symmetrization and removal of traces discussed above; This analysis involves a spectral decomposition of the ratio, which allows us to isolate the matrix element of the ground state: To account for the influence of the first excited state, we expand the expression, obtaining the leading contribution from excited states where ∆E represents the energy difference between the first excited state and the ground state (∆E = E 1 − E 0 ).In principle, one would aim to consider large values of T and τ to approach the limit defined in Equation (7).However, it is important to note that as T increases, so does the statistical noise.
The constants in the numerator, R 1 , R 2 , are dependent on the specific operator O X , and their values influence the extent of excited-state contamination in the matrix element.Smaller values of these constants or the presence of certain symmetries can lead to reduced excitedstate contamination in the final result.
In the sum of ratios excited-state contamination is exponentially suppressed with T compared to T /2 for the ratios themselves [25,26].Increasing τ skip reduces excited-state contamination.
Following the proportionality relation of the ratios and desired matrix element (2), we can extract the latter by use of a finite difference.Neglecting excited states, one finds Due to the available data we use a combination of δ /a ∈ {1, 2, 3} depending on whether a neighbour T +δ is available.
The analysis is outlined as follows: Estimation of Ratios: First, we calculate the ratios R(T, τ ) and ratio sums S(T, τ skip ) for each operator.In the unpolarized (V) case we use further, in the polarized (A) case we use and finally for the transversity (T) case we use These have been carefully chosen to have nonzero kinematic factors, compare Equation ( 2), and to be linearly independent.
Matrix Element Extraction: In the next step we extract matrix elements M using two different methods.Method 1 : We extract the slope via finite differences at a specific source-sink separation T = T ′ , compare (12).Method 2 : We obtain the matrix element from a simultaneous (over all source-sink separations) and fully correlated fit to the 2-state form, equation (10).A matrix element obtained through either method is denoted as M| T ′ ,m , where m represents the extraction method.For the second method the T ′ index can be ignored.
From fitting the C 2pt (τ ) we can obtain the groundstate energy E 0 which is used to calculate the kinematic factor.After this, we calculate the unrenormalized moment as To simplify the following equations, we define a compound index j = O X , p, m that runs over all operators and momenta with nonzero kinematic factors as well as the different methods to obtain the matrix element.
Renormalization Factors: We determine the renormalization factors in RI-(S)MOM and match them to MS(2 GeV); for details see appendix C.This allows us to express the renormalized moment as Moment of PDF: To obtain the second moments of PDFs, we define the central value as the weighted average of the different results: Here T j plat denotes the smallest source-sink separation such that X j (T ′ ) agree for all T ′ ≥ T j plat .Naturally, the sum over T ′ does not apply for the second method, where we fit the 2-state function, as there is no T ′ to consider.The weights W j (T ′ ) ∝ 1 /σ 2 j (T ′ ) are normalised in such a way that weights associated to sum ratios sum to 1 /2 as do the weights for the 2-state fit.The used variances are estimated via bootstrap over X j (T ′ ) and the errors of the renormalization constants are propagated.
Systematic Error Estimation: Finally, we estimate a systematic error, constraining the uncertainty from excited state contamination, by calculating the weighted standard deviation over the different results: Again the sum over T ′ is not applied for the 2-state fit.
Relation to Quark Spin-Orbit Correlations: The longitudinal quark spin-orbit correlation L q ℓ S q ℓ in the proton (where the subscript ℓ denotes alignment with the direction of motion of the proton) is related to the generalized transverse momentum-dependent parton distribution (GTMD) G q 11 [27] as in Equation ( 16), which in turn can be related to the generalized parton distributions (GPDs) H q , H q , E q T and H q T [28,29] as in Equation (17), where all distribution functions are quoted according to the nomenclature of [30] and are taken in the forward limit; k T denotes the quark transverse momentum.H q is the standard chiral-even helicity GPD and H q is the standard chiral-even unpolarized GPD; E q T and H q T are chiral-odd GPDs.The longitudinal quark spin-orbit correlation has been evaluated according to Equation (16) in Ref. [31]; on the other hand, the results of the present work can be used complementarily to access the correlation via Equation (17), which can be viewed as the axial analogue of Ji's sum rule for orbital angular momentum: At the physical pion mass, the term proportional to m q /m N is negligible.In the forward limit, dx H q corresponds to the number of valence quarks, i.e., unity in the isovector, u − d quark case considered here.Therefore, to an excellent approximation, one has where dx x H q = ⟨x⟩ A in the forward limit has been identified.The results obtained in the following section will be used to quantify the longitudinal quark spin-orbit correlation and will also be confronted with the results of Ref. [31].
Simulation Parameters: We use a tree-level Symanzik-improved gauge action with 2+1 flavour treelevel improved Wilson Clover fermions coupling via 2level HEX-smearing.Detailed information about the simulation setup can be found in references [32][33][34].Key simulation parameters are summarized in Table I.Two ensembles, coarse and fine, at the physical pion mass are used.These ensembles correspond to lattice spacings of 0.1163(4) fm and 0.0926(6) fm, respectively.As described in [34], the smearing is done using Wuppertal smearing [35] q ∝ (1 + αH) N q with H being the nearestneighbor gauge-covariant hopping matrix -at α = 3 and N = 60, 100 for the coarse and fine ensemble, respectively.For each ensemble, two-point and three-point correlation functions are calculated.These calculations involve source-sink separations ranging from approximately 0.3 fm to 1.4 fm for the coarse ensemble and approximately 0.9 fm to 1.5 fm for the fine ensemble.Furthermore, we consider two different momenta: ⃗ p = (p x , 0, 0) with p x = 0, −2[ 2π /L] for the coarse ensemble, and with p x = 0, −1[ 2π /L] for the fine ensemble.

III. ESTIMATION OF MATRIX ELEMENTS
In Figure 2, we present results obtained from the coarse ensemble, using two different operators O X per channel, as shown in the upper and lower rows.Each column is dedicated to a particular channel: From top to bottom we display the operators 2. and 3. (unpolarized), 1. and 2. (polarized), and 2. and 3. (transversity).Different source-sink separations are represented by various colours, while momenta are distinguished using filled circles for zero momentum and unfilled squares for finite momentum; this is kept consistent throughout all figures.A plateau in these plots corresponds to the matrix element of the shown operator.To simplify comparison we directly translate this to the bare moment, by multiplying with the kinematic factor R(T, τ ) = 1 /K • R(T, τ ).It is worth noting that we exclude the largest source-sink separation from these plots due to its substantial statistical uncertainty.
These operator choices are intentionally selected to illustrate the extreme variability of the excited-state contamination.While the upper row has a clearly visible cosh behavior -as expected from the 2-state function (10) -the lower row remains perfectly flat within statistics.Moreover, we observe that the convergence in sourcesink separation is much faster for the lower row.For instance, in the lower row the plateau already converges after T /a = 3 while the upper row requires T /a ≥ 8 in these particular examples.This rapid convergence in the lower row is noteworthy, but it also comes with a drawback, which we observe in general, compare analysis summary plots in appendix A for the other operators: Operators that exhibit such flat behavior at small source-sink separations have a vanishing kinematic factor at zero momentum, making them computationally more challenging to handle.
In Figure 3, we present sum ratios, using the same operators as in Figure 2 but put into one subplot.The upper and lower row represent the coarse and fine ensemble, respectively.The value of τ skip = 1 is fixed as the slope of the summed ratios did not change for larger values.The presence of excited-state contamination is subtly hinted at by the slight curvature observed in the data, although it is considerably less pronounced compared to the ratios.
In Figure 4, we present the result of the matrix element extraction for the same operators as displayed in Figure 2. Similar plots for all used operators can be found in appendix A. We plot horizontal lines to represent the average (over T ′ ≥ T j plat ) slope of the summed ratios, .A larger and a smaller lattice spacing, labelled as "Coarse" and "Fine" respectively, are available.The ensembles were generated with a tree-level Symanzikimproved gauge action with 2+1 flavour tree-level improved Wilson Clover fermions coupled via 2-level HEX-smearing [32][33][34].Furthermore, the available source-sink separations (T ) and momenta (px) which are used in the calculation of the ratios, Equation (7) . Ratios, cf.Equation (7), for the coarse ensemble.Various source-sink separations T are represented by different colors, while the two momenta are distinguished using filled circles and unfilled squares.Each subplot corresponds to a different operator from the different channels organized by column.For the Unpolarized (V) case we display operators 2. and 3.; for the polarized (A) case we display operators 1. and 2.; and for the transversity (T) case we display operators 3. and 1. for the upper and lower panel, respectively.
divided by the kinematic factor.These slopes are extracted with the finite difference approach (12).As the matrix element is given by a plateau of the ratios, the expectation is that the plotted slope agrees at least with the central points τ ∼ 0 of large source-sink separations, which can be verified for all operators within uncertainty.Again, those operators which are already flat match this expectation for more points and at smaller source sink separation.
Following the axolotl-like shape of the ratios, the solid, i.e. zero momentum, and dashed, i.e. finite momentum, lines indicate the central value 2-state fit result, using the form (10).The area around these indicate statistical un-certainty obtained via fitting on each bootstrap sample.We use all data points that are covered by the best fit plot in a (T, τ )-simultaneous fit.This presents a fit interval in τ /a which has been chosen by minimizing a χ 2 /dof.The smaller source-sink separations for the coarse ensemble are excluded by this condition, as no points were left in the fit interval.
Considering all fits, values of χ 2 /dof range from 0.4 to 2.7.Correlations which go into these were estimated over the bootstrap samples of the included points and then kept fixed for the central value fit as well as the fits per bootstrap sample.
Notably, the values of the matrix element obtained Unpolarized (V)  from summed ratios and 2-state fits always agree within statistics.The latter has reduced statistical uncertainty.

IV. MOMENTS OF PDFS
In Figure 5, we illustrate the results for the renormalized moments, which are extracted from the summed ratios (shown in grey, defined in Equation ( 11)) and the 2-state fits (displayed in red, as defined in Equation (10)) 3 .The final average is denoted by the blue solid line, while its statistical uncertainty is indicated by the blue dot-dashed lines.The blue band represents the combined statistical and systematic uncertainty, as outlined in Equation (15), added in quadrature.The light gray points are not included in the average, in accordance with the T j plat constraint.To enhance the resolution of the relevant data points, the ordinate limit is truncated at 4σ and centered around the final average.The numerical values of the final averages can be found in Table II.
Comparing the two ensembles we find agreement within statistics indicating a flat continuum extrapolation.With only two points a reliable extrapolation is not possible.The best we can do is to interpret the data points as Gaussian distributions, with mean equaling the central value and width given by the uncertainties added in quadrature, and perform a Bayesian fit.The relevant scale of discretization effects [24,36] is aΛ QCD resulting in a term proportional to α s aΛ QCD .The operators themselves have tree-level quadratic discretization effects, resulting in the extrapolation We use Gaussian priors for the coefficients, p mi = N (0, 2) and no prior on the continuum value ⟨x⟩ ren cont .We approximate α s ≈ 0.3 which is sufficient due the fact that the coefficients m i are mainly constrained by the prior.The continuum-extrapolated results are likewise given in Table II.
Table II.Final averages for the second moments of PDFs in the unpolarized, polarized and transversity channels, compare Figure 5.For the coarse and fine ensemble results, the central value is obtain as a weighted average over the different operators, momenta, and extraction methods, cf.Equation (14).Further, the statistical uncertainty (first uncertainty) comes from a bootstrap over the original ensemble, while the systematic uncertainty (second uncertainty) is computed using the weighted standard deviation over the same set of results, cf.Equation (15).We extrapolated the two points to the continuum limit using a Bayesian fit approach assuming them to be independent and Our results are in good agreement, at the level of one to two standard deviations, with moments previously computed by other Lattice QCD collaborations [14-17, 37, 38].Moreover, confronting our moments with phenomenological extractions, the comparison is quite favorable in the case of the axial moment, with Ref. [39] giving ⟨x⟩ A = 0.190 ± 0.008.On the other hand, in the unpolarized case, a certain tension between lattice and phenomenological results remains, with the recent determination in Ref. [4], ⟨x⟩ V = 0.143 (5), differing from our result by about three standard deviations.

V. QUARK SPIN-ORBIT CORRELATION
With the results from Table II we can calculate the longitudinal quark spin-orbit correlation in the proton according to Equation (18).The obtained values can be found in Table III, along with the result obtained using the GTMD approach, Equation ( 16), in Ref. [31].The results are in good agreement.As discussed in more detail in Ref. [31], the magnitude of this direct correlation between the spin and the orbital angular momentum of a quark significantly exceeds the correlation induced by the quark being embedded in a polarized proton environment.There is, therefore, a strong direct dynamical coupling between quark orbital angular momentum and spin, reminiscent of the jj coupling scheme in atomic physics, rather than the Russell-Saunders coupling scheme.

VI. SUMMARY
In this study, we compute the second Mellin moment ⟨x⟩ of the unpolarized, polarized, and transversity parton distribution functions using two lattice QCD ensembles at the physical pion mass.Our approach involves extracting forward nucleon matrix elements at both zero and finite momentum, boosted in the x-direction.Through the finite momentum data, we identify operators that exhibit remarkably small excited-state contamination.Given the two ensembles a reliable continuum extrapolation is not accessible.Regardless, we apply a Bayesian fit, accepting a strong dependence on the choice of priors, to provide a continuum estimate.The resulting values are in agreement with both individual ensembles: ⟨x⟩ u + −d + = 0.200 (17), ⟨x⟩ ∆u − −∆d − = 0.213 (16), and ⟨x⟩ δu + −δd + = 0.219 (21).Furthermore, we extract the isovector longitudinal quark spin-orbit correlation in the proton using the moment of the polarized PDF, 2L q ℓ S q ℓ = −0.393(08).We find good agreement with earlier calculations based on GTMDs [31].3 (10).The final average is displayed as a blue solid line while its statistical uncertainty is indicated via the blue dot-dashed line.The blue band represents the statistical and systematic uncertainty, cf.Equation (15), added in quadrature.The light gray points are not included in the average as per the T j plat constraint, compare (14), to reduce excited-state effects.The ordinate limit is cut at 4σ around the final average to increase resolution of the relevant points.0.00 0.02 0.04 0.06 0.08 0.10 0.12 a[fm] 0.17

Transversity (T)
x ren a x ren cont (1 + m 1 s a QCD + m 2 (a QCD ) 2 ) Figure 6.Continuum extrapolation using a Bayesian fit with the model described in Equation (19).The limited amount of data makes this extrapolation strongly dependent on the chosen priors for the coefficients mi.Coming from a power counting these are expected to be of O(1), reflected in Gaussian priors of mean zero and variance 2. Resulting estimates are listed in Table II.

Appendix A: Results Per Operator
In this appendix we show the analysis summary resolved per operator.As before, coarse and fine ensemble results are displayed in the first and second row, respectively.Different colors represent different source sink separations and the horizontal dash-dotted, i.e. zero momentum, and dotted, i.e. finite momentum, lines represent the average (over T ′ ≥ T j plat ) slope of the summed ratios, divided by the kinematic factor.These slopes are extracted with the finite difference approach (12).The solid and dashed curves are the best-fit result of the 2state fit to (10), the surrounding band corresponds to the bootstrap uncertainty of the fit.Figure 7 displays the analysis of the operators corresponding to the unpolarized (vector) PDFs. Figure 8 displays the analysis of the operators corresponding to the polarized (axial) PDFs.Figures 9 and 10 display the analysis of the operators corresponding to the transversity (tensor) PDFs.As mentioned in Section III, agreement of the slope of summed ratios with the plateau region expected around τ = 0 is given for all operators within one sigma.Corresponding best 2-state fit lines are in perfect agreement with the data points.

Appendix B: Summary Plots
We present summary plots of the moments for the coarse 11 and fine 12 ensemble.The three channels, unpolarized (V), polarized (A), and transversity (T) are shown in the columns.Each result, i.e. the different operators and momenta, is displayed in the panels separated by the dotted and dashed lines.The solid black line separates the sum-ratio method, points in purple, and the 2-state fit method, points in red.For the sum-ratio method the different T ′ are spread across the abscissa.As a point of reference, the average over the points, as described in equation ( 14), is shown by the horizontal blue line, with the statistical uncertainty shown by the dotted dashed line and the combined uncertainty by the blue band.This corresponds to the blue line in figure 5.A (strong) dependence on the source-sink separation can be seen in the sum-ratio related points.
Analysis results of the ratios (points), slope of summed ratios (horizontal lines) and 2-state fit results (curves) for the operators 1., 2. and 3. corresponding to the unpolarized (vector) PDF.
Analysis results of the ratios (points), slope of summed ratios (horizontal lines) and 2-state fit results (curves) for the operators 1. and 2. corresponding to the transversity (tensor) PDF.Transversity . Analysis results of the ratios (points), slope of summed ratios (horizontal lines) and 2-state fit results (curves) for the operators 3. and 4. corresponding to the transversity (tensor) PDF.
Summary plot of the renormalized moment of PDF of the coarse ensemble resolved for each operator O X , represented by the corresponding ID and irrep, and momentum px.The purple points correspond to results obtained by the sum-ratio method evaluated at a source-sink separation T ′ .The red points are obtained using the two-state fit.The blue solid line corresponds to the overall average as described by equation ( 14  We determine renormalization factors for isovector vector, axial, and tensor one-derivative twist-two operators using the nonperturbative Rome-Southampton approach [46], in both RI ′ -MOM [46,47] and RI-SMOM schemes [48], and convert and evolve to the MS scheme at scale 2 GeV using perturbation theory.We label these renormalization factors Z ρ DV , Z ρ DA , and Z ρ DT for the onederivative vector, axial, and tensor operators, respectively, with ρ denoting the irreducible representation of the hypercubic group that takes on two possible values in each case.
We largely follow our earlier work that used operators with no derivatives [34].Our primary data are the Landau-gauge quark propagator, and the Landau-gauge Green's functions for operator O, (C2) Here O is an isovector quark bilinear with one derivative, yielding one Wick contraction: a connected diagram.We evaluate these objects using four-dimensional volume plane wave sources, yielding an effectively large sample size from the volume average.From these, we construct our main objects, the amputated Green's functions, Provided that O belongs to a definite irreducible representation of the hypercubic group, these renormalize diagonally: To avoid determining Z ψ directly, we will form ratios to determine Z O /Z V and take Z V computed from pion three-point functions in Ref. [34].

Conditions and matching
The RI ′ -MOM scheme uses kinematics p ′ = p, whereas RI-SMOM uses p 2 = (p ′ ) 2 = q 2 with q = p ′ − p.In both cases, the scale is defined as µ 2 = p 2 .For the vector current, we impose the conditions listed in Ref. [34] on Λ R Vµ to determine Z V /Z ψ .
For the one-derivative operators, we start with the continuum decomposition of the amputated Green's function Λ Oµν... (p ′ , p) into a sum of products of O(4)invariant functions Σ O,ρ (p 2 ).Our choice of decomposition, given below, is such that at tree level, Σ O (p 2 ) = δ i1 , and our renormalization conditions will all be of the form Basing this condition on a O(4)-invariant function computed within each irrep ensures that the ratio of renormalization factors for two different lattice irreps of the same continuum operator is scale and scheme-invariant.
The one-derivative vector operator is where S takes the symmetric traceless part of the tensor: Our decomposition for the RI ′ -MOM scheme is a scaled version of the one used by Gracey [49]: where here and below we neglect tensors of opposite chirality.For RI-SMOM, the derivative in the operators basis used by Gracey [50] did not yield a definite C-symmetry, unlike our operator containing ↔ D. This allows us to use half as many tensors as Gracey; see also [51,52].Defining p = (p ′ + p)/2, our tensors and their relation to Gracey's tensors P W2 (i) are the following: , (C10) , , where p 2 = 3 4 µ 2 and the square brackets denote antisymmetrization.The one-derivative axial operator, is related to the vector operator by chiral symmetry and its tensor structures correspond to those of the vector operator, multiplied by γ 5 .We use the four-loop anomalous dimension [53][54][55] 4 and three-loop matching to MS [50,52].
The one-derivative tensor operator is where the symmetrization and trace subtraction has the form [49] ST µνρ = 1 2 Tµνρ + Tµρν − 1 3 with Tµνρ = 1 2 (T µνρ − T νµρ ).Choosing to start by antisymmetrizing µν leaves us with only two tensor structures in the RI ′ -MOM scheme, compared with Gracey's three: For RI-SMOM, the supplementary data of Ref.
[57] lists 30 structures.First antisymmetrizing µν reduces this to 16 and charge conjugation further reduces the number to 8. We write the first six as The last two tensors involve γ 5 or the identity, and they have vanishing trace with each of the first six, so they can be neglected.We use the three-loop anomalous dimension [49], the three-loop matching from RI ′ -MOM [49], and the two-loop matching from RI-SMOM [57].
Our numerical setup follows Ref. [34], extended to include both sets of momenta on both ensembles.We use partially twisted boundary conditions, namely periodic in time for the valence quarks.The plane wave sources are given momenta either along the four-dimensional diagonal p . By contracting them in different combinations, we get data for both RI ′ -MOM kinematics, p ′ −p = 0, and RI-SMOM kinematics, p ′ − p = 2π L (0, 0, 0, ±2k) or ± 2π L (0, k, −k, 0).We used 54 gauge configurations from each ensemble; however, on one configuration of the coarse ensemble the valence twisted boundary condition yielded a near-singular Dirac operator and the multigrid solver was unable to converge.We omitted this configuration, using only 53 on the coarse ensemble.
Perturbatively matching from RI ′ -MOM or RI-SMOM to the MS scheme and evolving to scale 2 GeV does not eliminate dependence on the initial scale µ: there are remaining effects from lattice artifacts, truncation of the perturbative series, and nonperturbative contributions.To control these artifacts, we perform fits including terms polynomial in µ 2 and also a pole term, following Ref.[58].These fits have the form A + Bµ 2 + Cµ 4 + D/µ 2 , where the constant term A serves as our estimate of the renormalization factor ratio Z O /Z V .Results for this ratio, choosing one irrep for each of the three operator types, are shown in Fig. C 2. For each operator and scheme, we perform three fits: using the 4D data with two different fit ranges µ 2 ∈ [4, 20] and [10, 30] GeV 2 and using the 2D data with µ 2 ∈ [4, 15] GeV 2 .In some cases (particularly using the very precise RI-SMOM data), the fit quality is very poor and thus we scale the statistical uncertainty by χ 2 /dof whenever this is greater than one.Following the same prescription as in Ref. [34], we combine the results first within each scheme and then for our final result using both schemes, estimating the systematic uncertainty (which is dominant) from the scatter of results and conservatively taking the maximum statistical uncertainty.In all cases, there is good agreement between the two schemes.
Renomalization-group-invariant ratios of renormalization factors in different irreps Z in the RI-SMOM scheme using the 4D kinematics.Because in the continuum and infinite volume these ratios are independent of µ 2 , we can fit using much lower momenta and only avoid the region µ 2 < 1 GeV 2 due to finite-volume effects.We also omit the pole term, i.e. set D = 0.In all cases, we obtain results within a few percent of unity.
Our final values for the ratios of renormalization factors are given in       DT /ZV on the coarse (left) and fine (right) ensembles, determined using the RI ′ -MOM (green circles) and RI-SMOM (orange squares) intermediate schemes together with momenta along the four-dimensional diagonal (filled symbols) and two-dimensional diagonal (open symbols) and then matched to MS at scale 2 GeV.For most points, the statistical uncertainty is smaller than the plotted symbol.The solid curves are fits to the 4D data in the µ 2 range from 4 to 20 GeV 2 , the dashed curves in the range from 10 to 30 GeV 2 , and the dotted curves are fits to the 2D data in the range from 4 to 15 GeV 2 .The fit curves without the pole term are also plotted in the range 0 < µ 2 < 6 GeV 2 .To reduce clutter, uncertainties on the fits are not shown.The symbols filled with black near µ 2 = 0 provide the final estimat for each intermediate scheme; their outer (without end cap) and inner (with end cap) error bars show the total and statistical uncertainties.The filled dark gray diamonds are the final estimates that combine both schemes.DT , determined using the RI ′ -MOM (green circles) and RI-SMOM (orange squares) intermediate schemes together with momenta along the four-dimensional diagonal (filled symbols) and two-dimensional diagonal (open symbols) and then matched to MS at scale 2 GeV.For most points, the statistical uncertainty is smaller than the plotted symbol.The solid curves are fits to the 4D data in the µ 2 range from 1 to 8 GeV 2 and the dotted curves are fits to the 2D data in the range from 1 to 5 GeV 2 .To reduce clutter, uncertainties on the fits are not shown.The symbols filled with black near µ 2 = 0 provide the final estimate for each intermediate scheme; their outer (without end cap) and inner (with end cap) error bars show the total and statistical uncertainties.The filled dark gray diamonds are the final estimates that combine both schemes.

Figure 1 .
Figure 1.Graphical representation of C O X 3pt (T, τ ): a source nucleon inserted at time t = 0 and a sink nucleon removed at time t = T .A local leading twist operator (1) is inserted on a given time slice τ .The nucleons on the lattice are represented by interpolating operators χ (4) while O X is determined by finite differences connected with gauge links.

Figure 3 .
Figure 3. Ratio sums S(T, τ skip ) on the coarse and fine lattice, employing the same operators as in 2. Each S(T, τ skip ) is plotted at fixed τ skip = 1.As in Fig.2, different momenta are displayed with hollow and filled markers.

Figure 4 .
Figure 4. Extraction result for the matrix elements, cf.Equation(13), plotted on top of the original ratio data.The same operators and layout as in Figure2are used.Coarse and fine ensemble results are displayed in the first two and second two rows, respectively.Dot-dashed and dotted horizontal lines represent the average slope of the summed ratios divided by the kinematic factor.Solid and dashed lines represent the simultaneous and fully correlated central value fit to the ratios using the 2-state form(10).Surrounding colored areas represent bootstrap uncertainties.
Figure 11.Summary plot of the renormalized moment of PDF of the coarse ensemble resolved for each operator O X , represented by the corresponding ID and irrep, and momentum px.The purple points correspond to results obtained by the sum-ratio method evaluated at a source-sink separation T ′ .The red points are obtained using the two-state fit.The blue solid line corresponds to the overall average as described by equation (14) with the dashed line indicating statistical uncertainty and the blue area indicating statistical and systematic uncertainty added in quadrature.

Figure 12 .
Figure 12.Summary plot of the renormalized moment of PDF of the fine ensemble, similar to figure11

Table I .
Details of the used ensembles.The ensembles are at the physical pion mass, mπ ≈ m phys π , are displayed.
Gaussian distributed with mean equaling the central value and standard deviation coming from the combined statistical and systematical uncertainty, compare Figure IV.

Table III .
(18)ced isovector longitudinal quark spin-orbit correlation, estimated using the results for the polarized (A) moment shown in TableIIand relation(18).