Lattice results for the longitudinal spin structure and color forces on quarks in a nucleon

Using lattice QCD, we calculate the twist-2 contribution $a_2$ to the third Mellin moment of the spin structure functions $g_1$ and $g_2$ in the nucleon. In addition we evaluate the twist-3 contribution $d_2$. Our computations make use of $N_f=2+1$ gauge field ensembles generated by the Coordinated Lattice Simulations (CLS) effort. Neglecting quark-line disconnected contributions we obtain as our best estimates $a_2^{(p)}= 0.069(17)$, $d_2^{(p)}= 0.0105(68)$ and $a_2^{(n)}= 0.0068(88)$, $d_2^{(n)}= -0.0009(70)$ for the proton and the neutron, respectively, where we use the normalizations given in Eqs. (58) and (59). While the $a_2$ results have been converted to the $\overline{\mathrm{MS}}$ scheme using three-loop perturbation theory, the numbers for $d_2$ are given in the regularization independent momentum subtraction (RI$^\prime$-MOM) scheme, i.e., the conversion has been performed only in tree-level perturbation theory. The $d_2$ results can be interpreted as corresponding to a transverse color Lorentz force on a quark in a transversely polarized proton of size $F^{(u)} = 116(61)$ MeV/fm and $F^{(d)} = -38(66)$ MeV/fm for $u$ and $d$ quarks, respectively. The error estimates quoted include statistical and systematic uncertainties added in quadrature.

Using lattice QCD, we calculate the twist-2 contribution a2 to the third Mellin moment of the spin structure functions g1 and g2 in the nucleon. In addition we evaluate the twist-3 contribution d2. Our computations make use of N f = 2 + 1 gauge field ensembles generated by the Coordinated Lattice Simulations (CLS) effort. Neglecting quark-line disconnected contributions we obtain as our best estimates a (p) 2 = 0.069 (17), d = −0.0009 (70) for the proton and the neutron, respectively, where we use the normalizations given in Eqs. (58) and (59). While the a2 results have been converted to the MS scheme using three-loop perturbation theory, the numbers for d2 are given in the regularization independent momentum subtraction (RI -MOM) scheme, i.e., the conversion has been performed only in tree-level perturbation theory. The d2 results can be interpreted as corresponding to a transverse color Lorentz force on a quark in a transversely polarized proton of size F (u) = 116(61) MeV/fm and F (d) = −38 (66) MeV/fm for u and d quarks, respectively. The error estimates quoted include statistical and systematic uncertainties added in quadrature.

