New insights on proton structure from lattice QCD: the twist-3 parton distribution function $g_T(x)$

In this work, we present the first-ever calculation of the isovector flavor combination of the twist-3 parton distribution function $g_T(x)$ for the proton from lattice QCD. We use an ensemble of gauge configurations with two degenerate light, a strange and a charm quark ($N_f=2+1+1$) of maximally twisted mass fermions with a clover improvement. The lattice has a spatial extent of 3 fm, lattice spacing of 0.093 fm, and reproduces a pion mass of $260$ MeV. We use the quasi-distribution approach and employ three values of the proton momentum boost, 0.83 GeV, 1.25 GeV and 1.67 GeV. We use a source-sink separation of 1.12 fm to suppress excited-states contamination. The lattice data are renormalized non-perturbatively. We calculate the appropriate matching within Large Momentum Effective Theory, which is applied to the lattice data in order to obtain $g_T$. The final distribution is presented in the $\overline{\rm MS}$ scheme at a scale of 2 GeV. We also calculate the helicity distribution $g_1$ to test the Wandzura-Wilczek approximation for $g_T$. We find that the approximation works well for a large range in $x$. This work demonstrates the feasibility of accessing twist-3 parton distribution functions from novel methods within lattice QCD, and can provide important insights on the structure of hadrons.

