Lattice QCD constraints on the parton distribution functions of ${}^3\text{He}$

The fraction of the longitudinal momentum of ${}^3\text{He}$ that is carried by the isovector combination of $u$ and $d$ quarks is determined using lattice QCD for the first time. The ratio of this combination to that in the constituent nucleons is found to be consistent with unity at the few-percent level from calculations with quark masses corresponding to $m_\pi\sim 800$ MeV, extrapolated to the physical quark masses. This constraint is consistent with, and significantly more precise than, determinations from global nuclear parton distribution function fits. Including the lattice QCD determination of the momentum fraction in the nNNPDF global fitting framework results in the uncertainty on the isovector momentum fraction ratio being reduced by a factor of 2.5, and thereby enables a more precise extraction of the $u$ and $d$ parton distributions in ${}^3\text{He}$.

The fraction of the longitudinal momentum of 3 He that is carried by the isovector combination of u and d quarks is determined using lattice QCD for the first time. The ratio of this combination to that in the constituent nucleons is found to be consistent with unity at the few-percent level from calculations with quark masses corresponding to mπ ∼ 800 MeV, extrapolated to the physical quark masses. This constraint is consistent with, and significantly more precise than, determinations from global nuclear parton distribution function fits. Including the lattice QCD determination of the momentum fraction in the nNNPDF global fitting framework results in the uncertainty on the isovector momentum fraction ratio being reduced by a factor of 2.5, and thereby enables a more precise extraction of the u and d parton distributions in 3 He.
A central pillar of our understanding of the internal structure of strongly interacting hadronic and nuclear systems is knowledge of their partonic structure as accessed in deep-inelastic scattering (DIS) experiments and other hard processes. Since the 1960s, such experiments have revealed the longitudinal momentum distributions of quarks and gluons in a fast moving proton, known collectively as parton distribution functions (PDFs). The simplest PDFs, q(x, µ) (and g(x, µ)), describe the probability of a quark of flavor q (or gluon g) carrying a fraction x of the longitudinal momentum of the struck proton at a renormalization scale µ. In 1983, the European Muon Collaboration (EMC) [1] observed that the partonic structure of nuclei differs substantially from that of the constituent protons and neutrons, a landmark in the development of nuclear physics [2][3][4][5][6]. Since the DIS processes observed in the EMC experiments were at very high energy, and the binding energy of a nucleus is small in comparison to its mass, the appearance and size of the EMC effect was surprising at the time. Interest in the EMC effect has been rekindled by recent data from SLAC and JLab [7][8][9][10][11] on EMC ratios for light nuclei. Not only have these data provided precise determinations of the EMC effect for nuclei with small atomic number A, but they have revealed a correlation between the strength of the EMC effect and so called "short range correlations" [12,13].
In addition to experimental investigations, theoretical calculations of the partonic structure of hadrons and nuclei from the Standard Model can have important impact on our understanding of the structure of matter.
For example, Standard Model calculations of nuclear partonic structure would reveal the QCD origin of the EMC effect as well as aid in the flavor-separation of proton PDFs. Parton distributions are inherently rooted in the strong interaction dynamics of QCD and cannot be determined using perturbative methods. Since the seminal works of Refs. [14,15], lattice Quantum Chromodynamics (LQCD) calculations have addressed the simplest aspects of the parton distributions of the proton, notably determining the first few Mellin moments of the unpolarized, polarized, and transversity quark distributions [16], as well as their gluonic analogues [17][18][19][20][21]. Recently, efforts have been made to extend these studies to the full x-dependence of the proton PDFs [16,[22][23][24]. More complicated extensions of partonic structure, such as generalized parton distribution functions and transversemomentum dependent parton distribution functions of the proton, have also been studied using LQCD [25][26][27][28][29][30].
In this letter, the partonic structure of light nuclei is studied in LQCD for the first time through an investigation of the isovector quark momentum fractions (the first moments of the corresponding isovector PDFs) of the proton, diproton, and 3 He. At the heavier-than-physical quark masses used in this LQCD study, percent-level nuclear effects are resolved in the momentum fraction of 3 He. After an extrapolation to the physical quark masses, these calculations provide a constraint on the isovector momentum fraction that is used as an additional input into the nNNPDF2.0 [31] global nuclear PDF analysis framework. Since the isovector combination of nuclear PDFs is poorly determined from experiment, this LQCD constraint significantly reduces the uncertainties on the 3 He PDFs and thereby improves knowledge of nuclear structure.
LQCD methodology -The existence of strong interactions between quarks and gluons necessitates the use of LQCD for calculations of the partonic structure of nuclei. The calculations presented here are performed using a single ensemble of gauge-field configurations generated with a Lüscher-Weisz gauge action [32] with N f = 3 degenerate light-quark flavors with the clover-improved Wilson fermion action [33], and quark masses tuned to produce a pion mass of m π = 806 MeV. The lattice geometry is L 3 × T = 32 3 × 48, and the lattice spacing is determined to be a ∼ 0.145 fm from Υ spectroscopy [34]. This ensemble, and two others with different spacetime volumes have previously been used to study the spectrum [34,35] and properties [36][37][38][39][40][41][42][43][44][45][46] of light nuclei up to atomic number A = 4. The multi-volume spectroscopy studies show that the pp and 3 He states that are investigated here are bound systems with infinite volume energies below threshold. Consequently, matrix elements in these states are expected to receive only exponentially small finite volume effects, O(e −κL , e −mπL ), that will be neglected in this work [47][48][49][50][51][52][53].
The Mellin moments of the unpolarized isovector quark PDFs, q 3 (x, µ), are determined from matrix elements of twist-two operators as where p is the momentum of the state h, τ 3 is a Pauli matrix in flavor space, where D µ is the gauge covariant derivative, and {. . .} indicates symmetrization and trace-subtraction of the enclosed indices. The above operators are constructed to transform irreducibly under the Lorentz group, but the hypercubic spacetime lattice used in the LQCD calculations reduces these symmetries, in general inducing mixing between operators of different Lorentz spin. In particular, the two-index operators that determine the isovector quark momentum fraction, x  representation [54] are computed, namely where γ ν is also Euclidean. With a lattice regulator, this operator is discretized as a covariant finite difference whose form is given in the Supplementary Material. For both spin-zero and spin-half systems, spin-averaged in the latter case, matrix elements in states with zero three-momentum determine the momentum fraction as The renormalized operator in the modified minimal subtraction scheme (MS) is related to the bare lattice operator in Eq. (2) as where the renormalization coefficient Z RI MOM (µ 0 , a) is defined non-perturbatively in a regularizationindependent momentum-subtraction scheme [55] at a scale µ 0 and then matched to MS through the three-loop perturbative coefficient R MS/RI MOM (µ, µ 0 ) [56,57], as detailed in the Supplementary Material. For µ = 2 GeV, R MS/RI MOM (µ, µ 0 )Z RI MOM (µ 0 , a) = 0.89 (4).
The techniques needed to compute matrix elements of this operator are simple generalizations of those used for calculations of isovector matrix elements of quark currents using the compound-propagator background-field method introduced in Ref. [42] and further detailed in Refs. [43,44,58] and the Supplementary Material. Quark propagators and T -compound propagators are computed from an average of N src = 24 source points randomly distributed on N cfg = 2290 gauge-field configurations for N B = 5 different background field strengths. These compound propagators are then used to construct baryon two-point correlation functions, where λ is the T -background field strength, χ h is an interpolating field for states with the quantum numbers of the hadron or nucleus h, and spinor indices on the interpolating operators are suppressed. Correlation functions are constructed from Gaussian-smeared source interpolating operators [59], while the sink interpolating operators are either smeared or point-like and the multibaryon contractions are performed using the techniques of Refs. [60]. This quantity contains responses to the field up to O(λ N Q ), with N Q being the number of valence quarks in the state. The linear response of this background-field two-point function, G h (t; λ)| O(λ) , is determined by the matrix element of T . This term can be extracted exactly from the computed set of fixedorder background field correlation functions with N Q field strengths [42,43].
Combining the linear response of the two-point correlation function with the zero-field correlation function, it is straightforward to show that the ratio is related to matrix elements of T through the spectral representation of each term in Eq. asymptoting as with exponentially vanishing contamination at early times that involves excited-state overlap factors and transition matrix elements. Ground-state matrix elements are extracted from R h (t), and systematic fitting uncertainties are estimated, using a procedure for sampling from all possible fit ranges and models analogous to the procedure described for twopoint correlation functions in Ref. [61]. In summary, in analyzing R h (t) to extract the momentum fractions, the full t dependence that results from the spectral decomposition of each term in Eq. (5) is fit, and combined fits to two-and three-point correlation functions are used to constrain the relevant energies, overlap factors, and matrix elements. All possible choices of fit ranges and up  (14) 1.028 (15) TABLE I. The isovector quark momentum fractions in p, pp and 3 He, calculated at mπ = 806 MeV in MS-scheme at µ = 2 GeV. The first uncertainty combines LQCD statistical and systematic uncertainties and the second uncertainty is from operator renormalization. The correlated ratios of the isovector momentum fraction in nuclei to those in the constituent nucleons, in which the renormalization constants and their uncertainties cancel, are also given.
to 4 states contributing to the spectral decompositions are considered using a model selection process described in the Supplementary Material. A weighted average over fits from all acceptable fit ranges is used to define groundstate energy results, including systematic uncertainties from fit range and model variation. Results are shown in Fig. 1 for the proton, diproton and 3 He.
Results and Discussion -The extracted values of the isovector quark momentum fractions for p, pp, 3 He at quark masses corresponding to m π = m K = 806 MeV are shown in Tab. I and displayed graphically in Fig. 2. The uncertainties are separated into those from the LQCD calculation of the bare matrix elements, and the (larger) uncertainty from the renormalization and matching to the MS scheme. The proton isovector momentum fraction is consistent with other LQCD extractions at similar values of the quark masses [62] given the different renormalization procedures and lattice spacings. The pp and 3 He momentum fractions are determined with O(5%) uncertainties and are found to be approximately consistent with those of the constituent nucleons. The ratios of the nuclear momentum fractions to that of the proton are independent of operator renormalization to O(α s ), and are determined at few-percent precision even for 3 He.
In Refs. [63][64][65], nuclear effective field theory (EFT) was used to study nuclear effects in PDF moments. In particular, it was shown that the leading source of such effects is the two-nucleon correlations that couple to the twist-two operators defining the PDF moments. In terms of the parameters defined in that work, nuclear effects in the isovector momentum fraction are encapsulated in the low energy constant (LEC) α 3,2 and nuclear factor G 3 ( 3 He); their product is bounded as α 3,2 G 3 ( 3 He) = 0.0018 (14) at µ = 2 GeV from the numerical calculations presented here (see the Supplementary Material for details). While the quark momentum fractions themselves have nonanalytic dependence on the quark masses [66][67][68], this two-body LEC is expected to be relatively insensitive to variation of the quark masses, as seen for the the analogous two-body contribution in the np → dγ [39] and pp → de + ν e [42,69] processes. This relative mass-independence assumption allows an extrapolation to the physical quark masses: a naive estimate is given by taking the central value determined at m π = 806 MeV and inflating the uncertainty by 50% to account for possible quark-mass dependence as well as the effects of the nonzero lattice spacing and finite volume (this uncertainty is estimated based on the mass dependence seen for the analogous two-body LECs in Refs. [39,42,69]). This extrapolated value can be combined with the physical value of the nucleon momentum fraction, x (7) at µ = 2 GeV from the nNNPDF2.0 analysis [31], to determine the isovector momentum fraction ratio 3 x (26) at the physical quark masses (see the Supplementary Material for more details).
It is interesting to compare the LQCD results for the momentum fractions and their ratios to phenomenology. In particular, the isovector momentum fractions determined here provide valuable information that is complementary to experimental constraints on the nuclear modification of PDFs; almost all information on the nuclear modification of partonic structure has been obtained for the ratio of isoscalar-corrected F 2 structure functions of nuclei to that of the deuteron [3,5,6]. Additional constraints are especially valuable in the context of the intriguing question as to whether there is flavordependence to the EMC effect. Such flavor dependence has been conjectured in models of QCD [70][71][72][73][74][75] and in EFT [63][64][65] and is included in recent data-driven analyses of experimental results [76,77] and provides a potential explanation of the NuTeV anomaly in sin 2 θ W [78]. sumption that the nuclear effects vary slowly with A). In this way, correlations between the 3 He and proton PDFs are accounted for. For the isovector combination, the 68% confidence interval is 3 x u−d | nNNPDF2.0 = 1.007 (63). In the nNNPDF approach, it is also straightforward to impose the LQCD constraint on the nuclear PDFs by reweighting the Monte Carlo replicas as discussed in Ref. [79]; the combined confidence region is shown in Fig. 3. The 68% confidence interval reduces to 3 x (25). Fig. 4 compares the ratio of the isovector PDF for 3 He to that of the constituent nucleons, with and without the imposition of the LQCD constraint. As can be seen from the reduced uncertainties in Figs. 3 and 4, LQCD calculations such as those presented here, as well as new experimental constraints [80,81], can significantly improve our knowledge of the flavor dependence of nuclear PDFs.
Summary -In this work, the isovector momentum fractions of the proton, diproton and 3 He systems have been determined using LQCD, complementing a previous study of the gluon momentum fraction on the same ensemble [45]. These calculations were performed at a single set of unphysical SU(3)-symmetric values for the quark masses corresponding to m π = 806 MeV, and in a single lattice volume and at a single lattice spacing. Bearing these caveats in mind, the isovector nuclear momentum fractions were calculated precisely and found to 3 (x) of the nNNPDF2.0 isovector PDF in 3 He to that in the proton [31], as well as the same distribution with the LQCD moment constraint imposed into the global analysis as described in the text. 68% confidence intervals are shown.
be similar to that of the proton. In particular, the ratios x While in its early stages, this work emphasizes the utility of LQCD in constraining less well-measured aspects of partonic structure in an analogous way to how LQCD inputs have been used to constrain the proton transversity PDFs [82]. Future calculations at the physical quark masses will consider higher moments of nuclear PDFs (or even directly study their x dependence) for a wider range of nuclei and provide a complete flavor decomposition. Calculations will also quantitatively address the full set of systematic uncertainties.
Acknowledgements  [83], Qlua [84], QUDA [85,86], QDP-JIT [87] and QPhiX [88] software libraries were used in data production and analysis. WD, DJM and PES acknowledge support from the U.S. DOE grant DE-SC0011090. WD is also supported within the framework of the TMD Topical Collaboration of the U.S. DOE