I. INTRODUCTION
For a number of reasons hadron spin structure has attracted intense interest for more than two decades with no sign of attenuation. Quite to the contrary, CEBAF@12GeV [1] and the EIC [2] promise to bring such investigations to a higher level both with respect to the precision and the variety of observables investigated. This provides a strong motivation also to update theory predictions for the relevant spin-dependent quantities. The central goal of the JLAB and BNL programs is to better understand the structure of hadrons. This includes multiparton correlations, which are parametrized by higher-twist coefficients.
The most prominent such example is the matrix element d 2 , the third Mellin moment (x 2 ) of the twist-3 contribution to the helicity structure function g 2 (x, Q 2 ) of deep-inelastic longitudinally polarized lepton-nucleon scattering. As this corresponds to the lowest-dimensional nontrivial chiral-even twist-3 matrix element, d 2 is of particular theoretical and phenomenological interest. For instance, the same correlations of quarks and gluons constitute the leading contribution to the Qiu-Sterman distributions [3], which play a central role in the collinear factorization of single spin asymmetries. These distributions represent the limit of vanishing impact parameter b ⊥ = 0 of the Sivers functions [4], i.e., of the transverse-momentum dependent parton distribution functions f ⊥ 1T (x, b ⊥ , Q 2 ) that describe the distribution of an unpolarized quark inside a transversely polarized nucleon. The measurement of the Sivers distributions in polarized semi-inclusive deep-inelastic scattering and in Drell-Yan experiments is one of the main goals of the experimental programs at JLAB and the EIC. For more details, see, e.g., Refs. [5,6].
Neglecting the twist-3 contributions, g 2 can be ob-tained from the helicity structure function g 1 (x, Q 2 ). This involves invoking the well-known Wandzura-Wilczek relation [7,8]. However, a remarkable property of g 2 is that the twist-3 contribution is not power suppressed in 1/Q, relative to its twist-2 part. Nevertheless, its determination from longitudinally polarized deep-inelastic scattering experiments alone still represents a serious challenge and new high precision measurements are planned at JLAB and the EIC. In order to match the expected statistical precision of the planned experiments, a much improved theoretical understanding of higher-twist contributions is needed and d 2 is the ideal starting point. The matrix element d 2 has another very interesting phenomenological interpretation: As was argued in Refs. [9,10] it is related to the average transverse color Lorentz force acting on a quark q in a nucleon which moves in the z direction and is transversely polarized. More explicitly, the y-component of the color Lorentz force is given by Here the four-momentum p of the nucleon state |p, s has been chosen as (p µ ) = p + (1, 0, 0, 1)/ √ 2, the spin vector s is normalized according to s 2 = −m 2 N , G µν denotes the color field strength tensor and g is the strong coupling constant. Using the lattice results for the proton presented in Ref. [11] (coauthored by some of us) estimates for this force were published some time ago in Ref. [9]. It appears, however, that these estimates arXiv:2111.08306v2 [hep-lat] 9 Mar 2022 were affected by a misunderstanding of the respective conventions and by a sign error noted later in Ref. [10]. The corrected numbers differ by a factor − 1 2 from those given in Ref. [9] and read F u,y (0) = 50 ± 60 MeV/fm , Unfortunately, the errors, which are purely statistical, are very large. Systematic uncertainties were not estimated, in particular, those arising from finite lattice spacing. Meanwhile several experiments have extracted estimates of d (p) 2 and d (n) 2 [12][13][14][15][16][17][18][19], where the superscript indicates proton or neutron, respectively. These estimates are found to be quite small compared to various model predictions but compatible with the old lattice results (considering the large error bars), see, e.g., Fig. 2 of Ref. [19]. (Actually, the lattice results in this figure should also have been divided by 2.) We remark that with the natural energy scale for the force F being Λ 2 QCD one would not expect d 2 to be much smaller than the central values of this early lattice calculation given in Eq. (2). So there is hope that with a moderate reduction of the lattice uncertainties, this time also including systematics, one may be able to demonstrate that d 2 and thus the average color force F is different from zero.
Let us stress that the experimental and lattice investigations of d are only meant to be the starting point of much broader investigations. For example, it was also argued in Ref. [10] that there exists an analogous relationship between generalized parton distributions and force distributions F i λ λ ( b ⊥ ) in the transverse plane: Here λ and λ denote the nucleon polarization and ∆ ⊥ is the transverse momentum conjugate to the impact parameter b ⊥ . Another interesting result was derived in the very recent paper [20], where QCD factorization for quasidistributions was analyzed up to twist-3. Approaches based on so-called quasi-and pseudodistribution functions have gained prominence in lattice QCD calculations of hadron structure observables, due to their prospect of providing information that goes beyond the computation of Mellin moments of (generalized) parton distribution functions, distribution amplitudes etc., see Refs. [21][22][23][24][25][26][27] and references therein.
In Ref. [20] it was shown that for quasidistributions the Wandzura-Wilczek relation [7,8] is modified such that twist-2 and twist-3 contributions stay mixed, making their separate determination on the lattice far more difficult. Knowing d 2 from a direct lattice calculation would obviously help to unravel the different contributions.

