Nucleon Helicity and Transversity Parton Distributions from Lattice QCD

We present the first lattice-QCD calculation of the isovector polarized parton distribution functions (both helicity and transversity) using the large-momentum effective field theory (LaMET) approach for direct Bjorken-$x$ dependence. We first review the detailed steps of the procedure in the unpolarized case, then generalize to the helicity and transversity cases. We also derive a new mass-correction formulation for all three cases. We then compare the effects of each finite-momentum correction using lattice data calculated at $M_\pi\approx 310$ MeV. Finally, we discuss the implications of these results for the poorly known antiquark structure and predict the sea-flavor asymmetry in the transversely polarized nucleon.

present the lattice results in Sec. V, and discuss the implications of these results in Sec. VI, focusing on the less-known antiquark distribution. The details of the finite-momentum corrections are given in the Appendices.

II. REVIEW OF THE LAMET APPROACH
The original definition of the parton distribution function (PDF) of a hadron with momentum P = (P 0 , 0, 0, P z ) depends on the correlator of the quark bilinear operator defined on a lightcone: h(ξλ · P ) ≡ 1 2λ · P P ψ (0)λ · γΓ (0, ξλ) ψ(ξλ) P , where ψ is the quark field operator, λ = (1, 0, 0, −1)/ √ 2 is a lightlike vector with λ 2 = 0, γ is the Dirac matrix and the gauge link is Γ (ζλ, ηλ) ≡ exp ig ζ η dρ λ · A(ρλ) (2) with g the strong coupling constant and A the gauge field. h (0) is normalized to the total quark number of the hadron. The physical PDF q(x) of the hadron is the Fourier transform of h: where µ is the renormalization scale. This definition is invariant under a boost along the z-direction. In particular, it is valid in the rest frame where P z = 0. Under the operator product expansion, with D the covariant derivative and higher-twist terms neglected. The operator can be written as where {...} indicates symmetrization of the enclosed indices. The tensor λ µ1 · · · λ µn is symmetric and hence only the symmetric part of O will contribute to the PDF. Furthermore, since λ 2 = 0, the tensor λ µ1 · · · λ µn is automatically traceless. For example, for n = 2 we can write λ µ1 λ µ2 = λ µ1 λ µ2 − g µ1µ2 λ 2 /4 + g µ1µ2 λ 2 /4, where the first term is symmetric and traceless while the second term is the trace term. It is clear that the trace part vanishes when λ 2 = 0. Therefore, only the symmetric and traceless part of O {µ1···µn} contributes in Eq. 5. This symmetric and traceless operator is a twist-2 operator whose matrix element is related to moment of the parton distribution a n = 1 −1 dx x n−1 q(x) and q(−x) = −q(x) with P O {µ1···µn} − traces P = 2a n (P µ1 · · · P µn − traces) .
One can check easily that with Eqs. 5-7, Eq. 3 is indeed satisfied. The definition in Eq. 3 cannot be used to compute the PDF in Euclidean space. Spacetime points on a lightcone in Minkowski space satisfy the equation t 2 − r 2 = 0. This equation becomes −t 2 E − r 2 = 0 in Euclidean space, which is only satisfied by the point at the origin. Therefore, the quark bilinear operator defined on a lightcone in Minkowski space becomes a local operator in Euclidean space, which is not desirable; Eq. 3 yields q(x) ∝ δ(x) in this case.
In principle, one can use the twist-2 operators in Eq. 7 to recover the PDF, provided all the moments of the PDF can be computed. However, in practice one can only obtain the first three moments because of the difficulty to reliably subtract the power divergence arising from mixing of higher moments with lower ones.
The first step (a) is to start from the computation of Eq. 3 with λ = (0, 0, 0, −1) and a nonzero but finite P z . Note that now the quark bilinear is equal-time, and this quantity, referred to as a "quasi-distribution"q(x, Λ, P z ), can be computed on a Euclidean lattice with an ultraviolet (UV) cutoff Λ. With λ 2 = 0, O {µ1···µn} in Eq. 5 is symmetric but not traceless. However, its structure only differs from a twist-2 operator by trace terms. Therefore, by Eq. 7, its matrix element P O {µ1···µn} P is still related to the moment of the PDF a n plus trace terms. As discussed in Ref. [14], the quark-level trace contribution on the left-hand side of Eq. 7 is a twist-4 effect and is an O(Λ 2 QCD /P 2 z ) correction, while the trace contribution on the right-hand side is an O(M 2 /P 2 z ) correction with M the nucleon mass. There is also an O(α s ≡ g 2 /4π) quantum correction to the operator O {µ1···µn} which could depend on P z as well. Taking into account these P z -dependent corrections, one can take P z → ∞ and go from step (a) to step (b). From step (b) to (c), nothing needs to be done, since, as explained above, they are the same system viewed in different Lorentz frames.
Going from step (c) to step (d) is nontrivial, but one can instead go from (b), which is identical to (c), to (d) and use the boost invariance of λ 2 = 0 to bring (d) to the P z → ∞ frame. Now both (b) and (d) have the same hadronic state with P z → ∞ but with different quark bilinear operators: the one in (b) with λ 2 < 0 while the one in (d) with λ 2 = 0. From the discussion above, this difference yields O(Λ 2 QCD /P 2 z ) and O(M 2 /P 2 z ) corrections which vanish as P z → ∞. The remaining correction is the O(α s ) Wilson coefficient renormalization in the operator product expansion of the quark bilinear. It is governed by short-distance physics and is independent of the hadronic state.
The lightlike condition λ 2 = 0 in (d) implies that the quark bilinear operator is boost invariant, and hence, the Wilson coefficient does not depend on P z , but this is not the case for (b). This is because in the former case P z is taken above the UV cutoff and is no longer dynamical, while in the latter case P z is below the UV cutoff and is still dynamical. The former can be considered as an effective field theory of the latter which is analogous to the relation between the heavy-quark effective theory (HQET) and its full theory [15]. In HQET, the heavy-quark mass is taken above the UV cutoff and is no longer a dynamical quantity while in the full theory the heavy-quark mass is below the UV cutoff and is dynamical. The difference between the two theories are compensated by higher-dimensional operators or counterterms in the effective field theory which encode the short-distance physics of the full theory that is integrated out in the effective theory. This "matching" procedure can be implemented order by order in powers of α s in perturbation theory.
Summarizing the above discussion, the quasi-distributionq(x, Λ, P z ), which can be computed in Euclidean space with λ = (0, 0, 0, −1) and nucleon momentum P z , can be related to the P z -independent physical distribution q(y, µ) with λ = (1, 0, 0, −1)/ √ 2 through [14] q(x, Λ, where µ is the renormalization scale, usually in the MS scheme, Λ will be set by the lattice spacing, and Z is the kernel from the matching. Here, we have concentrated on the flavor nonsinglet case such that the mixing with the gluon PDF is not needed.

III. FINITE-Pz CORRECTIONS FOR UNPOLARIZED PDFS
In this section, we detail the procedure to implement the P z corrections needed to extract the physical q(x, µ) from the quasi-distributionq(x, Λ, P z ) computed from the lattice. We first explain each of the corrections and then summarize the procedure at the end of the section.

A. One-Loop Matching
In the limit P z → ∞, the matching becomes the most important P z correction. The factor Z ξ = x y , µ Pz , Λ Pz has been computed up to one loop in Ref. [16] using a momentum-cutoff regulator instead of a lattice regulator. Therefore, this Z factor is accurate up to the leading logarithm but not for the numerical constant. To determine this constant, a lattice perturbation theory calculation using the same lattice action is required. At tree level, the Z factor is just a delta function, such that Since the difference between q(x) and q(x) starts at the loop level, we can rewrite the above equation as with an error of O α 2 s [28]. Z(ξ) is a singular function of ξ with terms like a/ (1 − ξ) 2 and (b ln |1 − ξ| + c) /(1 − ξ) [16]. The value of |ξ| is not bounded by unity when P z is finite. This is because Z describes the evolution of partons, which includes the gluon emission process where a mother quark splits into a daughter quark and a gluon. When P z is finite, the gluon can travel in the opposite direction from the mother quark, which makes the momentum fraction of the daughter quark bigger than that of the mother quark; hence, ξ can be bigger than one. To show that these singular terms are not harmful, we use the fact that Z has the structure with the first term coming from gluon emission and the second term from the quark self-energy diagram. C = ∞ −∞ dξ Z (1) (ξ ) such that dξ Z(ξ) = 0 and particle-number conservation is satisfied. Using this, Eq. 13 becomes Then, when ξ = x/y is close to one, the double poles 1/ (1 − ξ) 2 in Z (1) (ξ) cancel out. The single pole (b ln |1 − ξ| + c) /(1 − ξ) is odd in (y − x), which is not an endpoint singularity in y because ξ is not bounded by unity. Therefore, the integral and q(x) are finite.
Another observation is that the double-pole structure in Z (1) makes q(x) behave like α s /x 2 at large x. This can be seen by approximating q(y) in the integral by q(y) with y bounded by unity. So even though q(x) is finite, its higher moments are not. This is consistent with the fact that moments ofq(x), q(x) and Z(x) satisfy where Note that for f = q this is identical to the usual moments since q(x) vanishes outside [−1, 1]. Eq. 16 can be obtained by taking moments of Eq. 10, and the divergence of higher moments of Z implies the divergences of higher moments ofq.
In this subsection, we show how the M 2 /P 2 z correction to all orders (denoted as M 2n /P 2n z ) can be computed exactly. We will make use of the property λ µ1 · · · λ µn P (µ1 · · · P µn) = λ (µ1 · · · λ µn) P µ1 · · · P µn , where (...) means the indices enclosed are symmetric and traceless. A useful identity is where i max = n−Mod[n,2] 2 and B n,0 = 1. The B coefficients can be determined by the tracelessness of λ (µ1 · · · λ µn) which implies g µ1µ2 P µ3 · · · P µn λ (µ1 · · · λ µn) = 0, The left-hand side of this equation is a polynomial of powers of λ 2 with each term involving at most two B coefficients. Then, the identity of Eq. 21 yields the following recurrence relation: To implement the M 2n /P 2n z correction, we first compute the ratio of the moments where C is the binomial function and c = −λ 2 M 2 /4 (λ · P ) 2 = M 2 /4P 2 z with λ µ = (0, 0, 0, −1) and λ · P = P z . As shown in Appendix B, the above factors can be converted to the following relation between PDFs where f ± = √ 1 + 4c ± 1. Unlike the mass-correction expression of Ref. [25], particle number is conserved in this expression.
This correction comes from the trace part on the left-hand side of Eq. 7, which is a twist-4 effect and can be implemented by adding aq twist-4 contribution toq, such that q(x, Λ, P z ) →q(x, Λ, P z ) +q twist-4 (x, Λ, P z ). (25) As derived in Appendix C,q where Γ 0 is the incomplete Gamma function and Instead of computing these corrections directly on the lattice, we only parametrize and fit them as a 1/P 2 z correction after we have removed other leading-P z corrections.

D. Summary of the Finite Pz Corrections
Here we summarize the procedure needed to implement finite-P z corrections, and use the unpolarized u(x) − d(x) PDF as an example. We focus on the flavor-nonsinglet PDF such that there is no mixing with the gluon PDF. The generalization to the polarized case will be shown in the next section.
We start with the computation of the equal-time correlator on a Euclidean lattice: which is the lattice version of Eq. 1 with λ µ = (0, 0, 0, −1) and U µ is a discrete gauge link in the µ direction. h lat (0) is the total quark number of the hadronic state. The quasi-PDF is the Fourier transform of h lat : The second step is to implement the one-loop matching to convert q(x, P z , Λ) in the lattice scheme to q I (x, P z , µ) in the MS scheme: We will use the Z (1) factor derived in Ref. [16] (also listed in Appendix A for completeness), which matches the momentum cutoff scheme with the MS scheme. This Z (1) factor is accurate up to the leading logarithm but not for the numerical constant. One should replace the momentum cutoff calculation with a lattice perturbation theory calculation to get the correct numerical constant in the future. The third step is to remove the O(M 2n /P 2n z ) correction from q I (x, P z , µ): This gives rise to q II (x, P z , µ) whose remaining O Λ 2 QCD /P 2 z correction can be removed by adding the twist-4 contributionq twist-4 defined in Eqs. [26][27] q(x, µ) = q II (x, P z , µ) +q twist-4 (x, P z , µ). (32) q twist-4 (x, P z , Λ) can be computed on the lattice and in principle another matching to the MS scheme is required to obtainq twist-4 (x, P z , µ). However, the difference is O (α sqtwist-4 ) and hence negligible. In this work, the effect of q twist-4 will be parametrized, fitted to data and extrapolated P z → ∞ to obtain the P z -independent left-hand side of Eq. 32.

IV. SPIN-POLARIZED PDFS
In this section, the finite-P z corrections for the longitudinally polarized PDF (helicity) and the transversely polarized PDF (transversity) is documented.

A. Helicity
The lattice definition of the helicity distribution is with where M S z = P 2 z + M 2 . The one-loop matching is where ∆Z (1) from the vertex correction is given in Appendix A, while the factor Z (1) from wavefunction renormalization of Eq. 30 is the same as in the unpolarized case. The symmetric operator in Eq. 34 is whose symmetric traceless version is a twist-2 operator with matrix element where ∆a n = dx x n−1 ∆q(x) and S is the polarization vector with S 2 = −1. Using Eq. 19, we havē wherec = −4c and we have used S · P = 0. The O(M 2n /P 2n z ) correction can be removed by (see Appendix B for the detailed derivation) where The remaining O Λ 2 QCD /P 2 z correction can be removed by adding the twist-4 contribution ∆q twist-4 Again, in principle ∆q twist-4 (see Eq. 87) might be computed on the lattice directly, but here we just parametrize and fit it.

B. Transversity
The lattice definition of the transversity distribution is with The one-loop matching is where δZ (1) is given in Appendix A, while the factor Z (1) is again the same as in the unpolarized case.
For the mass correction of the transversely polarized case, we need to compute t α λ (µ1 · · · λ µn) S [α P µ1] P µ2 · · · P µn where the vector t and S are transverse to P (t · P = S · P = 0) and S [α P µ1] = (S α P µ1 − S µ1 P α ) /2. The ratio factor between moments is given by where we have used t · P = 0 in the first equality and S · P = 0 in the second one and then taken t µ = (0, 1, 0, 0). The ratio factor K n is the same as the unpolarized case in Eq. 23. We therefore have The remaining O Λ 2 QCD /P 2 z correction can be removed by adding the twist-4 contribution δq twist-4 δq(x, µ) = δq II (x, P z , µ) + δq twist-4 (x, P z , µ).
Again, in principle δq twist-4 can be computed on the lattice directly, but here we just parametrize and fit it.

V. NUMERICAL RESULTS
In this paper, we report the results of a lattice-QCD calculation using clover valence fermions on an ensemble of gauge configurations with lattice spacing a = 0.12 fm, box size L ≈ 3 fm and pion mass M π ≈ 310 MeV with N f = 2 + 1 + 1 (degenerate up/down, strange and charm) flavors of highly improved staggered quarks (HISQ) [29] generated by MILC Collaboration [30]. The gauge links are hypercubic (HYP)-smeared [31] and then clover parameters are tuned to recover the lowest pion mass of the staggered quarks in the sea 1 . HYP smearing has been shown to significantly improve the discretization effects on operators and shift their corresponding renormalizations toward their tree-level values (near 1 for quark bilinear operators). The volume of this ensemble is large enough, M π L ≈ 4.5, that there is no visible finite-volume correction in current lattice-QCD calculations of nucleon matrix elements. The results shown in this work are done using correlators calculated from 3 source locations on 449 configurations.
On the lattice, we first calculate the time-independent, nonlocal (in space, chosen to be the z direction) correlators of a nucleon with finite-P z boosth where U z is a discrete gauge link in the z direction and P = {0, 0, P z } is the momentum of the nucleon. Γ = γ z , γ z γ 5 and σ xz γ 5 for the unpolarized, helicity, and transversity distributions, respectively. The spin direction is along the z-direction for helicity and x-direction for transversity. In this work, we are only studying isovector quantities (such as the up-down flavor asymmetry). To control the systematics due to contamination by nearby excited-state quantities, we make a simultaneous fit of the nucleon matrix element correlators, using two source-sink nucleon separations, 0.96 and 1.2 fm; the detailed procedure is described in Ref. [34] for the nucleon charges. Examining the individual fits to each source-sink nucleon separation, we do not see noticeable excited-state contamination for either separation at the current statistics. Figure 2 shows the bare lattice nucleon matrix elements at the three boost momenta used here: {1, 2, 3}2π/L, which correspond to nucleon momenta of 0.43, 0.86 and 1.29 GeV, respectively. We note that in all three cases, the matrix elements vanish when the link length reaches 10-12. The signal-to-noise ratios worsen as the nucleon is increasingly boosted, so to push this method forward, future studies should investigate methods for improving nucleon momentum sources. We then take the integrals to transform the lattice matrix elements as functions of spatial link length z into the quasi-distributions as functions of parton momentum fraction x: where x = k/P z , Λ is the renormalization scale set by the lattice spacing a and C Γ = (P z /M S z ) for helicity and C Γ = 1 for unpolarized and transversity PDFs. We have sampled δk as finely as 0.002 but have not observed any dependence in downstream results on the choice of interval used here. Since the matrix elements go to zero beyond about 12, the integral does not depend sensitively on the choice of maximum z. The normalization of the long-link operators is currently estimated through zeroth moment of the quark distribution, This choice reduces the systematic uncertainty arising from the matching and other systematics such as finite-volume effects and lattice discretization. Given that the lattice renormalization constants for most observables are close to 1 on this ensemble, we will get reasonable cancellation of the remaining factors. Similar normalizations apply to  the helicity and transversity. The normalization of each distribution is then set by multiplying in the corresponding vector, axial or tensor charge, as obtained on the same lattices by Ref. [34] using standard techniques. The isovector nucleon quark, helicity and transversity quasi-distributions are shown in Fig. 3, using in the same color scheme to indicate different boosted momenta. We see that our lattice-QCD result has nonzero values for q(x), ∆q(x) and δq(x) at x ≥ 1 and that it does not vanish until x ≈ 1.5. In all three cases, the smallest momentum has the widest distribution, spreading out to large positive and negative x, beyond |x| = 1. As we discussed after Eq. 13, when P z is finite, the range of |x| is not bounded by unity. But as the boosted momentum increases, the distribution sharpens and narrows, decreasing the contribution coming from the |x| > 1 regions, just what we would expect in the lightcone distribution. This is not hard to understand (as we discussed in our earlier work [24]): in the infinite-momentum frame, no constituents of the nucleon can carry more momentum than the nucleon as a whole. However, since the momentum in our calculation is finite, the PDF does not have to vanish at x = 1. The peak location for the density and helicity distributions remains roughly the same for P z = 2 and 3, but in the case of the transversity, the peak shifts toward x = 0 for P z = 3. Note that there is a substantial difference in magnitude between P z = 2 and 3, and an even more severe difference in shape between P z = 1 and the others. We note that since x is defined as k/P z and k is arbitrary, we can make k as small as desired to obtain small-x physics. However, the small-x region corresponds to long-distance physics, which requires longer physical links to probe. This is similar to the finite-volume effect commonly seen in LQCD calculations, except the large-z links are essential to obtain a reasonable description of the physics in this region.
To improve the quasi-distribution closer to the infinite-momentum frame (IMF) proton distribution functions, we follow the recipes described in Sec. III for the one-loop and mass corrections. The effects of the one-loop (with α s set to 0.2) and the final quark distribution (one-loop first, followed by mass correction) and original quasi-distribution are shown in Fig. 4 for P z = 2 and P z = 3. We found that corrections for P z = 1 distributions are poorly behaved due to the smallness of the boosted momentum; the results are ignored here. First, we compare the quasi-(green band) and one-loop-corrected (red band) distributions. For quark density, helicity and transversity distributions, we find a significant dip caused by the one-loop correction near x = 0. The depth of this dip increases as we increase the resolution in x, dx, used in the integral; this artifact may disappear with proper one-loop renormalization in the future calculations. We also observe a clear evidence of higher values of the peak in the positive-x area and pushing outward of the peak location of the distribution. In the large-x region, the distribution is pulling back, making it rarer for quarks to carry a large fraction of momentum as one approaches the IMF, which is what we expect. For the P z = 3 distribution, the magnitude of the changes due to the one-loop correction decrease, as expected. As we expand the reach of the lattice calculation to larger values of P z , the corrections will be even smaller. The pushing outward in the large-x region may be caused by the validity of the one-loop correction requiring larger momentum. Future calculations should be designed to study this further with larger momentum and higher statistics.
We then apply the mass-correction formula to the one-loop-corrected distribution, shown as blue bands in Fig. 4 for all distributions and both P z ∈ {2, 3}. The peaks are shifted toward x = 0, the distribution sharpens, and the large-x region distribution is suppressed further, as expected. In both the quark density and transversity distributions, the mass correction also reduces the depth of the dip caused by the one-loop correction formula, and the effect of the mass correction also diminishes for the P z = 3 case. However, for the helicity, the mass-correction causes a significant unphysical spike rising near x = 0 due to the singularity in the double-integral terms. We note that the peak significantly decreases between P z = 2 and P z = 3, and this should be reduced with larger P z data in the future. The height of the peak depends on the resolution of the integral, but has very small effect on the zeroth moment. In addition, the mass-correction formulae used in this paper differ from what we used in our earlier publication, Ref. [24]. x Extrap.  This change shifts the central value of the unpolarized and longitudinally polarized up-down quark asymmetry and increases the estimated errors. However, the results remain consistent within the given errors.
To further reduce the remaining O(Λ 2 QCD /P 2 z ) correction due to higher-twist operators, we extrapolate to infinite momentum using the form a + b/P 2 z at each x point. The resulting distribution, shown in Fig. 5, has |x| > 1 region within 2 sigma of zero; thus, we recover the correct support for the physical distribution within error. Note that the smallest reliable region of x is related to the largest momentum on available on the lattice O(1/a), which is roughly the inverse of length of the lattice volume in the link direction; therefore, we expect large systematic uncertainty in the region x ∈ [−0.08, 0.08]. In the case of quark density, there are also indications of momentum convergence within 2 sigma from P z = 2 and 3 data. In addition, the final extrapolated distribution (orange band) is consistent with the largest momentum distribution. However, for the polarized distributions, even larger P z calculations are needed to improve the convergence rate and reduce the uncertainty due to extrapolation, especially for the helicity.
There are many aspects that need to be improved to get the systematics under control, as indicated at various points in the earlier sections. The operator renormalization also needs to be determined to one-loop level or better in the future calculations. We intend in this work mainly to demonstrate that one can achieve light-cone quantities with reasonable accuracy using currently available computational resources, and it opens the door for many more lattice-QCD calculations on parton physics.

VI. DISCUSSION
In this section, we take the distributions from the previous section after all corrections, shown in Fig. 6, and discuss the physics implications. We will focus on the results new to this paper, mainly the isovector helicity and nucleon transversity distribution. We strongly believe that it is worth more lattice-QCD effort to improve our knowledge of the polarized PDFs, which still lack precision experimental data over most x regions, especially the antiquark distribution, which we emphasize in this section.
The left panel of Fig. 6 shows the helicity distribution x(∆u(x) − ∆d(x)), along with selected recent global analyses by JAM [35], CCVS09 [36], and NNPDFpol1.1 [37], whose nucleon isovector distribution uncertainties have been ignored. Also note that the plots now show the distribution multiplied by x, since this form is used in global-analysis parametrizations. We see more weight distributed in the large-x region, which could shift toward smaller x as we lower the quark masses. This is because lower quark mass increases the long-range correlations in ∆h lat (z), which in turn increases the small-x contribution in the Fourier transformation. Since the increasing small-x distribution will decrease the large-x distribution due to charge conservation, we expect this inconsistency to reduce as we go to smaller quark mass. We see there are noticeable differences between the extracted polarized PDFs depending on the experimental cuts, theory inputs, parametrization, and so on. For example, JAM excludes SIDIS data, leaving the sign of the light antiquark determined by the valence and the magnitude determined from sum rules. DSSV also relies on assumptions such as SU(3) symmetry to constrain the analysis and adds a very small symmetry-breaking term. A direct lattice study of hyperon axial couplings [38] suggested that SU(3) breaking is roughly 20% at the physical point, bigger than these assumptions. Similar assumptions also made in the NNPDFpol1.1 [37]. These assumptions are unavoidable due to the difficulties of getting constraint data from polarized experiments. Future experiments with neutral-and charged-current DIS (such as at EIC) will provide useful measurements to constrain our understanding of the antiquark helicity distribution.
Our result for antiquark helicity favor more polarized up quark than down flavor 2 , which is consistent with the current PDF analysis and model calculations, such as chiral quark soliton model (χQSM) [39]. This was first pointed out in our earlier paper [24] (using a different mass-correction formulation), which concentrated on the sea flavor asymmetry in the unpolarized distribution; it was also noted in preliminary studies in conference proceedings [22,23,40]. The sea flavor asymmetry was confirmed in the full analysis of the Run-9 data by both STAR [26] and PHENIX [27] collaborations. RHIC experiments on longitudinal single-spin asymmetry and parity-violating W production at RHIC might shed more light on the polarized sea distribution [41].
We see a moderate polarized total sea asymmetry, 1 0.08 ∆u(x) − ∆d(x) = 0.14(9), which is smaller than the previous determination [24] but still consistent within errors. The update is due to the application of the masscorrection formula of Eq. 39 instead of Eq. 38. The latter requires transforming back and forth between the PDF and moments, which introduces oscillatory artifacts. Most QCD models predict smaller polarized sea asymmetry; for example, see the recent review article by Chang and Peng, in particular Table 5 of Ref. [42]. χQSM, a large-N c model, gives rather different results by predicting a large polarized sea asymmetry: 0.31. Unfortunately, our current statistical error does not help rule out many models yet based on the total sea asymmetry. On the experimental global analysis side, the total polarized sea asymmetry estimated by DSSV09 is consistent with zero within 2 sigma, and the central value is also smaller (≈ 0.07) than the unpolarized case. Current results for STAR [26] and PHENIX [27] in the middle-x range do not clarify what the total asymmetry would be. The upcoming RHIC data from Run-13 with significantly improved statistics may shed some light on this matter. The upcoming Fermilab Drell-Yan experiments (E1027/E1039) can also provide precise experimental input on the polarized sea asymmetry magnitude.
The transversity distribution is the least known PDF among the three PDF structures studied; there is much less information available due to the difficulties in experiments. There have been a few attempts to extract the transversity distribution, but they suffer from fundamental defects. Ref. [43] makes various assumptions such as the evolution form and that there is no antiquark contribution. Ref. [44] uses dihadron fragmentation functions using data from HERMES and COMPASS analysis of pion-pair production in DIS off a transversely polarized target for two combinations of "valence" (q +q) helicity distribution. They have a proper Q 2 evolution but still rely on assumptions such as the Soffer inequality. Kang et al. [45] has improved evolutions implemented in their analysis, but they also make the assumption that the sea asymmetry is zero. The distribution for the positive x goes quickly to zero, likely due to lack of data.
Our transversity result is shown in the the right panel of Fig. 6, along with an estimate from a QCD model, χQSM [39] and the latest transversity fit from Ref. [45] 3 . Surprisingly, our result is rather similar to χQSM within 90% confidence, but with slower descent to zero in the x ≈ 1 region, similar to the quark distribution. This can be, again, due to the heavier pion mass used in the calculation, as well as the need to push for even larger momenta. In contrast, the phenomenological results from Ref. [45] fall faster as x approaches near 1.
Our result favors δd(x) > δu(x) with total sea asymmetry 0. 10(8), whose central value is still larger than most model predictions (for example, χQSM estimates 0.082 asymmetry) and in contradiction to the assumption that the antiquark is consistent with zero in some transversity extractions using experimental data [43][44][45][46]. One interesting thing to note is that the central values of the lattice determination of the tensor charge g T (that is, +1 −1 δu(x) − δd(x)) extrapolated to the continuum limit from various groups are consistently higher than the phenomenological ones who assume zero total sea asymmetry in transversity; see the summary plot Fig. 10 in Ref. [32]. This may indicate nonzero sea contribution with the same sign as our prediction here, or missing larger-x data in containing their fit. It would be interesting to see whether such a nonzero sea asymmetry remains in the future high-statistics physical quark mass ensemble; it is certainly contrary to traditional expectation. Improved phenomenological analysis with new experimental data would also help to narrow the phenomenological uncertainties and explore the discrepancy.
The cleanest measurement of the transversity would have both a polarized beam and polarized target, but given the limited setups available, once again, more data are needed. PHENIX and STAR will be able to help give more insight into this quantity. Planned experiments, such as SoLID at Jefferson Lab, can provide good transversity measurements for a wide range of positive x. The Drell-Yan experiment at FNAL (E1027+E1039) can in principle extract sea-asymmetry information in the near future to settle the size of the total transversely polarized sea.

Acknowledgments
The LQCD calculations were performed using the Chroma software suite [47]. 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, and on Hyak clusters at the University of Washington managed by UW Information Technology, using hardware awarded by NSF grant PHY-09227700. We thank MILC Collaboration for sharing the lattices used to perform this study. This work was partially supported by the U. In this Appendix, we list the one-loop matching factors used throughout this paper. As the UV cutoff in a practical calculation is finite, we use the results of Eqs. 5, 6, 21 and 24 in Ref. [16].
In the unpolarized case, the matching factors are given as follows: to the asymmetric errorband reported in the components; the error given here might be larger than if derived using the correlations of the original analysis. Also, note that the scale is set at 10 GeV 2 ; however, there is only a small difference in their central values between lower-Q 2 scale and 10 GeV 2 . where Λ(x) = Λ 2 + x 2 P 2 z . We have not taken the Λ xP z limit, because they could be the same order on the lattice.
For the helicity distribution, we have and for the transversity distribution, we have where the quark self-energy contribution is the same as Eq. 52.
We have derived the mass correction up to O(r 3 ). Although in the present work we did not implement O(r 3 ) correction due to its computational complexity, the above result will be useful for future improvements with more computational resources.
Then we can define the trace term as a twist-4 PDF that needs to be subtracted: where Γ 0 is the incomplete Gamma function satisfying 1 0 dt t e ix/t = Γ 0 (−ix) .