Hadron spin structure from lattice QCD

We present results on the spin carried by quarks and gluons in the nucleon using lattice QCD simulations with physical values of the light, strange and charm quark masses. We also discuss selective results on the x-dependence of parton distribution functions computed in lattice QCD employing the quasi-parton distribution approach and the large momentum effective theory.


Introduction
Quantum Chromodynamics (QCD) is the fundamental theory of the strong interactions that describes a plethora of complex nuclear physics phenomena including giving rise to hadron states, formed by quarks and gluons. The most fundamental hadron is the proton, which, being a stable particle, has served for many years as a laboratory for studying the complex dynamics of QCD, both theoretically and experimentally. A key question that has been posed for many years is how the spin of the proton arises from its constituents, the quarks and the gluons, collectively called partons. The successful quark model, that describes well properties of the low-lying hadrons, predicted that all the spin is carried by the three valence quarks. A major surprise came from the measurements of the European Muon Collaboration (EMC) [1,2] that determine the proton spin-dependent structure function down to momentum fraction x = 0.01. The conclusion was that only about half of the proton spin is carried by the valence quarks. This came to be known as the proton spin puzzle. It triggered a series of precise measurements by the Spin Muon Collaboration (SMC) in 1992-1996 [3] and by COMPASS [4] since 2002. Recent experiments using polarized deep inelastic lepton-nucleon scattering (DIS) processes indeed confirmed that only about 25-30% [5][6][7][8][9][10] of the nucleon spin comes from the valence quarks.
While experiments play a crucial role in the understanding of the origin of the proton spin, a theoretical approach to compute the various contributions is equally important. However, due to the nonperturbative nature of QCD, computing the parton contributions to the spin is difficult. See Ref. [11] for a recent theoretical overview on the nucleon spin. Lattice QCD provides the initio non-perturbative framework that is suitable to address the key question of how the nucleon spin and momentum are distributed among its constituents using directly the QCD Lagrangian. Tremendous progress has been made in simulating lattice QCD in recent years. State-of-the-art simulations are currently being performed with dynamical up, down and strange quarks with masses tuned to their physical values. A subset of simulations also include a dynamical charm quark with mass fixed to approximately its physical value. A particularly suitable discretization scheme for studying hadron structure is the so called twisted mass fermion (TMF) action employed by the Extended Twisted Mass Collaboration (ETMC). This progress was made possible using efficient algorithms and in particular multigrid solvers [12] that were developed for twisted mass fermions [13]. In Fig. 1 we show the ensembles generated by ETMC as well as the landscape of currently available gauge ensembles worldwide. As can be seen, ETMC and a number of other collaborations have now ensembles generated with physical values of the pion mass, referred to as physical point ensembles. We also note that currently all ensembles are generated with lattice spacing a > 0.05 fm due to the so-called critical slow down of simulations. New approaches that include machine learning are being explored for addressing this issue. Fig. 1. Left: Ensembles of TMF configurations generated with a mass degenerate doublet of up and down quarks, referred to as light quarks and denoted as N f = 2 ensembles, as well as with a strange and charm quark (N f = 2 + 1 + 1), as the mass of the light quarks is varied and as function of the lattice spacing a. The mass of the strange and charm quarks are tuned to approximately their physical values. The size of the symbol is proportional to the spatial extent of the lattice. Right: A summary of zero-temperature simulations that include clover, twisted mass, staggered and domain wall fermions by various lattice QCD collaborations, as indicated in the legend.
Although lattice QCD provides an exact formalism for solving QCD, one needs a careful study of systematic errors. In the past, the absence of simulations at physical values of the light quark mass necessitated a chiral extrapolation. For the pion sector such an extrapolation using e.g. NLO SU(2) chiral perturbation theory works well for pion masses m π < ∼ 250 MeV [14]. In the nucleon sector, chiral extrapolation is more problematic and introduces an uncontrolled systematic error [15]. Nowadays, with simulations with physical pion mass this systematic error can be eliminated. Therefore, in what follows we will focus on results obtained using simulation generated with physical pion mass, since we will focus on nucleon properties. Other systematic effects that need to be investigated are: • Discretization effects: Since the computation is done at finite lattice spacing a one needs to take the continuum limit. This requires simulations with at least three values of a, ideally at fixed quark masses and volume.
• Volume effects: A numerical evaluation is necessarily done using a finite volume. At least three volumes would be required ideally at fixed quark masses and a to take the infinite volume limit.
• Renormalization: Matrix elements computed on the lattice must be properly renormalized in order to compare with what is measured in the laboratory. In state-of-the-art lattice QCD computations, renormalization is carried out nonperturbatively. However, for gluon and flavor singlet quantities there is mixing and nonperturbative renormalization is still difficult and a subject of on-going study.
• Extraction of nucleon matrix elements: Extracting the ground state from a tower of higher excited states necessitates taking the large Euclidean time limit that gives rise to large gauge noise that makes difficult the identification of ground state properties and requires very large statistics. The development of noise reduction techniques is thus an important ongoing effort.
The work-flow of a typical lattice QCD computation for baryon structure is shown schematically in Fig. 2. After defining the theory on a discrete finite 4-dimensional Euclidean lattice, one generates