II. OPE AND RENORMALIZATION IN THE CONTINUUM
A leading-order OPE (operator product expansion) analysis with massless quarks shows that the moments of g 1 and g 2 can be written as [28] where q runs over the light quark flavors with charges Q q and µ denotes the renormalization scale. Equations (6) and (7) hold for even n, with n ≥ 0 for the former and n ≥ 2 for the latter. The Wilson coefficients E (q) 1,n and E (q) 2,n depend on the ratio of scales µ 2 /Q 2 and the running coupling constant g(µ), To the best of our knowledge, the loop corrections for E 2,n have not yet been calculated, while they are known up to two-loop order for E 1,n [29]. The first-order corrections are flavor independent, The reduced matrix elements a (q) n (µ) and d (q) n (µ) are defined as [28] p, s|O p, s|O 5(q) [σ{µ1]µ2···µn} |p, s in terms of matrix elements of the local operators (12) in the nucleon state |p, s . Here D and the symbol {· · · } ([· · · ]) indicates symmetrization (antisymmetrization) of the enclosed indices. The operator in Eq. (10) has twist two, whereas the operator in Eq. (11) has twist three. As far as the Wilson coefficients may be considered as flavor independent, we can define a n and d n for the nucleon as Remarkably, in the moments (7) of g 2 the twist-3 matrix elements d (q) n (µ) are not suppressed relative to the twist-2 matrix elements a (q) n (µ). Note that our definitions of a 2 and d 2 have been taken from Ref. [28]. In many publications alternative definitions are employed, where a alt 2 = a 2 /2 and d alt 2 = d 2 /2.
Utilizing the equations of motion of massless QCD and the relation [D µ , D ν ] = −igG µν , the twist-3 operators O

5(q)
[σ{µ1]µ2···µn} can be rewritten in a manifestly interaction-dependent form. For n = 2 one finds whereG µν = 1 2 µνρσ G ρσ is the dual gluon field strength tensor and the totally antisymmetric tensor is such that 0123 = 1. Therefore we can define the reduced matrix element d 2 in the chiral limit also by (see, e.g., Ref. [30]) The Wilson coefficients (8) can be computed in perturbation theory, while the nucleon matrix elements a The renormalization of the operators which contribute to the moments of g 2 has been studied by several authors in continuum perturbation theory [31][32][33][34][35][36][37][38]. For example, in Refs. [37,38] the following operators are considered for n = 2 in the flavor-nonsinglet sector: The gluon field strength tensor G αβ could alternatively be expressed in terms of a commutator of two covariant derivatives. As Eqs. (11) and (12) show, the matrix element d 2 corresponds to the nucleon matrix elements of the renormalized operators (17). The operators (17)- (20) are linearly dependent: In the massless case, this relation leads to Eq. (15) upon application of the equations of motion.
Calculating the quark-quark-gluon three-point functions with a single insertion of each of these operators in oneloop perturbation theory, one sees that also a gauge-variant operator has to be taken into account in the process of renormalization: Of course, in physical matrix elements neither R eq nor R eq1 will contribute. They show up, however, in offshell vertex functions and influence the renormalization factors.