MATRIX ELEMENT CALCULATION
The Euclidean finite difference form of the twist-two operator that determines the quark momentum fraction is The combination T ( ) = 1 irreducible representation of the hypercubic group [54], is used in this work. To compute matrix elements of this operator, the compound propagator technique is generalized from that previously used in Refs. [42,46]. The operator insertion point is used as a sequential source, and three-point correlation functions are formed by first calculating a (smeared) point-to-all quark propagator extending from the hadronic/nuclear source to the operator insertion point and subsequently calculating additional quark propagators from the operator insertion point to the sink. Since the finite difference form of the operator contains shifts, three different sequential inversions are utilized. Taking the appropriate linear combinations of displaced sources to implement T ( ) results in fixed-order background-field compound propagators that include the operator insertions throughout spacetime. These compound propagators are used to construct the two-point correlation functions in Eq. (4). The background-field two-point correlation functions G h (t; λ) have a spectral representation as a sum of exponentials. This, in turn, determines the full t-dependent form of the ratio of two-point correlation functions in zero and non-zero background fields, R h (t) (defined in Eq. (5)) in terms of the eigenenergies, interpolating operator overlap factors, and ground-and excited-state matrix elements of the lattice operator in Eq. (S1). The ground-state matrix elements of interest for each nuclear system can thus be extracted by fitting LQCD results for R h (t) to the form arising from the spectral representations.
Zero background field two-point correlation function have the spectral representation whereT ( ) = x T ( ) (x, 0) and thermal effects from the finite temporal extent of the lattice have been neglected, see Refs. [42,43,58] for further discussions. The spectral representation for R ss h (t) follows from inserting Eqs. (S2) and (S3) into Eq. (5) and can be expressed as  are used to fix the free parameters { n T ( ) n , E n , Z s n (Z s n ) * , N ss n }, for ss ∈ {SS, SP}. It is noteworthy that if this spectral representation is truncated at N states = 1, then the second sum in Eq. (S4) vanishes and R ss h (t) is independent of N ss 0 . For N states > 1, R ss h (t) is similarly independent from a linear combination of the N ss n that can be used to eliminate one redundant parameter, chosen here to be N ss 0 . Further, if the sums and finite-difference derivatives in Eq. (S3) and Eq. (S4) were replaced by continuum integrals and derivatives, then R ss h (t) would be independent of N ss n , suggesting that fits might not be sensitive to N ss n in practice if lattice artifacts are small. In this work, fits are performed both using Eq. (S4) and also neglecting these lattice artifacts by setting N ss n = 0 in Eq. (S4). The Aikiake Information Criterion (AIC) is used to select whether fits with or without these lattice artifacts are preferred for each choice of N states and fit range that is considered. In all cases, fits without these lattice artifacts are preferred by the AIC and used for subsequent analysis.
Care must be taken in order to account for systematic uncertainties associated with the fit-range choice and  S4)). In this work, an automated fitting procedure analogous to the procedure described in Ref. [61] is used to assess these systematic fitting uncertainties. Results for SS and SP interpolating operators are fit simultaneously and maximum source-sink separations {t SS max , t SP max } are chosen so that the signal-to-noise ratio of each correlation function is always greater than a tolerance of 0.5. A minimum separation t min min = 4 is required for the spectral representation to be well-defined for three-point correlation functions for the action and operator considered here. Within these constraints, N max = 200 choices of minimum separations for both sources (t SS min , t SP min ) are chosen randomly. For each fit range, a one-state fit is performed simultaneously for both sources. A two-state fit is subsequently performed, and if the two-state fit results in an AIC score of ∆AIC≤ −0.5, then the two-state fit is accepted. This procedure is repeated until an N state -state fit is rejected, at which point an (N state −1)-state fit is used. For the correlation functions analyzed in this work, this procedure results in accepted fits with N state ∈ {1, 2, 3}. Correlated χ 2 -minimization is used for all fits with covariance matrices regularized using optimal shrinkage [89,90]. The χ 2 depends linearly on the overlap factors, so these parameters are solved for using variable-projection methods [91,92]. The excited-state energy gaps are then determined using nonlinear optimization (Nelder-Mead and gradient-based Newton solvers from the Julia package Optim are used and verified to reproduce the same minimum energy gaps within an absolute tolerance of 10 −5 ). The fit is then repeated for N boot = 200 bootstrap samples, and the marginalized parameter uncertainty on the ground-state matrix element is determined using rank-based bootstrap confidence interval estimation [93] to add robustness to outlier bootstrap samples. Various checks on the fit result are then performed: an uncorrelated fit must reproduce the result within a tolerance of 5-sigma, the bootstrap median must reproduce the mean within a tolerance of 2-sigma, and the goodness-of-fit must be below a tolerance of χ 2 /N dof < 2. All fit results passing these criteria are included in an ensemble of acceptable fits. A weighted average of these acceptable fits is used to determine the final central value and total systematic uncertainty, where the (in principle arbitrary -see Ref. [94] for a discussion of this in a Bayesian framework) weights are taken to be the ratio of the p-value to the variance of each fit as in Ref. [90]. Final uncertainties are obtained by adding in quadrature the statistical uncertainty of the highest weight fit with the systematic uncertainty obtained from the weighted average. The fitting procedure is fully specified by the parameters and tolerance values above, as well as by a random seed for bootstrap resampling that is fixed to allow correlated ratios of matrix elements for different hadrons to be formed. Fig. S1 shows summaries of the fits of R 3pt h for each state studied in this work using the approach described above. It is also convenient to define the further ratio for a nuclear state h: which determines the ratio of the nuclear momentum fractions to that of the proton. Fig. S2 shows results for R h (t) including results obtained by fitting R h (t) and R p (t) independently as described above, using correlated bootstrap resampling to determine ratios of the ground-state matrix elements from all successful pairs of fit ranges, and taking a weighted average of the results. in principle independent of the matching scale p 2 = µ 2 0 , in practice there are discretization artifacts which arise as contamination in the form of dependence on the hypercubic invariants p [2n] ≡ µ p 2n µ , resulting in the classic jellyfishbone structure shown in Fig. S3.p µ can be expanded in a Taylor series in p [2n] , so either p [2n] orp [2n] may be used to perform the fit; p [2n] is used for this analysis, but consistent results are obtained by fitting top [2n] .
To obtain the renormalization factor, the hypercubic artifacts which depend on p [2n] for n > 1 are fit using the one-window-fit approach, detailed in Ref. [95]. Running terms and remaining artifacts which depend only on p 2 are then fit separately. To allow a quantification of systematic uncertainties, the fits are performed over a range of windows, with a range of functional forms. In particular, artifacts of the form are considered, where the H i denote hypercubic terms and the R i are contributions from running. Terms are grouped by their order in a 2 and by whether or not they contain logarithmic corrections. The full functional form is truncated at order a 4 . Adding further logarithmic terms into the functional form does not increase the fit quality. A fit form F , constructed from these components, is chosen for a specific data window to maximize the goodness of fit while preventing overfitting using the following procedure. Given a window of momenta, the first fit form is initialized to be F (1) = 0. Given a fit form F (n) , the subsequent form F (n+1) is determined by considering all possible forms which can be built from F (n) using only the terms X j which are not currently present in F (n) , and X ∈ {H, R} is as appropriate for the fit to the hypercubic or running artifacts. A fit form F u−d at quark masses corresponding to m π = 806 MeV is extrapolated to the physical masses using the assumption of weak mass dependence of the short-distance two-body counterterms in nuclear EFT. The isovector twist-two operators in Eq. (1) match on to hadronic operators in nuclear EFT as [63] O µ1...µn → x n (p) u−d v µ1 . . . v µn N † τ 3 N (1 + α 3,n N † N ) + . . .
where N is a nucleon field, v is the nuclear velocity, and the ellipsis denotes higher order nucleonic and pionic operators. Defining the nuclear factor G 3 (N, Z) = N, Z|N † τ 3 N N † N |N, Z , the two-nucleon counterterm α 3,2 relates the nuclear and proton momentum fractions as [63] This LEC can be determined from the LQCD results for 3 He most precisely by re-expressing it in terms of the quantities in Table I as Since renormalization effects cancel in the ratio, it is more precisely determined than the individual momentum fractions themselves and, with a naive error propagation, computing α 3,2 G 3 ( 3 He) via Eq. (S21) rather than Eq. (S20) achieves smaller uncertainties. The matching of the LEC, and extrapolation to physical quark masses proceeds in the following steps: