Twist-3 Generalized Parton Distribution for the Proton from Basis Light-Front Quantization

We investigate the twist-3 generalized parton distributions (GPDs) for the valence quarks of the proton within the basis light-front quantization (BLFQ) framework. We first solve for the mass spectra and light-front waved functions (LFWFs) in the leading Fock sector using an effective Hamiltonian. Using the LFWFs we then calculate the twist-3 GPDs via the overlap representation. By taking the forward limit, we also get the twist-3 parton distribution functions (PDFs), and discuss their properties. Our prediction for the twist-3 scalar PDF agrees well with the CLAS experimental extractions.

The GPDs are functions of the longitudinal momentum fraction carried by the struck parton (x) [6,[22][23][24][25], the longitudinal momentum transfer, known as skewness (ξ), and the square of total momentum transfer to the hadrons (−t), providing a three-dimensional picture of the hadron structure.The information encoded in the GPDs is richer than the ordinary parton distribution functions (PDFs) since PDFs only contain onedimensional information.In the forward limit, the GPDs reduce to the ordinary PDFs.The GPDs can also be connected with the form factors (FFs), the orbital angular momentum (OAM), the charge distributions, etc., [6,[26][27][28].
Most of the GPD-related studies concentrate on the leading twist [1,29,30].The GPDs at sub-leading twist are suppressed and less known, but in fact these are not negligible.There are several reasons for the importance of the twist-3 GPDs.First, there are relations between the quark orbital angular momentum inside the nucleon and twist-3 GPDs [6,26,31].Second, some studies show that from twist-3 GPDs, we can obtain information about the average transverse color Lorentz force acting on quarks [32,33].Third, the twist-3 DVCS amplitude can be expressed in terms of the twist-3 GPDs through the Compton form factors [34].
In this work, we calculate the twist-3 GPDs of the proton using the light-front wave functions (LFWFs) obtained by diagonalizing an effective Hamiltonian of the proton within the BLFQ framework.The basis truncations and the Fock sector truncations are adopted to enable practical numerical calculations.We also present the twist-3 PDFs obtained by taking the forward limit of the twist-3 GPDs.
The organization of this work is as follows.We briefly introduce the BLFQ framework in Sec.II.Then we give a simple introduction to the GPDs in Sec.III, where we also show the overlap representation of the twist-3 GPDs.In Sec.IV, we present all the numerical results including the twist-3 GPDs and PDFs, and their properties including sum rules.At the end, we summarize this work and elaborate on future prospects in Sec.V.
The main idea of BLFQ is to obtain the mass spectrum and bound state wave functions simultaneously by solving the eigenvalue problem where P − is the light-front Hamiltonian of the system, P − ψ is the eigenvalue that represents the lightfront energy of the system, and |ψ⟩ is the eigenvector of the system that encodes the structural information of the bound state.The notation for an arbitrary four-vector in the light-cone coordinate that we adopt is a = (a + , a 1 , a 2 , a − ).The + and − components are defined by a ± = a 0 ± a 3 , and the transverse directions are ⃗ a ⊥ = (a 1 , a 2 ).
The invariant mass is related to the light-front energy according to where P + represents the longitudinal momentum and P ⊥ represents the transverse momentum of the system.BLFQ uses a Fock-space expansion for hadronic states.At fixed light-front time (x + = x 0 + x 3 = 0), one can expand the proton state as where • • • represents all other possible parton combinations that can be found inside the proton.This work takes only the leading Fock sector into consideration.Within each Fock space, the proton eigensolution |ψ(P, Λ)⟩ can be expanded in terms of n-parton states |p + i , ⃗ p ⊥ i , λ i ⟩ as: where is the transverse relative momentum, and λ i is the light-cone helicity for the i-th parton.Those n-parton states are normalized as In this work, the light-front effective Hamiltonian designed for the proton in leading Fock sector is [44] where m i is the mass of i-th constituent, the subscript i, j are the indexes for particles in the Fock sector, V conf i,j is the confining potential, and is the one-gluon exchange (OGE) interaction.The light-front wavefunction ψ n (x i , ⃗ k ⊥ i , λ i ) defined in Eq. ( 4) can be obtained from diagonalizing the light-front Hamiltonian.The specific form of the confining potential is [44,55] where κ defines the strength of the confining potential, ⃗ r ij⊥ = √ x i x j (⃗ r i⊥ − ⃗ r j⊥ ) is the relative coordinate.The schematic form of the OGE potential [55] where C F = −2/3 is the color factor, α s is the coupling constant, Q 2 ij = −q 2 is the average momentum transfer squared carried by the exchanged gluon, and u s (p) represents the spinor which is the solution of free Dirac equation with momentum p and spin s.The OGE interaction plays an important role in the dynamical spin structure of LFWFs [43], which allows us to calculate spin dependent observables.
For the longitudinal component of the basis state, we choose a plane-wave state confined in a longitudinal box as with anti-periodic boundary condition for fermions.2L is the length of the confining longitudinal box, and k is the quantum number that represents the longitudinal degree of freedom.The longitudinal momentum is given by p + = 2πk/2L, with k = {1/2, 3/2, 5/2, • • • }.Two-dimensional harmonic oscillator (2D-HO) states are adopted for the transverse components as is the associated Laguerre polynomial, and φ = arg(⃗ p ⊥ ).The last degree of freedom is the light-cone helicity state in the spin space represented by λ.Then we have a complete set of single-particle quantum numbers that represent a Fock-particle state, {k, n, m, λ}.These single-particle states are orthonormal.
Within the BLFQ framework, we introduce a Focksector truncation and two cutoffs for practical calculations.Only the first Fock sector |qqq⟩ is taken into consideration in this work.The two cutoffs are represented by N max and K. N max is the cutoff in the total energy of the 2D-HO basis states in the transverse direction, given by i (2n i + |m i | + 1) ≤ N max , and K = Σ i k i represents the cutoff in the longitudinal direction.
This work uses single-particle states to construct the BLFQ basis.It has an advantage for retaining the correct Fermion statistics for quarks.However, the many-particle basis therefore incorporates the transverse center-of-mass (CM) motion which is entangled with intrinsic motion.It is then necessary to add a constraint term into the effective Hamiltonian to enforce factorization of LFWFs into a product of internal motion and CM motion components.By removing the transverse CM motion component, one obtains a boost-invariant LFWF.
The CM motion is governed by Here λ L is the Lagrange multiplier, 2b 2 is the zero-point energy and I is a unity operator.By setting λ L sufficiently large, it is possible to shift the excited states of CM motion to higher energy and ensure that low-lying states are all in the ground state of CM motion.
The expansion of the proton state in terms of BLFQ many-particle basis states is (13) where the completeness |x, n, m, λ⟩⟨x, n, m, λ| = 1 is used, and the ψ Λi xi,ni,mi,λi = ⟨x i , n i , m i , λ i |P, Λ⟩ are the wave functions in BLFQ.

III. GENERALIZED PARTON DISTRIBUTION
In this section, we will present the definition of twist-3 GPDs, and their overlap representations using LFWFs.The GPDs are functions of three variables, x representing the longitudinal fraction, ξ for the skewness, and −t signifying the momentum transfer squared.The GPDs are defined as off-forward matrix elements of a bilocal operator as where P , P ′ represent the momenta of the initial and final proton respectively, Λ, Λ ′ represent the light-front helicities of the initial and final proton respectively, and W (y, x) ≡ P exp(ig is the gauge link that ensures that the bilocal operator remains gauge invariant.Since we are working in the light-cone gauge (where A + = 0), the gauge link is then unity.Γ is one of the sixteen Dirac gamma matrices.We choose a symmetric frame throughout this work: where P = (P + P ′ )/2 is the average momentum, ∆ is the momentum transfer, ξ = −∆ + /2 P + is the skewness, M is the proton mass, and t ≡ −∆ 2 is the momentum transfer squared.Usually there are two ξ-dependent regions in the DVCS process for quarks, one is −ξ < x < ξ, called the Efremov-Radyushkin-Brodsky-Lepage (ERBL) region, the other is ξ < |x| < 1, called the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) region.This work focuses on the zero-skewness limit, i.e. ξ = 0, so only the DGLAP region applies.The DGLAP region describes a quark scattered off the proton, absorbing the virtual photon and immediately radiating a real photon, then returning to form a recoiled proton.This is a n → n diagonal (parton-number conserved) process.
With the parameterization taken from Ref. [64], which includes both chiral-even and chiral-odd GPDs, one finds that the sixteen sub-leading twist GPDs are defined as where , and ϵ ij ⊥ = ϵ −+ij with the antisymmetric Levi-Civita tensor ϵ −+12 = 1.Here i, j can only be transverse indice 1, 2.
A different parameterization of GPDs with vector and axial vector is introduced in Ref. [65] F µ =ū(P ′ ) P µ γ F µ =ū(P ′ ) P µ γ + γ 5 By using relations based on the Dirac equation [65][66][67][68], the two types of GPDs above can actually be related to each other according to the equations in Appendix A. In addition, there is also another parameterization defined in Ref. [69], and their relations to the GPDs in this work can also be found in Appendix A. There are also some other parameterizations of GPDs [70][71][72][73], but we will not illustrate them here.

A. Twist-3 GPDs
With overlap representation and the LFWFs obtained from solving BLFQ with the specific implementations descibed above, we calculate all the twist-3 GPDs and display the results in Figs.
(1)- (5).These overlap expressions in Eqs. ( 28)- (35) are slightly different between i = 1 and i = 2, but we find that the numerical results are the same, which demonstrates rotational symmetry in the transverse plane within our approach.Thus we only show one of the two choices for both u and d quarks for simplicity.Since ξ is fixed at zero in this work, we plot the GPDs as functions of x and −t.In this zeroskewness limit, eight of the twist-3 GPDs, H 2T , E 2T , H2T , Ẽ′ 2T , H2 , H ′ 2 , E ′ 2 and Ẽ′ 2 , are consistent with zero in our calculations to within our numerical uncertainty.This conclusion is consistent with their properties that we will discuss later.
Ẽ2T is the only γ ⊥ structure GPD that survives in the zero-skewness limit, and its 3D structures are shown in Fig. 1.The distributions for both d and u quarks have peaks at small x and small −t region with opposite sign.This is due to the appearance of x and ∆ in the denominator of the expression of Ẽ2T .The distribution of Ẽ2T falls very fast in the small x region at small −t followed by the appearance of a rounded peak moving to higher x with increasing −t.Now referring to Fig. 2, we first note that for the γ ⊥ γ 5 related GPDs, only Ẽ′ 2T is zero.H ′ 2T behaves like Ẽ2T .

E ′
2T possesses a dipole structure in the x direction and, accordingly, has a peak in the small −t region at around x = 0.2 for both the u and the d quark as seen in Fig.

Similar structures appear in the H′
2T , where a dip is observed near the origin.We tested their structures with different parameter sets, and similar behavior persists.We suspect that this is a model-dependent feature related to the BLFQ calculations.
As depicted in Fig. 3, the behaviors of H 2 closely coincide for u and d quarks.However, for E 2 , a distinct pattern emerges, exhibiting a peak near the origin for the d quark and a distinct valley at around x = 0.2 for the u quark.
For Ẽ2 and H′ 2 in Figs. 4 and 5, they exhibit comparable trends with varying magnitudes.This alignment can be readily confirmed by examining the expressions for Ẽ2 and H′ 2 in Eqs. ( 21) and (23).
Key features of the GPDs pictured above are highlighted through 2D forms in Figs. 6 and 7. From Fig. 6, where −t is fixed at 0.25 GeV 2 and 1.44 GeV 2 , one finds that the peaks are moving towards higher x as −t increases.The 2D plots with fixed x at 2.5/16.5 as shown in Fig. 7, display the smooth trend towards zero as −t increases.
As mentioned before, the transverse symmetry between γ 1 and γ 2 has been verified in our calculations to within numerical precision.We choose −t = 1.44 GeV 2 for Ẽ2T , H ′ 2T , E ′ 2T and H′ 2T in Fig. 8 to illustrate this conclusion.Clearly, γ 1 related distributions (solid lines) and γ 2 related distributions (dashed lines) coincide with each other very well.

B. Sum Rules
Sum rules for the twist-3 GPDs G i and Gi have been derived in Ref. [74].For x 0 -moments, where i = 1, 2, 3, 4. For x 1 -moments, 1 where F 1 (−t) and F 2 (−t) are the Dirac and the Pauli form factors, G P (−t) and G A (−t) are the pseudo-scalar and the axial-vector form factors, G M (−t) and G E (−t) are the magnetic and the electric form factors respectively.From definitions, we have Note that opposite sign of the right hand side (r.h.s) of Eq. (47) after taking forward limit is exactly the socalled kinetic OAM defined in Ref. [6].This implies that twist-3 GPDs also have physical relevance that is not negligible.44)-( 53) in terms of GPDs in Eqs. ( 18) and ( 19).This work mainly focuses on the zero-skewness limit, ξ = 0, so some of the sum rules above, which contain 1/ξ or ∂/∂ξ will not be discussed here.The twist-2 GPDs used here have been calculated with the same parameter set and truncation [44,53].Note that one has to consider nonzero skewness, ξ ̸ = 0 to compute the GPD Ẽ.Using the properties of GPDs under discrete symmetries (see Appendix B), some of the GPDs are already zero at both the theoretical and the numerical level, making the sum rules (44) (i = 1, 3), (45) (i = 3) and ( 52) automatically satisfied.The differences between the left hand side (l.h.s) and the r.h.s of Eqs.(45) (i = 2, 4) ( 51) and ( 53) are not small enough to be treated as the numerical error.However, they all exhibit a common trend that they decrease as −t increases.We suspect that the deviation is due to the fact that we have only retained the leading Fock-sector and will discuss those results in more details later when we extend to the higher Fock-sectors to include contributions arising from a dynamical gluon.

C. PDF Limit
The PDF was originally introduced by Feynman to describe the deep inelastic lepton scattering process, providing information on the hadron structure in the longitudinal direction.The PDFs can be interpreted as parton densities corresponding to a specified longitudinal momentum fraction x.The PDFs are also defined by a quark-quark density matrix [75], By taking different Γ structures, we could get the PDFs for all twist.This work only focuses on twist-3 PDFs, which are given by where the indice i in Eq. ( 63) can only be 1 or 2.
It can be immediately found that some pairs of the GPDs and the PDFs are connected to each other by taking the forward limit Since ∆ appears in the denominator of the expressions of H ′ 2T , we try to perform the limit numerically to obtain the twist-3 PDFs.Fortunately, we find that our results are numerically converging when −t is smaller than 10 −7 GeV 2 .So we choose −t = 10 −20 GeV 2 as the numerical point representing the forward limit, and the results are shown in Fig. 9.
Our results satisfy the sum rule [75] xe where m q is the quark mass and N q is the number of valence quarks of flavor q.More interesting things happen to g T (x).At twist-2 level, the similar gamma structure parameterised as g 1 (x) is the so-called helicity PDF that measures the quark helicity distribution in a longitudinally polarized nucleon.At the twist-3 level, Eqs. ( 66) and (32) show that g T (x) measures the transverse spin distribution for quarks.Here there is also a well known sum rule called the Burkhardt-Cottingham sum rule [76], where g 2 (x) = g T (x) − g 1 (x).This sum rule means that the contributions of quark spin to the proton spin with different polarizations are the same.But this sum rule is derived from three-dimensional rotational symmetry, which is broken by the truncation of Fock sector and is therefore not satisfied in this work.We hope that with the inclusion of higher Fock sectors in the future this will be improved.As for h L (x), it contributes to a polarized Drell-Yan process.Both g T (x) and h L (x) are sensitive to quark-gluon interaction [75].
Recently, the point-by-point extraction of the twist-3 scalar PDF e(x) through the analysis of both CLAS and CLAS12 data for dihadron production in semi-inclusive DIS off of a proton target has been reported in Ref. [77].We also notice that there are many model calculations of e(x) for which we provide a comparison below.The proton flavor combination quantity e v is defined by [77][78][79]   In Fig. 10, we present the e v for both theoretical model calculations and experimental extractions.We utilize the Higher Order Perturbative Parton Evolution toolkit (HOPPET) to numerically solve the NNLO DGLAP equations [80], and evolve BLFQ results from the initial scale µ 2 0 = 0.195 ± 0.02 GeV 2 [44] to the scale at Q 2 = 1.0 GeV 2 and Q 2 = 1.5 GeV 2 .The error bands in our evolved distributions are due to the 10% uncertain-ties in the initial scale.We find that our results show good agreement with other model calculations and the experimental data in most regions.In the small x region the experimental data rise sharply, whereas our results rise more slowly.It should also be noted that there is large uncertainty in the experimental data in the low x region.Despite that, our preliminary results are generally consistent with most of the model calculations and the experimental data.

V. CONCLUSION
In this paper, we present all the twist-3 GPDs in the zero-skewness limit of the proton within the theoretical framework of Basis Light-front Quantization (BLFQ).The LFWFs are obtained by diagonalizing the proton light-front Hamiltonian using BLFQ.Then the LFWFs are used to calculate twist-3 GPDs via the overlap representation.We have calculated all well-defined twist-3 GPDs, both chiral-odd and chiral-even, with a truncation to leading Fock sector and numerical cut offs.We also calculate twist-3 PDFs and evolve the spin-independent PDF e(x) to a higher scale that we then compare with other theoretical model calculations and with experimental data.We find there is general agreement among the selected models and rough agreement between models and experiment considering the large experimental uncertainties.
In the future, we shall calculate twist-3 GPDs with higher Fock sectors, and study the relation between twist-3 GPDs and orbital angular momentum distribution.The study of nonzero-skewness is also of interest because it can be connected to DVCS twist-3 cross section [34] and can be measured through the experiments at the EIC and EicC [81].Results at twist-3 and non-zero skewness will deepen our understanding of the proton structure, and further, help to solve the proton spin puzzle.
where n, m are the radial and angular quantum numbers, respectively, ρ ≡ |p|/b is a dimensionless argument, b is the basis scale parameter which has mass dimension, L |m| n

FIG. 1 :
FIG. 1: Twist-3 quark GPDs in the proton associated with Γ = γ ⊥ ; (a) and (b) are for the down and up quarks on the quark level, respectively.The flavor level distributions are given by X u f lavor = 2X u quark and X d f lavor = X d quark , where X stands for all the GPDs.The GPDs are evaluated with N max = 10 and K = 16.5.

1 − 1 dxxFIG. 2 :
FIG. 2: Twist-3 quark GPDs in the proton associated with Γ = γ ⊥ γ 5 ; {(a) (c) (e)} and {(b) (d) (f)} are for the down and up quarks on the quark level, respectively.The flavor level distributions are given by X u f lavor = 2X u quark and X d f lavor = X d quark , where X stands for all the GPDs.The GPDs are evaluated with N max = 10 and K = 16.5.

FIG. 3 :
FIG. 3: Twist-3 quark GPDs in the proton associated with Γ = 1; {(a) (c)} and {(b) (d)} are for the down and up quarks on the quark level, respectively.The flavor level distributions are given by X u f lavor = 2X u quark and X d f lavor = X d quark , where X stands for all the GPDs.The GPDs are evaluated with N max = 10 and K = 16.5.

FIG. 4 :
FIG. 4: Twist-3 quark GPDs in the proton associated with Γ = γ 5 ; (a) and (b) are for the down and up quarks on the quark level, respectively.The flavor level distributions are given by X u f lavor = 2X u quark and X d f lavor = X d quark , where X stands for all the GPDs.The GPDs are evaluated with N max = 10 and K = 16.5.

FIG. 5 :
FIG. 5: Twist-3 quark GPDs in the proton associated with Γ = iσ +− γ 5 ; (a) and (b) are for the down and up quarks on the quark level, respectively.The flavor level distributions are given by X u f lavor = 2X u quark and X d f lavor = X d quark , where X stands for all the GPDs.The GPDs are evaluated with N max = 10 and K = 16.5.

2 FIG. 6 :
FIG.6:The upper row and the lower row are the twist-3 GPDs with fixed −t = 0.25 GeV 2 and fixed −t = 1.44 GeV 2 on the quark level respectively.The flavor level distributions are given by X u f lavor = 2X u quark and X d f lavor = X d quark .We use the same color to represent the same GPD for each figure, and the solid lines with fill markers and the dashed lines with open markers to represent the down and up quarks respectively.

5 FIG. 7 :
FIG. 7: The left figure and the right figure are the twist-3 GPDs with fixed x = 2.5/16.5 on the quark level.The flavor level distributions are given by X u f lavor = 2X u quark and X d f lavor = X d quark .We use one color to represent the same GPD for each figure, and the solid lines with fill markers and the dashed lines with open markers to represent the down and up quarks respectively.

FIG. 9 :|FIG. 10 :
FIG. 9: Twist-3 PDFs in the proton on the flavor level with different colour curves; the fill markers and the open markers are for the down and up quarks, respectively.The PDFs are evaluated with N max = 10 and K = 16.5.

gic
Priority Research Program of the Chinese Academy of Sciences, Grant No. XDB34000000.J. P. V. is supported by the US Department of Energy under Grant Nos.DE-SC0023692 and DE-SC0023707.A portion of the computational resources were also provided by Gansu Computing Center.This research is supported by Gansu International Collaboration and Talents Recruitment Base of Particle Physics (2023-2027) and the International Partnership Program of Chinese Academy of Sciences, Grant No.016GJHZ2022103FN.