III. LATTICE OPERATORS AND RENORMALIZATION
In the following, we use Euclidean notation. For our lattice evaluation of the reduced matrix elements a 2 and d 2 we construct discretized versions of the relevant operators. In the process of the renormalization of these operators, operator mixing requires particular attention because the discrete symmetry group H(4) of a hypercubic lattice is less restrictive than the continuous symmetry group O(4) of Euclidean spacetime.
In the case of the twist-2 matrix element a 2 we use the four-dimensional multiplet of operators spanned by where The operators (23) transform according to the representation τ (4) 3 of the hypercubic group H(4) [39][40][41]. They have mass dimension five and charge conjugation parity +1. These properties ensure that they do not mix with any other gauge-invariant operators of the same or lower dimension. Therefore, they are multiplicatively renormalizable. We take the corresponding renormalization factor from Table XI of Ref. [42].
For the evaluation of the twist-3 matrix element d 2 we use multiplets of operators with charge conjugation parity +1 which transform under the hypercubic group according to the representation τ The lattice operators (25) are Euclidean counterparts of the Minkowski operators (17), while the operators (27) correspond to the operators (18) with the field strength expressed in terms of a commutator of two covariant derivatives. The operators (26) are analogous to (19). Under renormalization all three multiplets are expected to mix with each other. In the continuum, R m disappears in the chiral limit. On the lattice, the explicit breaking of chiral symmetry caused by Wilson-type fermions persists even for massless quarks. Therefore O which should be small. The same holds for lattice counterparts of R eq and R eq1 . Hence, in a first approximation we take into account only the multiplets (25) and (26). Their renormalization and mixing can be treated along the lines of Ref. [42], provided one multiplies the operators (26) with 1/a. Both operator multiplets then have dimension five.
The renormalized operators of the multiplet (25) are now given by Here we stick to the notation of Ref. [42], wherê Z mm (µ, a) denotes the renormalization and mixing matrix in the nonperturbative scheme used on the lattice. For this scheme we choose the RI -MOM scheme, i.e., the operators are taken at vanishing momentum. In order to suppress powerlike lattice artifacts as far as possible the external quark momenta are chosen as with the renormalization scale µ. Presently we cannot convert the coefficientsẐ mm (µ, a) (and hence the renormalized operators) to the MS scheme, because the required perturbative calculations in the continuum are not yet available. Our procedure accounts for the mixing with lower-dimensional operators caused by the explicit breakdown of chiral symmetry in our simulations, but further mixing effects are still neglected.
In the chiral limit, the matrix element d 2 is multiplicatively renormalizable [31]. Rewriting Eq. (28) as we see that O will have a multiplicative dependence on µ if the ratioẐ 12 (µ, a)/Ẑ 11 (µ, a) does not depend on µ. It turns out that this requirement is better fulfilled when we use instead of the 2×2 matrixẐ(µ, a) the matrix Z(µ, a) constructed in the following way, cf. Ref. [43]. We computeẐ(µ, a)Ẑ −1 (µ 0 , a) for the renormalization scales µ of interest and a reference scale µ 0 chosen as µ 0 = 2 GeV. Within our approximations, this matrix should have a continuum limit, which we evaluate by fitting the lattice spacing dependence with a quadratic polynomial in a 2 . Denoting the result R(µ, µ 0 ), we definẽ To improve on this would require the consideration of quark-quark-gluon matrix elements instead of quarkquark matrix elements.

A. Lattice setup
To compute the reduced matrix elements a (10) and (11) for n = 2 we analyze a subset of the lattice gauge ensembles generated within the Coordinated Lattice Simulations (CLS) effort [44]. The ensembles have been generated using a tree-level Symanzik improved gauge action with N f = 2+1 flavors of nonperturbatively O(a) improved Wilson (clover) fermions. Near zero modes of the Wilson-Dirac operator are avoided by applying twisted-mass determinant reweighting to achieve stable Monte Carlo sampling [45]. Furthermore we improve the overlap of the interpolating currents at the source/sink time slice using Wuppertal smeared quarks [46] in the source/sink interpolators with APE smoothed spatial gauge links [47].
In most of our simulations we use open boundary conditions in time. Especially for the very fine lattices this avoids freezing of the topological charge and large autocorrelation times [45,48]. In order to suppress the distortions caused by the loss of translation invariance in time we restrict our measurements to regions with sufficiently large distances from the temporal boundaries, see, e.g., Refs. [42,49]. Only a few of the coarser lattices have been simulated using (anti)periodic boundary conditions. An overview of the gauge ensembles used in this work is given in Table I. They have been generated along three different trajectories in the quark mass plane, which are indicated in the last column of the table. Along the trajectory labeled by 'trm', the trace of the quark mass matrix is held constant, approximately equal to its physical value [44]. Along the trajectory labeled by 'msc', the renormalized strange quark mass is set to its physical value [50], and the symmetric line with equal masses of the light quarks and the strange quark is labeled by 'symm'. A general explanation of this strategy can be found in [50]. In summary we use six different lattice spacings ranging from about 0.039 fm up to 0.098 fm and m π goes down from ∼ 420 MeV to ∼ 220 MeV. With linear spatial lattice extents Lm π = aN s m π between 3.8 and 6.4, finite volume effects are expected to be moderate. Removing the leading O(a) discretization effects in the relevant matrix elements, requires Symanzik improvement of the corresponding operators. This has not been implemented in our study.
The extraction of the reduced matrix elements a  I. CLS gauge ensembles analyzed in this work. The ensemble identifier is given in the first column, followed by the inverse gauge coupling β and the lattice size. The column 'bc' indicates whether the boundary condition in time was open (o) or (anti)periodic (p). In the next columns we give the lattice spacing a, the pion mass mπ and the product of the spatial lattice size L = aNs with the pion mass. The column 't/a' contains the list of source-sink distances analyzed on this lattice. The subscript #meas specifies how many measurements have been performed for the respective source-sink distance. In physical units these distances roughly correspond to 0.9 fm, 1.0 fm and 1.2 fm. The number of analyzed configurations is given in the column 'n cnfgs ', and 'Traj.' specifies the trajectory in the quark mass plane to which the ensemble belongs. source method [59], extensively applying the so-called coherent sink method used by the LHPC Collaboration in Ref. [60]. All the computations are performed using the Chroma software package [61] and additional libraries implemented by our group.