Spin structure in lattice QCD
In lattice QCD studies, one typically starts from the decomposition of the energy-momentum tensor (EMT) into traceless and trace parts [16]: The traceless part,T αν , decomposes into two scheme and scale dependent terms, given byT where µ is the energy scale and curly braces mean symmetrization over indices and subtraction of trace. These are related, respectively, to the quark and gluon total angular momentum.

Nucleon matrix elements
In order to compute the nucleon spin, we need to evaluate the nucleon matrix elements of the EMT. They can be decomposed into generalized form factors (GFFs) that depend only on the momentum transfer squared q 2 . In Minkowski space we have [17] where u N is the nucleon spinor with initial (final) momentum p(p ′ ) and spin s(s ′ ), P = (p ′ + p)/2 is the total momentum and q = p ′ − p the momentum transfer. A q,g 20 (q 2 ), B q,g 20 (q 2 ) and C q,g 20 (q 2 ) are the three GFFs. In the forward limit, A q,g 20 (0) gives the quark and gluon average momentum fraction x q,g . Summing over all quark and gluon contributions gives the momentum sum x q + x g = 1. Furthermore, the total spin carried by a quark is given by [18]. The lattice operatorT µν q,g must be renormalized before we can connect it to physical matrix elements. For the flavor nonsinglet combination, this operator is multiplicatively renormalized. The standard approach to compute the nonsinglet renormalization function Z q nonperturbatively is to use an analysis within a regularization indepedent (RI) scheme. The flavor singlet operator, however, mixes with the gluon operator and one needs to compute a 2 × 2 matrix element. The mixing coefficients (Z qg , Z gq ) are thus needed to extract the momentum fraction and spin carried by quarks and gluons. The 2 × 2 mixing matrix needed for the proper renormalization procedure of the momentum fractions is given by In the above equations, x q + is understood to be the flavor singlet combination that sums the up, down, strange and charm quark contributions. The subscript B represents the bare matrix elements. We note that a complete calculation of the 2 × 2 mixing matrix would require the solution of a system of four coupled renormalization conditions that involve vertex functions of both gluon and quark EMT operators. One can impose decoupled renormalization conditions for the nonperturbative calculation of the diagonal elements, namely Z gg and Z qq , since the mixing coefficients Z gq and Z qg are small in one-loop perturbation theory [19,20]. The computation of Z qq proceeds analogously to the non-singlet case. The computation of the renormalization function Z gg of the gluon operator is more challenging due to the high noise-to-signal ratio appearing in the calculation of gluonic quantities. To this end, some equivalent renormalization prescriptions have been proposed to reduce the statistical uncertainties in the extraction of Z gg . Details on the different renormalization approaches of the gluonic operator are given in Refs. [20][21][22]. The off-diagonal elements so far are only computed in perturbation theory.