Introduction: Parton distribution functions (PDFs) are fundamental quantities characterizing the quark and gluon structure of hadrons such as the proton [1]. They can be classified according to their twist, which describes the order in 1/Q 2 , where Q is the momentum transfer, at which the relevant operator is suppressed in a hard scattering process. The leading-twist (twist-2) PDFs are O(1) in 1/Q 2 and thus, they describe the most important contributions. They can be considered probability densities for finding, inside the hadron, a parton which carries the fraction x of the hadron's momentum. Twist-2 PDFs have so far attracted most of the attention. But highertwist (such as twist-3) distributions are also important, as they are not necessarily smaller than twist-2 PDFs. While they don't have a density interpretation, twist-3 PDFs contain information about quark-gluon-quark correlations [2,3], and as such characterize the structure of hadrons in a new way.
Twist-3 PDFs appear in QCD factorization theorems for a variety of hard scattering processes, and have interesting connections with transverse momentum dependent parton distributions, thus offering important insights into the latter, see, e.g., Refs. [4][5][6]. Additionally, the extraction of the g 1 (x) twist-2 structure function is dependent of the twist-3 g 2 (x) structure function. Indeed, recent measurements were reported in Refs. [7,8] for the matrix elements of g 2 (x). However, pinning down the corre-sponding twist-3 PDFs from experiment is generally difficult, because the related observables are often kinematically suppressed. Nevertheless, measurements related to twist-3 PDFs are part of the Jefferson Laboratory ongoing 12 GeV program and will be important also for the planned Electron-Ion Collider [9,10].
In this work, we present the first-ever calculation of the isovector flavor combination of g u−d x) from first principles using lattice QCD. For simplicity in the notation, we omit the superscript u − d in the remainder of this paper. We make use of the socalled quasi-PDF approach suggested by X. Ji [11,12]. While standard (light-cone) PDFs are defined through light-cone correlation functions, quasi-PDFs and related quantities [13][14][15] are given by spatial correlation functions accessible in lattice QCD, which gave rise to an intensive surge of studies, see, e.g., Refs.  and the recent reviews of Refs. [45,46]. Because quasi-PDFs and light-cone PDFs have the same infrared physics [15,[47][48][49], they can be related using perturbative QCD, in a procedure called matching, see Refs. [31,35,[50][51][52][53] for related recent work. The matching equations are known only for twist-2 operators, and within this work we extract the one-loop matching kernel for g T .
Our first-principles calculation also allows us to address the validity of the Wandzura-Wilczek (WW) approximation for g T (x) [54]. Generally, the Mellin moments of g T (x) receive contributions from twist-2 operators and twist-3 operators (whose moments we denote by d n ). Therefore, g T (x) can be written as g W W T (x) + g twist−3 T (x), with g twist−3 T (x) the contribution from twist-3 operators [6]. In the WW approximation, one sets d n = 0, implying that g T (x) is fully determined by the twist-2 operators which define g 1 (x). Thus, the study of the WW approximation gives direct information about the importance of twist-3 operators. Here, we present the first check in lattice QCD of how relevant the twist-3 operators are for the x-dependence of g T (x).
Methodology: The calculation is based on matrix elements of a non-local operator, with space-like separated fermion fields, which are connected via a straight Wilson line (WL) of length z. The operator has a Dirac structure γ j γ 5 , and the matrix element is defined in position (z) space as The proton is boosted in a spatial direction, and the quasi-distributions approach requires that it is in the same direction as the WL, i.e. P = (iE, 0, 0, P 3 ). To obtain the twist-3 distribution, γ j must be γ x or γ y , each requiring a parity projector (1+γ 0 )iγ 5 γ j /4. In this work, we average over the two operators to increase the statistical accuracy.
For the proper evaluation of g T , one must extract the ground-state contribution from M g T . This is achieved by a large time separation between the initial (source) and final (sink) state of the proton, T sink , as well as by a current insertion which is away from the source and the sink. Once these conditions are satisfied, we identify the ground state using a fit to a constant (plateau region). The desired quantity, F g T , is extracted based on the continuum decomposition: in Euclidean space. The kinematic factor is obtained based on the normalization conventions on the lattice. E is the energy of the proton, E = m 2 + P 2 3 , m is its mass, and Z g T is the renormalization function, and it is also calculated in this work.
The so-called quasi-distribution, g T , is defined as the Fourier transform of F g T (P 3 , z) over z. It is, thus, given in the momentum representation, x, where Λ∼1/a is a UV cut-off. Our definition of g T is such that its lowest x-moment is independent of P 3 , see also Ref. [55]. As the momentum P 3 increases, the quasidistribution g T can be matched to the light-cone distribution g T using a perturbative formula obtained within Large Momentum Effective Theory [11,12].
Computational setup: In this work, we use one N f = 2 + 1 + 1 ensemble of twisted mass fermions [56,57] with clover improvement [58] and Iwasaki gluons [59]. The lattice spacing is 0.093 fm, its volume is 32 3 ×64 (L ≈ 3 fm), and the pion mass is around 260 MeV. We focus on the isovector combination, which receives contributions only from the connected diagram. T sink is taken to be above 1 fm (T sink = 1.12 fm), for which excited-states contamination is assumed to be suppressed for the values of P 3 we employ [35]. We apply stout smearing [60] to the links of the operator, which is known to reduce statistical uncertainties in matrix elements of non-local [35], and gluonic [61,62] operators. In this work, M g T is calculated for three values of momentum boost, P 3 = 0.83, 1.25, 1.67 GeV. The statistical uncertainties increase with the momentum, and therefore, the number of measurements must increase by around one order of magnitude with each additional momentum unit to achieve similar statistical errors. We use the momentum smearing method [63] on the proton interpolating field, which offers a better overlap of the interpolator and the ground state. The momentum smearing parameter has been tuned following the procedure described in [35], and leads to a significant reduction of statistical uncertainties. We achieve similar accuracy for each boost with 1552, 11696, and 105216 measurements at P 3 = 0.83, 1.25, and 1.67 GeV, respectively. Using the same ensemble, simulation parameters, and statistics, we also obtain the leading-twist helicity PDF, g 1 (x).
The dependence of the bare F g T on P 3 is shown in Fig. 1 for each z/a value. The top (bottom) plot shows the real (imaginary) part for z ≥ 0 and the real (imag-inary) part is an even (odd) function of z. For both parts, we find convergence between the two largest momenta, within statistical uncertainties. However, there is no guarantee that the convergence will hold for g T (x), as the matching formula depends on P 3 .

Renormalization:
One of the crucial aspects of the calculation is the renormalization of the bare matrix elements. Non-local operators containing a WL require an evolved renormalization procedure, in contrast to the local fermion operators. The presence of the WL is associated with a power divergence with respect to the lattice spacing [64,65]. Such divergence must be removed along with all logarithmic divergences so that physical meaning can be attributed to lattice data. Considering that this is the first study of twist-3 distributions, we focus on the extraction of the matrix elements, and we do not take into account any mixing with other twist-3 operators (e.g., with quark-gluon-quark operators). Estimating such mixing would require further evolution of lattice calculations for twist-3 operators. Therefore, the renormalization is taken to be multiplicative.
We employ the renormalization procedure developed for straight-WL non-local operators [21,66], and is also used for twist-2 distributions (see, e.g., Ref. [35]). We calculate non-perturbatively the renormalization functions in the RI scheme [67] at each value of z/a separately. We use a set of five N f = 4 ensembles [68] produced specifically for the renormalization functions of the ensemble used in this work. We eliminate possible systematic uncertainties in the renormalization functions by an advanced program, in which we: a. perform a chiral extrapolation on the five ensembles, b. use a wide range of RI renormalization scales and fit the MS estimates to eliminate any dependence on the initial scale, c. remove discretization effects utilizing results in lattice perturbation theory [69]. The renormalization factors are complex functions due to the presence of the WL and thus, the renormalized matrix elements are obtained from the complex multiplication Z MS g T (z) · M g T (P, z). The renor-malized matrix elements are given in the modified MSscheme (MMS) [35] at the scale of 2 GeV.
Reconstruction of x-dependence: The lattice calculation provides determinations of F g T (P 3 , z) for discrete values of z ≤ z max , with z max ∼ L/2. Thus, Eq. (3) needs to be discretized and becomes subject to an ill-defined inverse problem, as discussed in Ref. [30]. One of the methods advocated to solve this problem is the Backus-Gilbert method [70], which maximizes the stability of the solution with respect to the statistical variation of the data. Thus, it provides a model-independent assumption allowing one to obtain a unique reconstructed quasi-distribution from the available set of matrix element evaluations. We employ the Backus-Gilbert method for the results presented here.
Matching to light-cone g T (x): Another novel aspect of this work is the calculation of the matching formula, which connects g T (x) to the light-cone g T (x). C is the matching kernel, which is calculated within one-loop perturbation theory in momentum space. The matching kernel is not available in the literature for twist-3 distributions, and is calculated here for g T for the first time. We explore two schemes for the matching, which use the same bare matrix elements, but the renormalization functions are converted to different schemes. The first one uses Z MS , while the second one uses the so-called modified MS-scheme, Z MMS . The latter preserves the normalization, unlike the MS scheme, through an extra renormalization in the "unphysical" |ξ| > 1 region. More details on the perturbative calculation for the matching kernel can be found in a follow-up publication [71]. For the results presented here, we employ the MMS scheme for the quasi-distributions, for which C takes the form: Note that the light-cone g T (x) results are always in the MS-scheme, regardless of the scheme used for g T (x).
Results on g T (x): The various steps described above are combined to provide the final estimates for the twist-3 distribution g T (x). The renormalized ground-state contributions to the matrix elements, F g T , are transformed to x-space using the Backus-Gilbert method, and then matched using Eq. (4). In Fig. 2, we plot the dependence of g T (x) on the proton momentum P 3 for the quark (x > 0), and anti-quark (x < 0) regions. With red, green and blue bands we show the distributions for P 3 = 0.83, 1.25, and 1.67 GeV, respectively. The widths of the bands represent statistical uncertainties. We find that the distribution in the region x < 0.4 becomes slightly more narrow as the momentum increases. In addition, we find convergence between the two largest momenta. Regarding anti-quarks, we observe similar functional forms for these momenta. Since g T is a sub-leading contribution, it is interesting to compare it with the leading-twist PDF. To this end, we calculate the twist-2 helicity PDF g 1 (x) on the same ensemble and using the same values for P 3 and T sink . In Fig. 3, we show the results for the highest momentum, P 3 = 1.67 GeV. It is interesting to observe that g 1 and g T mostly differ in the positive-x region. While g 1 has a smaller value as compared to g T in the low-x region, it has a much smaller slope at x ≈ 0.1 − 0.3. As a consequence, it becomes dominant in the region x 0.2. The two distributions are in agreement in the anti-quark region within uncertainties.
The light-cone PDFs g T and g 1 (and the corresponding quasi-PDFs) are connected via the Burkhardt-Cottingham sum rule [72], which serves as an important check for the lattice data. We find, without any additional input, that Eq. (6) gives 0.01 (20) and therefore the sum rule is satisfied.

Wandzura-Wilczek approximation:
The lattice data of this work for g T and g 1 may be used to test the WW approximation [54]. We present here, for the first time, a check of the full x-dependence of the WW approximation in lattice QCD. As already discussed above, in this approximation, the twist-3 g T (x) is fully determined by the twist-2 g 1 (x) (and denoted by g WW T ) , We evaluate g WW T using the lattice data for a wide range of x. The resulting x-dependence can be compared to the data for g T (x), as shown in Fig. 4 for P 3 = 1.67 GeV. The focus is on the quark region (x > 0), which is less susceptible to systematic uncertainties, as compared to the anti-quark region. We find that for a considerable x-range, our numerical results for g T (x) are consistent with g WW T (x). However, given the uncertainties of the final distributions, a violation of the WW approximation is still possible at the level of up to 40% for x 0.4. Interestingly, a check of the WW approximation based on experimental data leads to a similar possible violation at the level of 15 − 40%, depending on x [6]. It is also notable to mention that while the slopes of g T and g 1 differ (see Fig. 3), the slopes of g T and g WW T are the same up to x ≈ 0.4. The difference of g T and g WW T for large x could be due to systematic uncertainties, yet to be investigated. However, it may also indicate larger violations of the WW approximation in this region.
As an additional consistency check, we calculate the r.h.s. of Eq. (7) using g 1 from global fits by the NNPDF [73] and JAM [74] collaborations. We find good agreement with lattice-extracted g WW T up to x ≈ 0.3. Above this x value, the discrepancy again indicates possible systematic effects.

Summary and prospects:
We presented a pioneering ab initio calculation of the proton twist-3 distribution g T (x), using numerical sim- ulations of lattice QCD, within the quasi-distribution method. The work comprised multiple non-trivial steps, i.e. calculation of matrix elements of fast-moving protons and non-local operators in position space, elimination of divergences, reconstruction of the x-dependence, as well as matching to the light-cone distribution. For the reconstruction of the quasi-distribution, we used the Backus-Gilbert method, which improves the results by providing a unique solution to the inverse problem. Another novel result of this work is the calculation of the matching kernel for the case of g T . We emphasize that obtaining the matching formula for twist-3 operators requires a very delicate process as compared to the twist-2 case, and details will be presented in a forthcoming article [71].
The light-cone g T was obtained for three values of the momentum boost, P 3 = 0.83, 1.25, 1.67 GeV, and is presented in Fig. 2. We found that g T decays much faster than the leading-twist g 1 , as shown in Fig. 3. One important aspect of this work is the implementation of the Wandzura-Wilczek approximation, using both lattice data and data from global fits. We find g T consistent with its WW approximation for x 0.4, but within uncertainties one can not exclude its violation at the level of up to 40%, which is consistent with earlier studies based on experimental data. A possibly larger violation is conceivable at larger x according to our results. Nevertheless, careful investigation of systematic uncertainties is needed for more precise quantitative statements, particularly at high-x. The role of systematics in this region is confirmed by our consistency check comparing lattice-extracted g WW T with the ones from global fits, where agreement is observed for x 0.3.
We are considering several directions to extend this calculation. Detailed investigations are required to quantify systematic uncertainties, such as excited states con-tamination, reconstruction of x-dependence, finite volume and discretization effects. The latter two require a minimum of two and three ensembles, respectively. Also, simulations with quark masses fixed to their physical values (physical point) are now feasible within the available computational resources. We will move in that direction once systematic uncertainties for twist-3 distributions are understood.
Finally, the possible breaking of the WW approximation at large x signals a sizeable contribution from the d n terms. The connection between lattice estimates and results from experiments and phenomenology is more immediate for these quantities, because of recent measurements (see, e.g., Ref. [7]) of d 2 . The latter has a semiclassical interpretation of the average transverse force acting on the struck quark in a transversely polarized proton in DIS, right after it has been hit by the virtual photon [75]. Our work shows that such a caliber lattice calculation is not in a distant future, and a systematic study of the x-dependence of the twist-3 contributions in the proton, which is currently unknown both experimentally and theoretically, is within reach. de-sc0020405. Computations for this work were carried out in part on facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy. This research was supported in part by PLGrid Infrastructure (Prometheus supercomputer at AGH Cyfronet in Cracow). Computations were also partially performed at the Poznan Supercomputing and Networking Center (Eagle supercomputer), the Interdisciplinary Centre for Mathematical and Computational Modelling of the Warsaw University (Okeanos supercomputer) and at the Academic Computer Centre in Gdańsk (Tryton supercomputer). The gauge configurations have been generated by the Extended Twisted Mass Collaboration on the KNL (A2) Partition of Marconi at CINECA, through the Prace project Pra13 3304 "SIM-PHYS".