B. Correlation functions
In Sec. II we used the OPE to relate the moments of the structure functions g i (x, Q 2 ) to the reduced matrix elements a  10) and (11). The corresponding matrix elements are extracted on the lattice from two-and three-point functions of the form The initial (final) momentum is denoted by p ( p ). The quantities of interest in this work allow us to restrict the kinematics to the forward limit, thus we use p = p from now on. The nucleon is created by the interpolating current N at the source time slice t src = 0 and annihilated at the sink time slice t. In the case of the three-point correlation function an additional local current O is inserted at the time slice τ with t > τ > 0. The nucleon interpolating current is defined by where C is the charge conjugation matrix and the quark fields are smeared separately in all spatial directions using the techniques mentioned in Sec. IV A. Furthermore, we define the positive parity projector P + = (1 + γ 4 )/2, and Γ = P + (−i)γ j γ 5 corresponds to the difference between the two spin projections with respect to the direction j = 1, 2, 3. When evaluating the three-point functions we consider quark-line connected diagrams only. Calculating the quark-line disconnected diagrams is computationally very expensive, but probably only of secondary importance for the physical quantities such as the color Lorentz force on a quark in a nucleon. However, we should keep in mind that, strictly speaking, only flavor-nonsinglet quantities like d are free of quark-line disconnected contributions.
The correlation functions are related to matrix elements by inserting complete sets of energy eigenstates. In the limit of large Euclidean times t, τ and t − τ excited states are exponentially suppressed and the correlation functions can be approximated by the ground-state contribution, where |N p σ denotes a nucleon state with spin projection σ and momentum p. The overlap matrix elements can be written as in terms of the momentum and smearing-dependent overlap factors Z p and the nucleon spinor u p,β σ . Similarly, the matrix elements of O can be expressed in the form Using the spinor identity σ u p σ u p σ = E p γ 4 −i p· γ+m N , we rewrite (35) and (36) as where The relations between the ground-state matrix elements N p σ |O|N p σ and the reduced matrix elements are given in Sec. II. However, in addition to the ground-state contributions we have to take into account possible excited states in Eqs. (39) and (40). An analysis of the first excited-state contribution is given in the next subsection.

C. Excited-state contributions
In the two-and three-point functions (39) and (40), respectively, the signal-to-noise ratio decreases exponentially with the source-sink separation in time. However, for small time distances between the operators we still find significant excited-state contributions. We take these contributions into account by including excitedstate terms in the spectral decomposition of the correlation functions. Our ansatz reads Here ∆E p denotes the energy difference between the first excited state and the ground state, which is taken to be the same in both correlators. The amplitudes of the excited-state contributions in the two-and three-point functions are denoted by A and B, respectively. All amplitudes depend on the smearing and on the momentum of the interpolating currents at the source and the sink, while B 10 , B 01 and B 11 also depend on the inserted current O and the spin matrix Γ. Since we only consider the forward limit, we may set B 10 = B 01 .