Extraction of nucleon matrix elements
In order to compute the nucleon matrix elements of the previous section, one needs to evaluate two-and three-point functions in Euclidean space. The nucleon two-point function is given by where J N (t, x) is the nucleon interpolator, x 0 is the initial lattice site at which states with the quantum numbers of the nucleon are created, referred to as source position and x s is the site where they are annihilated, referred to as sink. The three-point function is given by The Euclidean momentum transfer squared is given Q 2 = −(p ′ − p) 2 and Γ ρ is the projector acting on the spin indices. We consider Γ 0 = 1 2 (1 + γ 0 ) and Γ k = iΓ 0 γ 5 γ k taking the non-relativistic representation for γ µ .
In order to extract the desired nucleon matrix element, we construct appropriate combinations of three-to two-point functions, which in the large Euclidean time limit, cancel the time dependence arising from the time propagation and the overlap terms between the interpolating field and the nucleon state. An optimal choice that benefits from correlations is the ratio [23,24].
The sink and insertion time separations t s and t ins are taken relative to the source time t 0 . Taking the limits (t s − t ins ) ≫ a and t ins ≫ a in Eq.(8), filters out lowest nucleon state. When this happens, the ratio becomes independent of time and yields the desired nucleon matrix element. In practice, (t s − t ins ) and t ins cannot be taken arbitrarily large, since the signal-to-noise ratio decays exponentially with the sink-source time separation. Therefore, one needs to take (t s − t ins ) and t ins large enough so that the nucleon state dominates in the ratio. To identify when this happens is a delicate process and we use the following three methods: Plateau method: If contributions from excited states are sufficiently suppressed, the ratio of Eq. (8) can be written as where ∆E is the energy gap between the nucleon state and the first excited state. The first timeindependent term is the matrix element we interested in. Thus, we seek to identify nucleon state dominance by looking for a range of values of t ins for which the ratio of Eq. (8) is time-independent (plateau region). We fit the ratio to a constant within the plateau region and seek to see convergence in the extracted fit values as we increase t s . If such a convergence can be demonstrated, then the desired nucleon matrix element can be extracted. Summation method: One can sum over t ins the ratio of Eq. (8) to obtain Assuming the nucleon state dominates over excited state contributions, the desired matrix element given by Π µν (Γ ρ ; p ′ , p) is extracted from the slope of a linear fit with respect to t s . As in the case of the plateau method, we probe convergence by increasing the lower value of t s used in the fit until the resulting values for the slope converge. While both plateau and summation methods assume that the ground state dominates, the exponential suppression of excited states in the summation is faster and approximately corresponds to using twice the sink-source time separation t s in the plateau method. Two-state fit method: One can explicitly include the contributions from the first excited state and fit. The expressions and details on the analysis are given in Ref. [20].

Results on the spin decomposition of the nucleon
After determining the nucleon matrix element from the ratio of Eq. (8 we can extract nucleon charges, moments and GFFs. In Fig. 3, we show our results for the proton average momentum fraction for the up, down, strange and charm quarks, for the gluons as well as their sum. The up quark makes the largest quark contribution of about 35% and it is twice as big as that of the down quark. The strange quark contributes significantly smaller, namely about 5% and the charm contributes about 2%. The gluon has a significant contribution of about 45%. Summing all the contributions results to q=u,d,s,c confirming the expected momentum sum. Fig. 3, showing connected and disconnected contributions, demonstrates that disconnected contributions are crucial and if excluded would result to a significant underestimation of the momentum sum. Fig. 3. The decomposition of the proton average momentum fraction x P (left) and spin J P . We show the contribution of the up (red bar), down (green bar), strange (blue bar) and charm (orange bar) quarks and their sum (purple bar), the gluon (cyan bar) and the total sum (grey bar). Note that what is shown is the contribution of both the quarks and antiquarks (q + = q +q). Whenever two overlapping bars appear the darker bar denotes the purely connected contribution while the light one is the total contribution, which includes disconnected taking into account also the mixing. The error bars on the only connected part are omitted while for the total are shown explicitly on the bars. The percentages written in the figure are for the total contribution. The dashed horizontal line is the momentum and spin sums. Results are given in MS scheme at 2 GeV.
The individual contributions to the proton spin are presented in Fig. 3. The major contribution comes from the up quark amounting to about 40% of the proton spin. The down, strange and charm quarks have relatively smaller contributions. All quark flavors together constitute to about 60% of the proton spin. The gluon contribution is as significant as that of the up quark, providing the missing piece to obtain J P = 94.6(14.2)(2.8)% of the proton spin, confirming indeed the spin sum. The  (28), which is indeed compatible with zero.
Since the quark contribution to the proton spin is computed, it is interesting to see how much comes from the intrinsic quark spin. In Fig. 4 we show our results for the intrinsic quark spin 1 2 A is the axial charge for each quark, the values of which are taken from Ref. [25]. The up quark has a large contribution of about 85% of the proton intrinsic spin. The down quark contributes about half compared to the up and with opposite sign. The strange and charm quarks also contribute negatively with the latter being about five times smaller than the former giving a 1% contribution. The total 1 2 ∆Σ q + is in agreement with the upper bound from COMPASS [26] Fig. 4. Results for the intrinsic quark spin 1 2 ∆Σ (left) and orbital angular momentum L (right). We given the contributions to the proton spin decomposed into the up (red bar), down (green bar), strange (blue bar), and charm (orange bar) quark parts. The total contribution of the four flavors is also shown (grey bar) [25]. The dashed horizontal line is the observed proton spin and the percentages are given relative to it. Results in MS scheme at 2 GeV.
Having both the quark total angular momentum and the quark intrinsic spin allows us to extract the orbital angular momentum L q + P . Our results are shown in Fig. 4. The orbital angular momentum of the up quark is negative reducing the total angular momentum contribution of the up quark to the proton spin. The contribution of the down quark to the orbital angular momentum is positive almost canceling the negative intrinsic spin contribution resulting to a relatively small positive contribution to the spin of the proton. For a direct calculation of L P using TMDs see Ref. [27].

Direct computation of x-dependencne of parton distributions
While for many years only GFFs were accessible in lattice QCD, a new approach was proposed that enables us to compute the x-dependence of parton distribution functions (PDFs) within lattice QCD taking advantage of the so-called large momentum effective theory (LaMET) [28]. The main idea is to compute spatial correlators in lattice QCD e.g. along z-axis and boost the nucleon state in the same direction to large momentum. After renormalization, a matching kernel computed perturbatively is used to recover the light-cone correlation matrix element. Here we will present results using the quasi-PDF approach for which there have been many studies. For recent reviews see Refs. [29][30][31][32].
The starting point is the renormalized space-like matrix elements for boosted nucleon states where W is a Wilson line. The quasi-PDFF Γ is then matched using the light-cone PDF F Γ (x, µ). We extract the unpolarized, helicity and transversity PDFs when we take for Γ = γ 0 , γ i γ 5 and σ i, j , i j, respectively. Results using an N f = 2 + 1 + 1 twisted mass fermion ensemble on the three isovector PDFs are shown in Fig. 5. For details we refer to Refs. [33,34]. 38 GeV (blue band) and compared to phenomenological fits obtained using SIDIS data (gray) and SIDIS data constrained using the nucleon tensor charge computed in lattice QCD [35].  [36]. We compare with results from phenomenological analyses by JAM17 (green [37] and NNPDF POL 1.1 (red) [38].
To determine the flavor decomposition of PDFs one needs to compute, besides the isovector, the isoscalar PDFs. The latter as well as the strange PDFs involve disconnected contributions. These are computed by extending the formalism for local operators to extended operators. The quark disconnected loop is given by Tr D −1 q (x ins ; x ins + z)γ 3 γ 5 W(x ins , x ins + z) , where we took the Wilson line along the z-direction. The helicity up, down and strange PDFs are shown in Fig. 6 [36,39]. The quasi-PDF approach can be extended to study generalized parton distributions (GDPs), where momentum transfer is involved. One computes the nucleon matrix element

Conclusions
Moments of PFDs can be extracted precisely making these quantities part of the precision era of lattice QCD. From these moments we can extract a lot of interesting physics and also reconstruct the PDFs, as was demonstrated for the pion [41]. New theoretical approaches (quasi-distributions, pseudo-distributions, current-current correlates, etc) now allow the direct computation of PDFs within lattice QCD yielding very promising results including calculations using simulations with physical pion mass. The calculation of sea quark contributions to PDFs is shown to also be feasible providing valuable input e.g. for the determination of the strange helicity PDF. Exploratory studies of isovector GPDs are also under way, see e.g. recent studies by ETMC [40,42,43]. The way forward includes studying the continuum limit and the volume dependence and developing noise reduction techniques to reach larger boosts. Extending these approaches to twist-3 PDFs [44], TMDs, and other hadrons is also in the pipeline.