D. Ratios
Instead of performing a simultaneous fit to the twoand three-point functions, we consider the two-point functions along with ratios of three-point functions divided by two-point functions: In such a ratio the leading-order time dependence and the overlap factors are eliminated and the ground-state contribution corresponds directly to the matrix element we are interested in. Taking into account excited-state contributions according to Eqs. (42) and (43), we would arrive at the fit ansatz We assume that the ground-state energies are well described by the continuum dispersion relation in our fitting analysis. This need not be the case for ∆E p since in general multihadron states may contribute, e.g., N π and N ππ states.
Unfortunately, our data do not allow us to determine B 11 . Therefore we omit this term as well as the analogous contribution ∝ A in the denominator from our fit function for the ratio R p O,Γ . However, when performing a simultaneous fit to the ratio (44) and the two-point function (42) to extract the reduced matrix elements, the excited-state contribution ∝ A is taken into account in the two-point function. To fix the nucleon mass in the fits we include additional two-point correlators for p = 0. The fit range of τ is restricted to the interval 2a < τ < t − 2a resulting in reasonable values of χ 2 /dof. We choose on-axis momenta p = ± e i 2π/L with i taken to be different from the polarization direction of the nucleon, which is determined by Γ. The final analysis utilizes the data for all available momenta, nucleon polarizations and source-sink distances.
As an example we show in Fig. 1 the ratio (44) for the bare operators (25) (averaged over all members of the multiplet) with flavor u on the ensemble J304. The curves represent our fit, the horizontal line and the corresponding error band show the result for the ground-state contribution.

V. RESULTS
We present our results for the light flavors u and d separately, where quark-line disconnected contributions have been neglected. However, we consider it to be unlikely that the latter would modify our numbers beyond the size of the other uncertainties. Superscripts (u) and (d) always refer to the u and d quarks in the proton, while (p) and (n) denote the matrix elements (14) for the proton and the neutron, respectively, Alternatively, one can write where d does not suffer from the omission of disconnected diagrams.
Approximate renormalization in the RI -MOM scheme, as described in Sec. III, is performed at some scale µ . We attempt a combined continuum and chiral extrapolation using the fit functions (53) -(55) given below and evolve the results to the scale µ = 2 GeV with the help of the one-loop formula for the flavor-nonsinglet operators, where with N c = 3 and N f = 3. As we neglect disconnected contributions, we use this value of B, which is strictly speaking only correct for the (u−d) part, also for (u+d).
Varying the intermediate scale µ should give us some measure of the uncertainty related to the renormalization. The central results, however, are all obtained at µ = µ, so they do not depend on the perturbative value of B. When constructing our fit functions we take into account that the leading discretization effects in the matrix elements in our simulations are O(a). For the continuum and chiral extrapolation of the d 2 and a 2 data, we consider the fit formulas where The difference between the corresponding results should provide an impression of the uncertainty inherent in the fit procedure. The gauge field ensembles used in the fits are collected in Table I. Note that on a few configurations in some of the coarser ensembles (H105, H106, H107, N450, B452 and N201) we encountered measurements that were separated by more than a hundred 68% confidence level intervals from the results obtained on the remaining configurations. The origin of these deviations is unclear and we excluded these configurations from further analysis. The fit results for d 2 are given in Table II for the case of the renormalization matrixZ, cf. Sec. III. The dependence on the intermediate scale µ is quite weak, due to the use ofZ. However, the dependence on the choice of the fit function is more pronounced, although all three functions yield reasonable fits with χ 2 /dof between 0.66 and 0.94. In Fig. 2 we plot our data along with the fit function (53) for µ = 2 GeV versus the lattice spacing a. The legend identifying the ensembles is given in Fig. 3. The curve shows the fit function evaluated with the physical values of m π and m K . The fitted coefficients C 3 and C 4 have been used to shift the data points vertically such that they correspond to the physical masses.
In order to visualize the mass dependence we show in Fig. 4 the results for d  Table II plotted against m π . The curve represents the function f 1 (0, m π , m phys K ) evaluated with the fitted parameters, and the data points have been shifted by subtracting f 1 (a, m π , m K ) − f 1 (0, m π , m phys K ). Compared to the a dependence, the dependence on m π turns out to be more moderate.
Results for the twist-2 matrix element a 2 are presented in Table III. The flavor dependence is of the same form as in the case of d 2 . Since the operators that we use for the determination of a 2 are multiplicatively renormalizable, we can follow the standard RI -MOM procedure to obtain values in the MS scheme at µ = 2 GeV. As in the case of d 2 , we obtain reasonable fits with all three fit functions (0.80 ≤ χ 2 /dof ≤ 0.92), but find some dependence of the results for a 2 on the fit function. Plots of our a 2 data and the fit function (53) are shown in Fig. 5, which is analogous to Fig. 2. Again, the dependence on the pion mass appears to be rather weak.
For our final values, collected in Table IV, we take the results given in the first line of Table II for d 2 and  Table III for a 2 . The errors in these tables are purely statistical, but there are several sources of systematic uncertainties. Comparing the results obtained by varying the fit function or the scale µ allows us to estimate the influence of the extrapolation method. Since the renormalization in the case of a 2 is less subtle than in the case of d 2 , we have refrained from varying the intermediate renormalization scale µ in the analysis for a 2 . As the estimate of the systematic error due to our extrapolation we take the maximum of the (absolute value of the) difference between the final value and the results obtained by means of the fit functions f 2 and f 3 . Unfortunately, we are not able to assess the effect of neglecting the operators (27) in the process of the renormalization of d 2 . This problem must be left for future investigations. Another source of error that cannot yet be quantified is the 2 . This extrapolation corresponds to the first line of Table II. Results belonging to the same value of a are horizontally shifted for visibility by a small amount such that the exact value of a lies in the center of the respective group. Within each of these groups the pion mass decreases from left to right. The results from the ensembles A654, H105 and D451 are not shown on the plots because of their large error bars.
omission of quark-line disconnected contributions, which leaves only a unaffected. However, we expect that this error is small compared to the other uncertainties.

VI. COMPARISON WITH EXPERIMENT
Several experiments have measured the structure functions g i (x, Q 2 ) for certain ranges of the variables x and Q 2 and attempted to determine the moments for the proton as well as for the neutron. Results are given for Q 2 up to 5 GeV 2 . As far as we can see, the Wilson coefficients are taken into account only in leading order. In this approximation the moments of the structure functions are related to the reduced matrix elements by and with µ 2 = Q 2 .  Final results for a2 (MS scheme) and for d2 (renormalizationZ), both at µ = 2 GeV. The first error is statistical, while the second error accounts for the uncertainty due to the combined chiral and continuum extrapolations. The error caused by the approximations in the renormalization procedure for d2 cannot be quantified. Let us begin with the results for a 2 . In Table V we collect the values given for 1 0 dx x 2 g 1 (x, Q 2 ) in the literature with statistical and systematic errors (if they are given separately) added in quadrature.
In order to compare the results given in Table V with our numbers we take into account the one-loop QCD corrections, which are flavor independent, i.e., we divide the experimental values of 1 0 dx x 2 g 1 (x, Q 2 ) by the Wilson coefficient (9) with µ 2 = Q 2 . Then we use the renormalization group to evolve the renormalization scale of the resulting matrix element a 2 to our value µ 2 = 4 GeV 2 . We employ the five-loop anomalous dimension of the relevant operator multiplet along with the five-loop β function. The details of the calculation are the same as in 2 (µ) Ref. [13] 5.0 0.0497(40) −0.0096(64) Ref. [62] 4.2 0.0427(32) -Ref. [62] 5.0 0.0342(70) -Ref. [18] 3.21 -0.0032(24) Ref. [18] 4.32 -0.0020 (25) Ref. [42]. The resulting values for a 2 (2 GeV) can be found in Table VI.
In Table VII we collect the values presented for d alt 2 = d 2 /2 in the literature with statistical and systematic errors (if they are given separately) added in quadrature. As in this case Wilson coefficients beyond tree level are not available, we just use the renormalization group to evolve the renormalization scale to our value µ 2 = 4 GeV 2 . The corresponding factors are calculated from Eqs. (51) and (52), and the resulting values for d 2 (2 GeV) can be found in Table VIII.
As in the analysis of the experimental data the variation of Q 2 in the respective experimental setup generally was not taken into account, applying the renormalization group running is not too well justified, but the effect is anyhow quite small compared to the experimental uncertainties. If, however, the sign change from negative numbers at the smaller values of Q 2 to positive numbers at larger Q 2 is taken seriously and interpreted as a "nontrivial scale dependence" [19], the perturbative renormalization group would not be applicable in this range of Q 2 .
In addition to these results from single experiments, there is also a global analysis of polarized inclusive deepinelastic scattering available [63]. Unfortunately, the resulting values for d   [12] we have chosen the "SLAC average".

VII. COMPARISON WITH OTHER LATTICE DETERMINATIONS
There are a few previous lattice investigations of a 2 and d 2 , with which we can compare our new results. For this purpose we consider the reduced matrix elements a (q) 2 and d (q) 2 with q = u, d. In Ref. [11] a continuum limit was not attempted. Instead, the results on the finest lattice in the chiral limit were considered as the best estimates obtained in this N f = 2 simulation. With the help of the perturbative running factor we evolve our a 2 results from µ 2 = 4 GeV 2 to the scale µ 2 = 5 GeV 2 used in Ref. [11]. This yields a  (50) given in Ref. [11]. The statistical errors of our present determination are significantly smaller than those quoted in Ref. [11], while the central values are in rough agreement with each other. Unfortunately, the uncertainties due to the combined chiral and continuum extrapolations, which could not be estimated in the previous study, are still rather large.
Values for x 2 ∆u = a (u) 2 /2 and x 2 ∆d = a (d) 2 /2 from another N f = 2 simulation are presented in Ref. [66]. They are obtained at a single lattice spacing a ≈ 0.1 fm with the help of perturbative renormalization. Since the results are given at µ 2 = 4 GeV 2 in the MS scheme, they can directly be compared with our numbers. From Ta is contained in Ref. [60]. This paper relies on simulations with N f = 2 + 1 flavors at a single lattice spacing a = 0.124 fm and employs a combination of perturbative and nonperturbative methods for the renormalization. Results are given at a renormalization scale µ 2 ≈ 2.5 GeV 2 for the generalized form factorsÃ u−d 10 (t) andÃ u−d 30 (t), which are related to x 2 ∆u−∆d = a  Table IV.
Results from the quasidistribution approach [67] lead to the conclusion that the Wandzura-Wilczek approximation for g T (x) = g 1 (x) + g 2 (x) works well at least up to x ≈ 0.4. Whether the deviations observed at higher values of x indicate nonvanishing twist-3 effects or have a different origin, remains to be seen.

VIII. CONCLUSIONS
In the present paper we computed the nucleon matrix elements a 2 and d 2 , which determine the x 2 moments of the spin structure functions g 1 and g 2 , in lattice QCD, thus improving on the earlier evaluation [11] (coauthored by some of us). In both determinations quark-line disconnected contributions were neglected.
For the twist-2 matrix element a 2 we found a (p) 2 = 0.069(5)(16) and a (n) 2 = 0.0068 (33)(82) at the renormalization scale µ = 2 GeV, see also Table IV. In both cases the second (systematic) error is considerably larger than the first (statistical) error. Our results are in broad agreement with phenomenology, as shown in Figs. 6 and 7.
The parameter d 2 quantifies a specific twist-3 quarkgluon correlation in the nucleon. Because it is experimentally accessible it became a much discussed test case for our understanding of hadron structure beyond twist 2. We observed a strong dependence on lattice spacing for d Ref. [13] Ref. [61] Ref. [61] this work Ref. [13] Ref. [18] Ref. [18] this work Ref. [13] Ref. [14] Ref. [61] Ref. [61] Ref. [16] Ref. [19] Ref. [19] this work  Ref. [12] Ref. [13] Ref. [14] Ref. [15] Ref. [17] Ref. [17] this work probably was somewhat accidental. In contrast, in this new lattice determination of d we take the lattice spacing dependence into account and only when doing this, the results agree well with experiment (and, therefore, also with the numbers given in Ref. [11]). Note also that the a dependence of d (n) 2 is far less pronounced. These results provide a showcase example justifying the CLS strategy to focus its resources on controlling the continuum limit. The a dependence of any observable of interest can be strong (as for d 2 ). What is the case has to be carefully evaluated for each specific quantity.
The final results can be found in Table IV. We obtained d With m 2 N ≈ 4.47 GeV/fm we obtain where we have added the two errors in quadrature.