Lattice QCD constraints on the parton distribution functions of 3 He
William Detmold, Marc Illa, David J. Murphy, Patrick Oare, Kostas Orginos, Phiala E. Shanahan, Michael L. Wagman, Frank Winter
MMIT-CTP/5234, ICCUB-20-019, FERMILAB-PUB-20-466-T
Lattice QCD constraints on the parton distribution functions of He William Detmold, Marc Illa, David J. Murphy, Patrick Oare, Kostas Orginos,
3, 4
Phiala E. Shanahan, Michael L. Wagman, and Frank Winter (NPLQCD collaboration) Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Departament de F´ısica Qu`antica i Astrof´ısica and Institut de Ci`encies del Cosmos,Universitat de Barcelona, Mart´ı Franqu`es 1, E08028-Spain Department of Physics, College of William and Mary, Williamsburg, VA 23187-8795, USA Jefferson Laboratory, 12000 Jefferson Avenue, Newport News, VA 23606, USA Fermi National Accelerator Laboratory, Batavia, IL 60510, USA
The fraction of the longitudinal momentum of He that is carried by the isovector combinationof u and d quarks is determined using lattice QCD for the first time. The ratio of this combinationto that in the constituent nucleons is found to be consistent with unity at the few-percent levelfrom calculations with quark masses corresponding to m π ∼
800 MeV, extrapolated to the physicalquark masses. This constraint is consistent with, and significantly more precise than, determinationsfrom global nuclear parton distribution function fits. Including the lattice QCD determination ofthe momentum fraction in the nNNPDF global fitting framework results in the uncertainty on theisovector momentum fraction ratio being reduced by a factor of 2.5, and thereby enables a moreprecise extraction of the u and d parton distributions in He.
A central pillar of our understanding of the internalstructure of strongly interacting hadronic and nuclearsystems is knowledge of their partonic structure as ac-cessed in deep-inelastic scattering (DIS) experiments andother hard processes. Since the 1960s, such experimentshave revealed the longitudinal momentum distributionsof quarks and gluons in a fast moving proton, known col-lectively as parton distribution functions (PDFs). Thesimplest PDFs, q ( x, µ ) (and g ( x, µ )), describe the prob-ability of a quark of flavor q (or gluon g ) carrying a frac-tion x of the longitudinal momentum of the struck pro-ton at a renormalization scale µ . In 1983, the EuropeanMuon Collaboration (EMC) [1] observed that the par-tonic structure of nuclei differs substantially from thatof the constituent protons and neutrons, a landmark inthe development of nuclear physics [2–6]. Since the DISprocesses observed in the EMC experiments were at veryhigh energy, and the binding energy of a nucleus is smallin comparison to its mass, the appearance and size of theEMC effect was surprising at the time. Interest in theEMC effect has been rekindled by recent data from SLACand JLab [7–11] on EMC ratios for light nuclei. Not onlyhave these data provided precise determinations of theEMC effect for nuclei with small atomic number A , butthey have revealed a correlation between the strength ofthe EMC effect and so called “short range correlations”[12, 13].In addition to experimental investigations, theoreticalcalculations of the partonic structure of hadrons and nu-clei from the Standard Model can have important im-pact on our understanding of the structure of matter. For example, Standard Model calculations of nuclear par-tonic structure would reveal the QCD origin of the EMCeffect as well as aid in the flavor-separation of protonPDFs. Parton distributions are inherently rooted in thestrong interaction dynamics of QCD and cannot be de-termined using perturbative methods. Since the seminalworks of Refs. [14, 15], lattice Quantum Chromodynam-ics (LQCD) calculations have addressed the simplest as-pects of the parton distributions of the proton, notablydetermining the first few Mellin moments of the unpo-larized, polarized, and transversity quark distributions[16], as well as their gluonic analogues [17–21]. Recently,efforts have been made to extend these studies to thefull x -dependence of the proton PDFs [16, 22–24]. Morecomplicated extensions of partonic structure, such asgeneralized parton distribution functions and transverse-momentum dependent parton distribution functions ofthe proton, have also been studied using LQCD [25–30].In this letter, the partonic structure of light nuclei isstudied in LQCD for the first time through an investiga-tion of the isovector quark momentum fractions (the firstmoments of the corresponding isovector PDFs) of theproton, diproton, and He. At the heavier-than-physicalquark masses used in this LQCD study, percent-level nu-clear effects are resolved in the momentum fraction of He. After an extrapolation to the physical quark masses,these calculations provide a constraint on the isovectormomentum fraction that is used as an additional inputinto the nNNPDF2.0 [31] global nuclear PDF analysisframework. Since the isovector combination of nuclearPDFs is poorly determined from experiment, this LQCD a r X i v : . [ h e p - l a t ] S e p constraint significantly reduces the uncertainties on the He PDFs and thereby improves knowledge of nuclearstructure.
LQCD methodology —
The existence of strong in-teractions between quarks and gluons necessitates the useof LQCD for calculations of the partonic structure of nu-clei. The calculations presented here are performed us-ing a single ensemble of gauge-field configurations gener-ated with a L¨uscher-Weisz gauge action [32] with N f = 3degenerate light-quark flavors with the clover-improvedWilson fermion action [33], and quark masses tuned toproduce a pion mass of m π = 806 MeV. The lattice ge-ometry is L × T = 32 ×
48, and the lattice spacing isdetermined to be a ∼ .
145 fm from Υ spectroscopy [34].This ensemble, and two others with different spacetimevolumes have previously been used to study the spectrum[34, 35] and properties [36–46] of light nuclei up to atomicnumber A = 4. The multi-volume spectroscopy studiesshow that the pp and He states that are investigatedhere are bound systems with infinite volume energies be-low threshold. Consequently, matrix elements in thesestates are expected to receive only exponentially smallfinite volume effects, O ( e − κL , e − m π L ), that will be ne-glected in this work [47–53].The Mellin moments of the unpolarized isovectorquark PDFs, q ( h )3 ( x, µ ) = u ( h ) ( x, µ ) − d ( h ) ( x, µ ), in ahadronic or nuclear state h , defined as (cid:104) x n (cid:105) ( h ) u − d ( µ ) ≡ (cid:82) − dx x n q ( h )3 ( x, µ ), are determined from matrix elementsof twist-two operators as (cid:104) h |O µ ...µ n | h (cid:105) ≡ (cid:104) h | qτ γ { µ ( i ←→ D µ ) . . . ( i ←→ D µ n } ) q | h (cid:105) = (cid:104) x n (cid:105) ( h ) u − d ( µ ) p { µ . . . p µ n } , (1)where p is the momentum of the state h , τ is a Paulimatrix in flavor space, ←→ D µ = ( −→ D µ − ←− D µ ) / D µ is the gauge covariant derivative, and { . . . } indicatessymmetrization and trace-subtraction of the enclosed in-dices. The above operators are constructed to transformirreducibly under the Lorentz group, but the hypercu-bic spacetime lattice used in the LQCD calculations re-duces these symmetries, in general inducing mixing be-tween operators of different Lorentz spin. In particu-lar, the two-index operators that determine the isovectorquark momentum fraction, (cid:104) x (cid:105) ( h ) u − d , subduce to operatorsin two different irreducible representations of the hyper-cubic group. In this work, matrix elements of an Eu-clidean operator in the τ (3)1 representation [54] are com-puted, namely T = 1 √ T − T ) , with T µν = qτ γ { µ ←→ D ν } q , (2)where γ ν is also Euclidean. With a lattice regulator,this operator is discretized as a covariant finite differ-ence 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 zerothree-momentum determine the momentum fraction as (cid:104) h | T | h (cid:105) = (cid:104) x (cid:105) ( h ) u − d M h / √ T (MS) ( µ ) = R MS / RI (cid:48) MOM ( µ, µ ) Z RI (cid:48) MOM ( µ , a ) T ( a ) , (3)where the renormalization coefficient Z RI (cid:48) MOM ( µ , a )is defined non-perturbatively in a regularization-independent momentum-subtraction scheme [55] at ascale µ and then matched to MS through the three-loopperturbative coefficient R MS / RI (cid:48) MOM ( µ, µ ) [56, 57], asdetailed in the Supplementary Material. For µ = 2 GeV, R MS / RI (cid:48) MOM ( µ, µ ) Z RI (cid:48) MOM ( µ , a ) = 0 . T -compound propagators are computedfrom an average of N src = 24 source points randomlydistributed on N cfg = 2290 gauge-field configurations for N B = 5 different background field strengths. These com-pound propagators are then used to construct baryontwo-point correlation functions, G h ( t ; λ ) = (cid:88) x (cid:68) (cid:12)(cid:12)(cid:12) χ h ( x , t ) χ † h (0) (cid:12)(cid:12)(cid:12) (cid:69) λ , (4)where λ is the T -background field strength, χ h is an in-terpolating field for states with the quantum numbers ofthe hadron or nucleus h , and spinor indices on the in-terpolating operators are suppressed. Correlation func-tions are constructed from Gaussian-smeared source in-terpolating operators [59], while the sink interpolatingoperators are either smeared or point-like and the multi-baryon contractions are performed using the techniquesof Refs. [60]. This quantity contains responses to thefield up to O ( λ N Q ), with N Q being the number of va-lence quarks in the state. The linear response of thisbackground-field two-point function, G h ( t ; λ ) | O ( λ ) , is de-termined by the matrix element of T . This term canbe extracted exactly from the computed set of fixed-order background field correlation functions with N Q field strengths [42, 43].Combining the linear response of the two-point corre-lation function with the zero-field correlation function, itis straightforward to show that the ratio R h ( t ) = G h ( t ; λ ) | O ( λ ) G h ( t ; 0) − G h ( t − a ; λ ) | O ( λ ) G h ( t − a ; 0) (5)is related to matrix elements of T through the spec-tral representation of each term in Eq. (5), in particular □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ FIG. 1. The effective matrix element, Eq. (5), associatedwith the isovector quark momentum fractions of the proton, pp and He. Blue (orange) points, labelled SS (SP), showresults for interpolating operators with smeared sources andsmeared (point-like) sinks. For each effective matrix element,points are shown for t ≤ t max , where t max is the minimum t where the signal-to-noise ratio of G h ( t + a ; λ ) | O ( λ ) is less than0.5. Colored bands show the highest weight fit to the com-bined dataset and the shaded gray bands show the weightedaverage of all accepted fits and the total statistical plus fittingsystematic uncertainties. asymptoting as R h ( t ) t →∞ −→ (cid:104) h | T | h (cid:105) , (6)with exponentially vanishing contamination at earlytimes that involves excited-state overlap factors and tran-sition 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 rangesand models analogous to the procedure described for two-point correlation functions in Ref. [61]. In summary, inanalyzing R h ( t ) to extract the momentum fractions, thefull t dependence that results from the spectral decom-position of each term in Eq. (5) is fit, and combined fitsto two- and three-point correlation functions are used toconstrain the relevant energies, overlap factors, and ma-trix elements. All possible choices of fit ranges and up □□ ○○ ▽▽ □○▽ ○○ ▽▽ FIG. 2. Left: Renormalized isovector momentum fractionsfor h ∈ { p, pp, He } at a scale of µ = 2 GeV. Right: Ratiosof the isovector nuclear momentum fractions to that of theconstituent nucleons. p pp He (cid:104) x (cid:105) ( h ) u − d . . . (cid:16) AZ − N (cid:17) (cid:104) x (cid:105) ( h ) u − d / (cid:104) x (cid:105) ( p ) u − d — 1 . . p , pp and He, calculated at m π = 806 MeV in MS-scheme at µ = 2GeV. The first uncertainty combines LQCD statistical andsystematic uncertainties and the second uncertainty is fromoperator renormalization. The correlated ratios of the isovec-tor momentum fraction in nuclei to those in the constituentnucleons, in which the renormalization constants and theiruncertainties cancel, are also given. to 4 states contributing to the spectral decompositionsare considered using a model selection process describedin the Supplementary Material. A weighted average overfits from all acceptable fit ranges is used to define ground-state energy results, including systematic uncertaintiesfrom fit range and model variation. Results are shown inFig. 1 for the proton, diproton and He.
Results and Discussion —
The extracted values ofthe isovector quark momentum fractions for p, pp, He atquark masses corresponding to m π = m K = 806 MeV areshown in Tab. I and displayed graphically in Fig. 2. Theuncertainties are separated into those from the LQCDcalculation of the bare matrix elements, and the (larger)uncertainty from the renormalization and matching tothe MS scheme. The proton isovector momentum frac-tion is consistent with other LQCD extractions at similarvalues of the quark masses [62] given the different renor-malization procedures and lattice spacings. The pp and He momentum fractions are determined with O (5%) un-certainties and are found to be approximately consistentwith those of the constituent nucleons. The ratios of thenuclear momentum fractions to that of the proton areindependent of operator renormalization to O ( α s ), andare determined at few-percent precision even for He.In Refs. [63–65], nuclear effective field theory (EFT)was used to study nuclear effects in PDF moments. Inparticular, it was shown that the leading source of sucheffects is the two-nucleon correlations that couple to thetwist-two operators defining the PDF moments. In termsof the parameters defined in that work, nuclear effectsin the isovector momentum fraction are encapsulated inthe low energy constant (LEC) α , and nuclear fac-tor G ( He); their product is bounded as α , G ( He) =0 . µ = 2 GeV from the numerical calcula-tions presented here (see the Supplementary Material fordetails). While the quark momentum fractions them-selves have nonanalytic dependence on the quark masses[66–68], this two-body LEC is expected to be relativelyinsensitive to variation of the quark masses, as seen forthe the analogous two-body contribution in the np → dγ [39] and pp → de + ν e [42, 69] processes. This relativemass-independence assumption allows an extrapolationto the physical quark masses: a naive estimate is givenby taking the central value determined at m π = 806MeV and inflating the uncertainty by 50% to accountfor possible quark-mass dependence as well as the effectsof the nonzero lattice spacing and finite volume (this un-certainty is estimated based on the mass dependence seenfor the analogous two-body LECs in Refs. [39, 42, 69]).This extrapolated value can be combined with the physi-cal value of the nucleon momentum fraction, (cid:104) x (cid:105) ( p ) u − d =0 . µ = 2 GeV from the nNNPDF2.0 analy-sis [31], to determine the isovector momentum fractionratio 3 (cid:104) x (cid:105) ( He) u − d / (cid:104) x (cid:105) ( p ) u − d | LQCD = 1 . F structurefunctions of nuclei to that of the deuteron [3, 5, 6]. Ad-ditional constraints are especially valuable in the contextof the intriguing question as to whether there is flavor-dependence to the EMC effect. Such flavor dependencehas been conjectured in models of QCD [70–75] and inEFT [63–65] and is included in recent data-driven analy-ses of experimental results [76, 77] and provides a poten-tial explanation of the NuTeV anomaly in sin θ W [78].Fig. 3 shows the constraint on the isovector momen-tum fraction ratio for He obtained from the resultspresented here, compared with the constraints on theisovector and isoscalar momentum fraction ratios fromthe recent nNNPDF2.0 [31] global nuclear PDF fits. ThenNNPDF2.0 ellipse is generated by combining the MonteCarlo replica sets for the bound proton PDFs in Heappropriately to form the PDFs of He (under the as-
FIG. 3. The ratio of the isovector momentum fractions of He and p determined in this work compared to constraints onthe isovector and isoscalar momentum fraction ratios from thenNNPDF2.0 [31] global analysis before and after the LQCDconstraint is imposed. Both axes are normalized to unity inthe absence of nuclear effects. The LQCD constraint on theisovector ratio at m π = 806 MeV is also displayed. In allcases, 68% confidence intervals are shown. sumption that the nuclear effects vary slowly with A ). Inthis way, correlations between the He and proton PDFsare accounted for. For the isovector combination, the68% confidence interval is 3 (cid:104) x (cid:105) ( He) u − d / (cid:104) x (cid:105) ( p ) u − d | nNNPDF2.0 =1 . (cid:104) x (cid:105) ( He) u − d / (cid:104) x (cid:105) ( p ) u − d | nNNPDF2.0+LQCD = 1 . He to thatof the constituent nucleons, with and without the impo-sition of the LQCD constraint. As can be seen from thereduced uncertainties in Figs. 3 and 4, LQCD calcula-tions such as those presented here, as well as new experi-mental constraints [80, 81], can significantly improve ourknowledge of the flavor dependence of nuclear PDFs.
Summary —
In this work, the isovector momentumfractions of the proton, diproton and He systems havebeen determined using LQCD, complementing a previ-ous study of the gluon momentum fraction on the sameensemble [45]. These calculations were performed at asingle set of unphysical SU(3)-symmetric values for thequark masses corresponding to m π = 806 MeV, and ina single lattice volume and at a single lattice spacing.Bearing these caveats in mind, the isovector nuclear mo-mentum fractions were calculated precisely and found to FIG. 4. The ratio R ( He) ( x ) = 3 q ( He)3 ( x ) /q ( p )3 ( x ) of thenNNPDF2.0 isovector PDF in He to that in the proton [31],as well as the same distribution with the LQCD moment con-straint imposed into the global analysis as described in thetext. 68% confidence intervals are shown. be similar to that of the proton. In particular, the ra-tios (cid:104) x (cid:105) ( pp ) u − d / (cid:104) x (cid:105) ( p ) u − d = 1 . (cid:104) x (cid:105) ( He) u − d / (cid:104) x (cid:105) ( p ) u − d =1 . He result to global analysesof nuclear PDFs, providing important constraints on theflavor decomposition of nuclear PDFs that are comple-mentary to those obtained from experiment.While in its early stages, this work emphasizes the util-ity of LQCD in constraining less well-measured aspectsof partonic structure in an analogous way to how LQCDinputs have been used to constrain the proton transver-sity PDFs [82]. Future calculations at the physical quarkmasses will consider higher moments of nuclear PDFs(or even directly study their x dependence) for a widerrange of nuclei and provide a complete flavor decompo-sition. Calculations will also quantitatively address thefull set of systematic uncertainties. Acknowledgements —
We thank Juan Rojo forcomments and substantial help with the comparison tothe nNNPDF2.0 global fits, and Silas Beane, Jiunn-WeiChen, Zohreh Davoudi, Assumpta Parre˜no, Martin Sav-age and Brian Tiburzi for insightful discussions. Thisresearch used resources of the Oak Ridge LeadershipComputing Facility at the Oak Ridge National Labo-ratory, which is supported by the Office of Science ofthe U.S. Department of Energy under Contract numberDE-AC05-00OR22725, as well as facilities of the USQCDCollaboration, which are funded by the Office of Sci-ence of the U.S. Department of Energy. This researchused resources of the National Energy Research Scien-tific Computing Center (NERSC), a U.S. Departmentof EnergyOffice of Science User Facility operatedunderContract No. DE-AC02-05CH11231. The Chroma [83],Qlua [84], QUDA [85, 86], QDP-JIT [87] and QPhiX [88]software libraries were used in data production and anal-ysis. WD, DJM and PES acknowledge support from the U.S. DOE grant DE-SC0011090. WD is also supportedwithin the framework of the TMD Topical Collabora-tion of the U.S. DOE Office of Nuclear Physics, and bythe SciDAC4 award DE-SC0018121. PES is additionallysupported by the National Science Foundation under CA-REER Award 1841699 and under EAGER grant 2035015,by the U.S. DOE Early Career Award DE-SC0021006,by a NEC research award, by the Carl G and ShirleySontheimer Research Fund. MI is supported by the Uni-versitat de Barcelona through the scholarship APIF, bythe Spanish Ministerio de Econom´ıa y Competitividad(MINECO) under the project No. MDM-2014-0369 ofICCUB (Unidad de Excelencia Mar´ıa de Maeztu) andwith additional European FEDER funds under the con-tract FIS2017-87534-P. KO was supported in part byU.S. DOE grant DE-FG02-04ER41302 and in part bythe Jefferson Science Associates, LLC under U.S. DOEContract DE-AC05-06OR23177. This manuscript hasbeen authored by Fermi Research Alliance, LLC underContract No. DE-AC02-07CH11359 with the U.S. De-partment of Energy, Office of Science, Office of High En-ergy Physics. The authors thank Robert Edwards, B´alintJo´o, and members of the NPLQCD collaboration for gen-erating and allowing access to the ensembles used in thisstudy, as well as for helpful discussions. [1] J. J. Aubert et al. (European Muon), Phys. Lett.
B123 ,275 (1983).[2] M. Arneodo, Phys. Rept. , 301 (1994).[3] D. F. Geesaman, K. Saito, and A. W. Thomas, Ann.Rev. Nucl. Part. Sci. , 337 (1995).[4] G. Piller and W. Weise, Phys. Rept. , 1 (2000).[5] P. R. Norton, Rept. Prog. Phys. , 1253 (2003).[6] O. Hen, D. W. Higinbotham, G. A. Miller, E. Piasetzky,and L. B. Weinstein, Int. J. Mod. Phys. E , 1330017(2013).[7] N. Fomin, J. Arrington, R. Asaturyan, F. Benmokhtar,W. Boeglin, P. Bosted, A. Bruell, M. H. S. Bukhari,M. E. Christy, E. Chudakov, et al. , Phys. Rev. Lett. ,092502 (2012).[8] O. Hen, E. Piasetzky, and L. B. Weinstein, Phys. Rev. C85 , 047301 (2012), arXiv:1202.3452 [nucl-ex].[9] L. L. Frankfurt, M. I. Strikman, D. B. Day, andM. Sargsian, Phys. Rev.
C48 , 2451 (1993).[10] K. S. Egiyan et al. (CLAS), Phys. Rev.
C68 , 014313(2003), arXiv:nucl-ex/0301008 [nucl-ex].[11] K. S. Egiyan et al. (CLAS Collaboration), Phys. Rev.Lett. , 082501 (2006).[12] L. B. Weinstein, E. Piasetzky, D. W. Higinbotham,J. Gomez, O. Hen, and R. Shneor, Phys. Rev. Lett. , 052301 (2011).[13] O. Hen, E. Piasetzky, and L. B. Weinstein, Phys. Rev.C , 047301 (2012).[14] G. Martinelli and C. T. Sachrajda, Nucl. Phys. B306 ,865 (1988).[15] G. Martinelli and C. T. Sachrajda, Nucl. Phys.
B316 ,355 (1989). [16] H.-W. Lin et al. , Prog. Part. Nucl. Phys. , 107 (2018),arXiv:1711.07916 [hep-ph].[17] R. Horsley, R. Millo, Y. Nakamura, H. Perlt, D. Pleiter,P. E. L. Rakow, G. Schierholz, A. Schiller, F. Winter,and J. M. Zanotti (UKQCD, QCDSF), Phys. Lett.
B714 ,312 (2012), arXiv:1205.6410 [hep-lat].[18] C. Alexandrou, M. Constantinou, K. Hadjiyiannakou,K. Jansen, H. Panagopoulos, and C. Wiese, (2016),arXiv:1611.06901 [hep-lat].[19] Y.-B. Yang, M. Gong, J. Liang, H.-W. Lin, K.-F. Liu,D. Pefkou, and P. Shanahan, Phys. Rev.
D98 , 074506(2018), arXiv:1805.00531 [hep-lat].[20] P. E. Shanahan and W. Detmold, Phys. Rev.
D99 ,014511 (2019), arXiv:1810.04626 [hep-lat].[21] P. E. Shanahan and W. Detmold, Phys. Rev. Lett. ,072003 (2019), arXiv:1810.07589 [nucl-th].[22] X. Ji, Phys. Rev. Lett. , 262002 (2013),arXiv:1305.1539 [hep-ph].[23] K. Cichy and M. Constantinou, Adv. High Energy Phys. , 3036904 (2019), arXiv:1811.07248 [hep-lat].[24] X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang, and Y. Zhao,(2020), arXiv:2004.03543 [hep-ph].[25] P. Hagler, Phys. Rept. , 49 (2010), arXiv:0912.5483[hep-lat].[26] B. Yoon, M. Engelhardt, R. Gupta, T. Bhattacharya,J. R. Green, B. U. Musch, J. W. Negele, A. V. Pochinsky,A. Schfer, and S. N. Syritsyn, Phys. Rev. D , 094508(2017), arXiv:1706.03406 [hep-lat].[27] P. Shanahan, M. Wagman, and Y. Zhao, Phys. Rev. D , 014511 (2020), arXiv:2003.06063 [hep-lat].[28] C. Alexandrou, K. Cichy, M. Constantinou, K. Had-jiyiannakou, K. Jansen, A. Scapellato, and F. Steffens,(2020), arXiv:2008.10573 [hep-lat].[29] H.-W. Lin, (2020), arXiv:2008.12474 [hep-ph].[30] Q.-A. Zhang et al. (Lattice Parton), (2020),arXiv:2005.14572 [hep-lat].[31] R. Abdul Khalek, J. J. Ethier, J. Rojo, and G. vanWeelden, (2020), arXiv:2006.14629 [hep-ph].[32] M. L¨uscher and P. Weisz, Commun. Math. Phys. , 59(1985), [Erratum: Commun. Math. Phys.98,433(1985)].[33] B. Sheikholeslami and R. Wohlert, Nucl. Phys. B259 ,572 (1985).[34] S. R. Beane, E. Chang, S. D. Cohen, W. Detmold, H. W.Lin, T. C. Luu, K. Orginos, A. Parre˜no, M. J. Sav-age, and A. Walker-Loud (NPLQCD), Phys. Rev.
D87 ,034506 (2013), arXiv:1206.5219 [hep-lat].[35] M. L. Wagman, F. Winter, E. Chang, Z. Davoudi,W. Detmold, K. Orginos, M. J. Savage, and P. E. Shana-han, Phys. Rev. D , 114510 (2017), arXiv:1706.06550[hep-lat].[36] S. R. Beane et al. , (2017), arXiv:1705.09239 [hep-lat].[37] S. R. Beane et al. (NPLQCD), Phys. Rev. C88 , 024003(2013), arXiv:1301.5790 [hep-lat].[38] S. R. Beane, E. Chang, S. Cohen, W. Detmold, H. W.Lin, K. Orginos, A. Parre˜no, M. J. Savage, andB. C. Tiburzi, Phys. Rev. Lett. , 252001 (2014),arXiv:1409.3556 [hep-lat].[39] S. R. Beane, E. Chang, W. Detmold, K. Orginos,A. Parreo, M. J. Savage, and B. C. Tiburzi (NPLQCD),Phys. Rev. Lett. , 132001 (2015), arXiv:1505.02422[hep-lat].[40] E. Chang, W. Detmold, K. Orginos, A. Parre˜no, M. J.Savage, B. C. Tiburzi, and S. R. Beane (NPLQCD),Phys. Rev.
D92 , 114502 (2015), arXiv:1506.05518 [hep- lat].[41] W. Detmold, K. Orginos, A. Parre˜no, M. J. Savage, B. C.Tiburzi, S. R. Beane, and E. Chang, Phys. Rev. Lett. , 112301 (2016), arXiv:1508.05884 [hep-lat].[42] M. J. Savage, P. E. Shanahan, B. C. Tiburzi, M. L. Wag-man, F. Winter, S. R. Beane, E. Chang, Z. Davoudi,W. Detmold, and K. Orginos, Phys. Rev. Lett. ,062002 (2017), arXiv:1610.04545 [hep-lat].[43] B. C. Tiburzi, M. L. Wagman, F. Winter, E. Chang,Z. Davoudi, W. Detmold, K. Orginos, M. J. Savage,and P. E. Shanahan, Phys. Rev.
D96 , 054505 (2017),arXiv:1702.02929 [hep-lat].[44] P. E. Shanahan, B. C. Tiburzi, M. L. Wagman, F. Win-ter, E. Chang, Z. Davoudi, W. Detmold, K. Orginos,and M. J. Savage, Phys. Rev. Lett. , 062003 (2017),arXiv:1701.03456 [hep-lat].[45] F. Winter, W. Detmold, A. S. Gambhir, K. Orginos,M. J. Savage, P. E. Shanahan, and M. L. Wagman, Phys.Rev.
D96 , 094512 (2017), arXiv:1709.00395 [hep-lat].[46] E. Chang, Z. Davoudi, W. Detmold, A. S. Gambhir,K. Orginos, M. J. Savage, P. E. Shanahan, M. L. Wag-man, and F. Winter (NPLQCD), Phys. Rev. Lett. ,152002 (2018), arXiv:1712.03221 [hep-lat].[47] M. L¨uscher, Commun. Math. Phys. , 177 (1986).[48] M. L¨uscher, Commun. Math. Phys. , 153 (1986).[49] M. L¨uscher, Nucl. Phys. B , 531 (1991).[50] S. Beane, P. Bedaque, A. Parre˜no, and M. Savage, Phys.Lett. B , 106 (2004), arXiv:hep-lat/0312004.[51] Z. Davoudi and M. J. Savage, Phys. Rev. D , 114502(2011), arXiv:1108.5371 [hep-lat].[52] S. Knig and D. Lee, Phys. Lett. B , 9 (2018),arXiv:1701.00279 [hep-lat].[53] R. A. Briceo, M. T. Hansen, and A. W. Jackura, Phys.Rev. D , 114505 (2019), arXiv:1909.10357 [hep-lat].[54] M. Gockeler, R. Horsley, E.-M. Ilgenfritz, H. Perlt,P. E. L. Rakow, G. Schierholz, and A. Schiller, Phys.Rev. D54 , 5705 (1996), arXiv:hep-lat/9602029 [hep-lat].[55] G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa, andA. Vladikas, Nucl. Phys. B , 81 (1995), arXiv:hep-lat/9411010.[56] J. Gracey, Nucl. Phys. B , 242 (2003), arXiv:hep-ph/0306163.[57] J. Gracey, Nucl. Phys. B , 247 (2003), arXiv:hep-ph/0304113.[58] C. Bouchard, C. C. Chang, T. Kurth, K. Orginos,and A. Walker-Loud, Phys. Rev.
D96 , 014504 (2017),arXiv:1612.06963 [hep-lat].[59] M. Albanese et al. (APE), Phys. Lett. B , 163 (1987).[60] W. Detmold and K. Orginos, Phys. Rev.
D87 , 114512(2013), arXiv:1207.1452 [hep-lat].[61] S. Beane et al. , (2020), arXiv:2003.12130 [hep-lat].[62] D. Dolgov et al. , Nucl. Phys. B Proc. Suppl. , 303(2001), arXiv:hep-lat/0011010.[63] J.-W. Chen and W. Detmold, Phys. Lett. B , 165(2005).[64] J.-W. Chen, W. Detmold, J. E. Lynn, and A. Schwenk,Phys. Rev. Lett. , 262502 (2017), arXiv:1607.03065[hep-ph].[65] J. Lynn, D. Lonardoni, J. Carlson, J. Chen, W. Detmold,S. Gandolfi, and A. Schwenk, J. Phys. G , 045109(2020), arXiv:1903.12587 [nucl-th].[66] W. Detmold, W. Melnitchouk, J. W. Negele, D. B. Ren-ner, and A. W. Thomas, Phys. Rev. Lett. , 172001(2001), arXiv:hep-lat/0103006. [67] D. Arndt and M. J. Savage, Nucl. Phys. A , 429(2002), arXiv:nucl-th/0105045.[68] J.-W. Chen and X.-d. Ji, Phys. Lett. B , 107 (2001),arXiv:hep-ph/0105197.[69] NPLQCD collaboration (2020) in progress.[70] T. Uchiyama and K. Saito, Phys. Rev. C38 , 2245 (1988).[71] I. R. Afnan, F. R. P. Bissey, J. Gomez, A. T. Katramatou,W. Melnitchouk, G. G. Petratos, and A. W. Thomas,Phys. Lett.
B493 , 36 (2000), arXiv:nucl-th/0006003[nucl-th].[72] K. Saito, C. Boros, K. Tsushima, F. R. P. Bissey, I. R. Af-nan, and A. W. Thomas, Phys. Lett.
B493 , 288 (2000),arXiv:nucl-th/0008017 [nucl-th].[73] W. Bentz, I. C. Cloet, J. T. Londergan, and A. W.Thomas, Phys. Lett.
B693 , 462 (2010), arXiv:0908.3198[nucl-th].[74] I. C. Cloet, W. Bentz, and A. W. Thomas, Phys. Rev.Lett. , 182301 (2012), arXiv:1202.6401 [nucl-th].[75] A. J. Tropiano, J. J. Ethier, W. Melnitchouk,and N. Sato, Phys. Rev.
C99 , 035201 (2019),arXiv:1811.07668 [nucl-th].[76] B. Schmookler et al. (CLAS), Nature , 354 (2019).[77] E. P. Segarra, A. Schmidt, D. W. Higinbotham, E. Pi-asetzky, M. Strikman, L. B. Weinstein, and O. Hen,(2019), arXiv:1908.02223 [nucl-th].[78] I. C. Cloet, W. Bentz, and A. W. Thomas, Phys. Rev.Lett. , 252301 (2009), arXiv:0901.3559 [nucl-th].[79] R. D. Ball, V. Bertone, F. Cerutti, L. Del Debbio,S. Forte, A. Guffanti, N. P. Hartland, J. I. Latorre,J. Rojo, and M. Ubiali, Nucl. Phys. B , 608 (2012),arXiv:1108.1758 [hep-ph].[80] J. Mousseau et al. (MINERvA), Phys. Rev.
D93 , 071101(2016), arXiv:1601.06313 [hep-ex].[81] G. G. Petratos et al. , Jefferson Lab PAC37 Proposal ex-periment E12-10-103, 2010.[82] H.-W. Lin, W. Melnitchouk, A. Prokudin, N. Sato,and H. Shows, Phys. Rev. Lett. , 152502 (2018),arXiv:1710.09858 [hep-ph].[83] R. G. Edwards and B. Jo´o (SciDAC, LHPC, UKQCD),
Lattice field theory. Proceedings, 22nd International Sym-posium, Lattice 2004, Batavia, USA, June 21-26, 2004 ,Nucl. Phys. Proc. Suppl. , 832 (2005), [,832(2004)],arXiv:hep-lat/0409003 [hep-lat].[84] A. Pochinsky, “Qlua, https://usqcd.lns.mit.edu/qlua.”.[85] M. Clark, R. Babich, K. Barros, R. Brower, andC. Rebbi, Comput. Phys. Commun. , 1517 (2010),arXiv:0911.3191 [hep-lat].[86] R. Babich, M. A. Clark, and B. Joo, in
SC 10 (Super-computing 2010) (2010) arXiv:1011.0024 [hep-lat].[87] F. Winter, PoS
LATTICE2013 , 042 (2014).[88] B. Jo´o, D. D. Kalamkar, T. Kurth, K. Vaidyanathan,and A. Walden, in
High Performance Computing , editedby M. Taufer, B. Mohr, and J. M. Kunkel (SpringerInternational Publishing, Cham, 2016) pp. 415–427.[89] O. Ledoit and M. Wolf, Journal of Multivariate Analysis , 365 (2004).[90] E. Rinaldi, S. Syritsyn, M. L. Wagman, M. I. Buchoff,C. Schroeder, and J. Wasem, Phys. Rev. D99 , 074510(2019), arXiv:1901.07519 [hep-lat].[91] G. Golub and V. Pereyra, Inverse Problems , R1(2003).[92] B. W. R. Dianne P. O’Leary, Computational Optimiza-tion and Applications , 579 (2013).[93] A. C. Davison and D. V. Hinkley, Bootstrap Methods and their Application , Cambridge Series in Statistical andProbabilistic Mathematics (Cambridge University Press,1997).[94] W. I. Jay and E. T. Neil, (2020), arXiv:2008.01069[stat.ME].[95] B. Blossier, P. Boucaud, M. Brinet, F. De Soto,Z. Liu, V. Morenas, O. Pene, K. Petrov, andJ. Rodriguez-Quintero, Phys. Rev. D , 074506 (2011),arXiv:1011.2414 [hep-ph]. Supplementary Material
MATRIX ELEMENT CALCULATION
The Euclidean finite difference form of the twist-two operator that determines the quark momentum fraction is T ( (cid:3) ) µν ( x ) = 14 a (cid:16) q ( x ) τ γ { ν [ U µ } ( x ) q ( x + ˆ µ ) − U † µ } ( x − ˆ µ ) q ( x − ˆ µ )] − [ q ( x + ˆ µ ) U †{ µ ( x ) − q ( x − ˆ µ ) U { µ ( x − ˆ µ )] τ γ ν } q ( x ) (cid:17) . (S1)The combination T ( (cid:3) ) = √ ( T ( (cid:3) )33 − T ( (cid:3) )44 ), which belongs to the τ (3)1 irreducible representation of the hypercubicgroup [54], is used in this work. To compute matrix elements of this operator, the compound propagator techniqueis generalized from that previously used in Refs. [42, 46]. The operator insertion point is used as a sequentialsource, and three-point correlation functions are formed by first calculating a (smeared) point-to-all quark propagatorextending from the hadronic/nuclear source to the operator insertion point and subsequently calculating additionalquark propagators from the operator insertion point to the sink. Since the finite difference form of the operator containsshifts, three different sequential inversions are utilized. Taking the appropriate linear combinations of displaced sourcesto implement T ( (cid:3) ) results in fixed-order background-field compound propagators that include the operator insertionsthroughout spacetime. These compound propagators are used to construct the two-point correlation functions inEq. (4). The background-field two-point correlation functions G h ( t ; λ ) have a spectral representation as a sum ofexponentials. This, in turn, determines the full t -dependent form of the ratio of two-point correlation functions inzero and non-zero background fields, R h ( t ) (defined in Eq. (5)) in terms of the eigenenergies, interpolating operatoroverlap factors, and ground- and excited-state matrix elements of the lattice operator in Eq. (S1). The ground-statematrix elements of interest for each nuclear system can thus be extracted by fitting LQCD results for R h ( t ) to theform arising from the spectral representations.Zero background field two-point correlation function have the spectral representation G ss (cid:48) h ( t ; λ = 0) = (cid:88) n Z sn ( Z s (cid:48) n ) ∗ e − E n t , (S2)where { s, s (cid:48) } ∈ { S, P } specifies the source and sink smearing, E n is the energy of the n -th energy eigenstate, and Z sn is an overlap factor defined by Z sn = √ V (cid:104) n | χ sh (0) | (cid:105) , where V = a (cid:81) µ L µ is the (dimensionful) lattice volume and L µ is the extent of the lattice geometry in the µ direction. The corresponding background-field two-point functionsat O ( λ ) have the spectral representation G ss (cid:48) h ( t ; λ ) (cid:12)(cid:12)(cid:12) O ( λ ) = a t (cid:88) τ =0 (cid:88) n,m Z sn ( Z s (cid:48) m ) ∗ e − E n τ (cid:68) n (cid:12)(cid:12)(cid:12) ˜ T ( (cid:3) ) (cid:12)(cid:12)(cid:12) m (cid:69) = (cid:88) n Z sn ( Z s (cid:48) n ) ∗ t e − E n t (cid:68) n (cid:12)(cid:12)(cid:12) ˜ T ( (cid:3) ) (cid:12)(cid:12)(cid:12) n (cid:69) + (cid:88) n (cid:88) m (cid:54) = n Z sn ( Z s (cid:48) m ) ∗ a (cid:18) e − E n t − e ( E n − E m ) a + e − E m t − e ( E m − E n ) a (cid:19) (cid:68) n (cid:12)(cid:12)(cid:12) ˜ T ( (cid:3) ) (cid:12)(cid:12)(cid:12) m (cid:69) (S3)where ˜ T ( (cid:3) ) = (cid:80) x T ( (cid:3) ) ( 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 (cid:48) h ( t ) follows from inserting Eqs. (S2)and (S3) into Eq. (5) and can be expressed as R ss (cid:48) h ( t ) = (cid:88) n (cid:68) n (cid:12)(cid:12)(cid:12) T ( (cid:3) ) (cid:12)(cid:12)(cid:12) n (cid:69) Z sn ( Z s (cid:48) n ) ∗ (cid:18) ( t + a ) e − E n ( t + a ) (cid:80) k Z sk ( Z s (cid:48) k ) ∗ e − E k ( t + a ) − t e − E n t (cid:80) k Z sk ( Z s (cid:48) k ) ∗ e − E k t (cid:19) + (cid:88) n N ss (cid:48) n (cid:18) e − E n ( t + a ) (cid:80) k Z sk ( Z s (cid:48) k ) ∗ e − E k ( t + a ) − e − E n t (cid:80) k Z sk ( Z s (cid:48) k ) ∗ e − E k t (cid:19) (S4)where N ss (cid:48) n = a (cid:80) m (cid:54) = n (cid:10) n (cid:12)(cid:12) T ( (cid:3) ) (cid:12)(cid:12) m (cid:11) Z sn ( Z s (cid:48) m ) ∗ / (1 − e ( E n − E m ) a ) + (cid:10) m (cid:12)(cid:12) T ( (cid:3) ) (cid:12)(cid:12) n (cid:11) Z sm ( Z s (cid:48) n ) ∗ / (1 − e ( E m − E n ) a ) involvesexcited-state transition matrix elements as well as combinations of overlap factors not determined from fits to Eq. (S2).In order to extract the ground-state matrix elements of interest in this work, combined fits to G ss (cid:48) h ( t ; 0) and R ss (cid:48) h ( t ) □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □□□□□□□□ □□□□ □□□□ □□□□□□ □□□□ □□ □□□□□□□□□□□□ □□ □□□□□□□□ □□□□ □□□□□□ □□ □□ □□□□□□□□ □□□□□□ □□□□□□ □□□□□□ □□□□□□ □□□□ □□□□ □□ □□□□□□ □□□□ □□ □□□□□□□□ □□□□□□ □□□□□□ □□□□□□□□□□ □□ □□□□□□ □□□□□□ □□ □□□□ □□ □□□□ □□□□ □□□□□□ □□ □□□□ □□ □□□□ □□□□ □□□□ □□□□ □□□□ □□ □□□□□□□□ □□ □□□□□□ □□□□□□□□ □□□□□□ □□□□□□ □□ □□ □□□□□□□□ □□ □□□□ □□□□□□ □□ □□ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □□ □□□□□□ □□□□□□ □□□□ □□□□□□□□ □□□□□□□□□□□□□□ □□□□□□□□ □□□□ □□□□□□ □□□□ □□□□ □□□□□□□□ □□□□ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □□□□ □□□□□□ □□□□□□ □□ □□ □□ □□□□□□ FIG. S1. Left: Effective matrix elements (Eq. (6) of the main text), with best fit result (dark grey band) and final fit result(light grey band). Right: The results of all fits included in the weighted average to obtain the systematic and statisticaluncertainties of the corresponding matrix element. are used to fix the free parameters { (cid:10) n (cid:12)(cid:12) T ( (cid:3) ) (cid:12)(cid:12) n (cid:11) , E n , Z sn ( Z s (cid:48) n ) ∗ , N ss (cid:48) n } , for ss (cid:48) ∈ { SS , SP } . It is noteworthy that ifthis spectral representation is truncated at N states = 1, then the second sum in Eq. (S4) vanishes and R ss (cid:48) h ( t ) isindependent of N ss (cid:48) . For N states > R ss (cid:48) h ( t ) is similarly independent from a linear combination of the N ss (cid:48) n thatcan be used to eliminate one redundant parameter, chosen here to be N ss (cid:48) . Further, if the sums and finite-differencederivatives in Eq. (S3) and Eq. (S4) were replaced by continuum integrals and derivatives, then R ss (cid:48) h ( t ) would beindependent of N ss (cid:48) n , suggesting that fits might not be sensitive to N ss (cid:48) n in practice if lattice artifacts are small. Inthis work, fits are performed both using Eq. (S4) and also neglecting these lattice artifacts by setting N ss (cid:48) n = 0 inEq. (S4). The Aikiake Information Criterion (AIC) is used to select whether fits with or without these lattice artifactsare preferred for each choice of N states and fit range that is considered. In all cases, fits without these lattice artifactsare 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 □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ FIG. S2. The ratios of effective matrix elements, Eq. (S5), for pp and He to that in the proton. The gray bands show thereconstructed value of the matrix element ratio from the fits to the individual nuclear states. the model selection (i.e., the number of states to retain in truncated approximations to the infinite sums of states inEqs. (S2)-(S4)). In this work, an automated fitting procedure analogous to the procedure described in Ref. [61] is usedto assess these systematic fitting uncertainties. Results for SS and SP interpolating operators are fit simultaneously andmaximum source-sink separations { t SSmax , t
SPmax } are chosen so that the signal-to-noise ratio of each correlation functionis always greater than a tolerance of 0.5. A minimum separation t minmin = 4 is required for the spectral representationto be well-defined for three-point correlation functions for the action and operator considered here. Within theseconstraints, N max = 200 choices of minimum separations for both sources ( t SSmin , t
SPmin ) are chosen randomly. For eachfit 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 ≤ − .
5, then the two-state fit is accepted. This procedure isrepeated until an N state -state fit is rejected, at which point an ( N state − N state ∈ { , , } . Correlated χ -minimization isused for all fits with covariance matrices regularized using optimal shrinkage [89, 90]. The χ depends linearly on theoverlap factors, so these parameters are solved for using variable-projection methods [91, 92]. The excited-state energygaps 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 toleranceof 10 − ). The fit is then repeated for N boot = 200 bootstrap samples, and the marginalized parameter uncertainty onthe ground-state matrix element is determined using rank-based bootstrap confidence interval estimation [93] to addrobustness to outlier bootstrap samples. Various checks on the fit result are then performed: an uncorrelated fit mustreproduce the result within a tolerance of 5-sigma, the bootstrap median must reproduce the mean within a toleranceof 2-sigma, and the goodness-of-fit must be below a tolerance of χ /N dof <
2. All fit results passing these criteria areincluded in an ensemble of acceptable fits. A weighted average of these acceptable fits is used to determine the finalcentral value and total systematic uncertainty, where the (in principle arbitrary - see Ref. [94] for a discussion of thisin 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 thesystematic uncertainty obtained from the weighted average. The fitting procedure is fully specified by the parametersand tolerance values above, as well as by a random seed for bootstrap resampling that is fixed to allow correlatedratios of matrix elements for different hadrons to be formed.Fig. S1 shows summaries of the fits of R 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 : R h ( t ) = R h ( t ) / R p ( t ) (S5)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 bootstrapresampling to determine ratios of the ground-state matrix elements from all successful pairs of fit ranges, and takinga weighted average of the results. NONPERTURBATIVE RENORMALIZATION
An ensemble of N meas = 101 field configurations is used to compute the nonperturbative renormalization of theoperator studied in this work. The parameters of the ensemble match those used for the nuclear matrix elementcalculations in the main text, except they have a smaller volume of L × T = 24 × T µν in Eq. (2) is computed in the RI (cid:48) MOM scheme [55] byequating the Landau-gauge dressed vertex functionΓ µν ( p ) ≡ N C Tr C (cid:2) S − ( p ) G µν ( p ) S − ( p ) (cid:3) (S6)with its tree level value at a fixed momentum p . Here, S is the quark propagator and G µν is the quark three-pointcorrelation function for the zero-momentum projected operator ˜ T ( (cid:3) ) µν , and the trace is over color degrees of freedom.Defining ˜ T ( (cid:3) ) µν = (cid:80) z T ( (cid:3) ) µν ( z ) = (cid:80) z,z (cid:48) q ( z ) J ( (cid:3) ) µν ( z, z (cid:48) ) q ( z (cid:48) ), this zero-momentum-projected three-point correlationfunction is G µν ( p ) = 1 V (cid:88) x,y,z,z (cid:48) e ip ( x − y ) S ( x, z ) J ( (cid:3) ) µν ( z, z (cid:48) ) S ( z (cid:48) , y ) , (S7)which is computed with the sequential source technique applied through the operator [15]. At tree-level, the vertexfunction is proportional to two tensor structures [56]: i Λ µν ( p ) ≡
12 (˜ p µ γ ν + ˜ p ν γ µ ) − / ˜ pδ µν , i Λ µν ( p ) ≡ ˜ p µ ˜ p ν ˜ p / ˜ p − / ˜ pδ µν , (S8)where ˜ p µ ≡ a sin (cid:0) a p µ (cid:1) is the lattice momentum corresponding to p µ . In the continuum, only Λ µν appears, howeverΛ µν enters as an O ( a ) correction [56]. Expressing Γ µν in the space spanned by { Λ (1) , Λ (2) } amounts to imposing therenormalization condition: Γ µν ( p ) | p = µ = Z q ( p ) (cid:2) Z − ( p )Λ µν ( p ) + Z − ( p )Λ µν ( p ) (cid:3)(cid:12)(cid:12) p = µ , (S9)where Z − a ( p ) for a ∈ { , } are operator renormalization coefficients and the quark field renormalization Z q is definedas: Z q ( p ) | p = µ = i p Tr (cid:2) S − ( p ) / ˜ p (cid:3) (cid:12)(cid:12)(cid:12)(cid:12) p = µ . (S10)The renormalization coefficient of primary interest is Z RI (cid:48) MOM ( µ , a ) ≡ Z ( p ) (cid:12)(cid:12) p = µ .Eq. (S9) is solved by constructing a linear functional on the space of Dirac and Lorentz matrices. A bilinear form (cid:104)· , ·(cid:105) is defined whose action is: (cid:104) λ , λ (cid:105) ≡ (cid:88) i ∈ τ (3)1 Tr D (cid:8) λ i λ i (cid:9) , (S11)where λ , λ are objects with two Dirac and two Lorentz indices and the sum runs over elements of a basis for theirreducible representation of the hypercubic group that contains T , as defined in Eq. (2). The functionals (cid:104) Λ , ·(cid:105) and (cid:104) Λ , ·(cid:105) are applied to Eq. (S9) to yield a system of equations: Z q ( p ) (cid:88) b ∈{ , } (cid:104) Λ a ( p ) , Λ b ( p ) (cid:105)Z − b ( p ) = (cid:104) Λ a ( p ) , Γ( p ) (cid:105) , (S12)which can be solved for Z , ( p ).An MS renormalization factor can be constructed from the RI (cid:48) MOM factor computed as described above as Z MS ( µ ) = R MS / RI (cid:48) MOM ( µ, µ ) Z RI (cid:48) MOM ( µ , a ) , (S13)where the matching factor R MS / RI (cid:48) MOM ( µ, µ ) has been computed to 3-loop order in lattice perturbation theory [56, 57](using the two-loop value produces a statistically indistinguishable result). While Z MS ( µ ) defined in Eq. (S13) is FIG. S3. Renormalization coefficients in the MS scheme, color-coded by the amount of hypercubic breaking at each point. Z MS is computed on each ensemble for each mode p µ with 0 . ≤ ( ap ) ≤
10, and averaged over hypercubic orbits. in principle independent of the matching scale p = µ , in practice there are discretization artifacts which arise ascontamination in the form of dependence on the hypercubic invariants p [2 n ] ≡ (cid:80) µ p nµ , resulting in the classic jellyfish-bone structure shown in Fig. S3. ˜ p µ can be expanded in a Taylor series in p [2 n ] , so either p [2 n ] or ˜ p [2 n ] may be usedto perform the fit; p [2 n ] is used for this analysis, but consistent results are obtained by fitting to ˜ p [2 n ] .To obtain the renormalization factor, the hypercubic artifacts which depend on p [2 n ] for n > p are then fit separately. To allow a quantification of systematic uncertainties, the fits are performed over a range ofwindows, with a range of functional forms. In particular, artifacts of the form { H , H , H } = (cid:40) c a p [4] p + c a p [6] ( p ) , c a p [4] p log( a p ) , c a p [4] + c a p [6] p + c a (cid:18) p [4] p (cid:19) (cid:41) (S14) { R , R , R , R } = (cid:26) d a p , d a p , d a p log( a p ) , d (cid:0) a p (cid:1) (cid:27) , (S15)are considered, where the H i denote hypercubic terms and the R i are contributions from running. Terms are groupedby their order in a and by whether or not they contain logarithmic corrections. The full functional form is truncatedat order a . 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 whilepreventing overfitting using the following procedure. Given a window of momenta, the first fit form is initialized tobe F (1) = 0. Given a fit form F ( n ) , the subsequent form F ( n +1) is determined by considering all possible forms F ( n +1) j = F ( n ) + X j , (S16)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 ( n +1) j is accepted if and only if A ( F ( n +1) j ) < A ( F nj ), where A ( F ) ≡ N param ( F ) + χ ( F ) is the value of the AIC of the fit F . If no forms are acceptedor there are no fit forms left, iteration stops and F ≡ F ( n ) . The subsequent fit form F ( n +1) is chosen out of theaccepted fit forms { F ( n +1) j } j by maximizing the p -value of the fit.This procedure is applied to a range of fitting windows. Windows are chosen with ( ap ) ∈ [( ap ) , ( ap ) ] and p [4] / ( p ) ≤ h by independently adjusting ( ap ) and ( ap ) between 0.5 and 10 in increments of 0.5, alwayskeeping a window size ( ap ) − ( ap ) ≥
5, and taking h ∈ { . , . , ..., . } . For each fit window f , the optimalfit form as defined above is accepted if and only if its p value p f ≥ .
01, and each accepted fit is given a weight: w f ≡ p f ( δ Z MS f ) − , (S17)where δ Z MS f denotes the statistical uncertainty on Z MS f , which is the fit result for Z MS that is obtained in fit f . Theweights are normalized so that the maximum weight is 1. Given the large number of accepted fits, the N w ≡ FIG. S4. N w = 100 highest-weight accepted fits, Z f , plotted as a function of their fit weight, w f . Individual fits are color-codedby the maximum amount of hypercubic breaking in the window, and the gray band shows the final weighted average. highest-weight fits are combined in a weighted average, with the result: Z MS ( µ = 2 GeV) = 0 . . (S18)The final result and uncertainties are consistent under variation of N w , with larger systematic uncertainties whenadditional fits with lower weights are included (e.g., with N w = 200, the result is 0.895(55)). As a consistency check,the computation was also performed on an ensemble with identical action and parameters with lattice dimensions16 ×
48. A result of Z MS ( µ = 2 GeV) = 0 . CONNECTION TO PHENOMENOLOGY: MOMENTUM FRACTION EXTRAPOLATION
In order to connect to phenomenology, the lattice result for (cid:104) x (cid:105) ( He) u − d at quark masses corresponding to m π = 806MeV is extrapolated to the physical masses using the assumption of weak mass dependence of the short-distancetwo-body counterterms in nuclear EFT. The isovector twist-two operators in Eq. (1) match on to hadronic operatorsin nuclear EFT as [63] O µ ...µ n → (cid:104) x n (cid:105) ( p ) u − d v µ . . . v µ n N † τ N (1 + α ,n N † N ) + . . . (S19)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 ( N, Z ) = (cid:104)
N, Z | N † τ N N † N | N, Z (cid:105) , the two-nucleon counterterm α , relates the nuclearand proton momentum fractions as [63] α , G ( N, Z ) ≡ (cid:104) x (cid:105) ( N,Z ) u − d − Z − NA (cid:104) x (cid:105) ( p ) u − d . (S20)This LEC can be determined from the LQCD results for He most precisely by re-expressing it in terms of thequantities in Table I as α , G ( He) = 13 (cid:104) x (cid:105) ( He) u − d (cid:104) x (cid:105) ( p ) u − d − (cid:104) x (cid:105) ( p ) u − d . (S21)Since renormalization effects cancel in the ratio, it is more precisely determined than the individual momentumfractions themselves and, with a naive error propagation, computing α , G ( 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:1. Determine the counterterm α , G ( He) at m π = 806 MeV via Eq. (S21), using the LQCD calculations of theratio of (cid:104) x (cid:105) ( He) u − d / (cid:104) x (cid:105) ( p ) u − d , and (cid:104) x (cid:105) ( p ) u − d itself.2. While the momentum fractions themselves have nonanalytic quark mass dependence, the counterterms α , andnuclear factors G ( N, Z ) are expected to have only mild quark mass dependence and so their their combinationis extrapolated to the physical quark masses by assuming the same central value and increasing the uncertaintyby 50%.3. The extrapolated value of α , G ( He) is combined with the value of (cid:104) x (cid:105) ( p ) u − d from phenomenology to produce aphysical-point value of (cid:104) x (cid:105) ( He) u − dd
N, Z | N † τ N N † N | N, Z (cid:105) , the two-nucleon counterterm α , relates the nuclearand proton momentum fractions as [63] α , G ( N, Z ) ≡ (cid:104) x (cid:105) ( N,Z ) u − d − Z − NA (cid:104) x (cid:105) ( p ) u − d . (S20)This LEC can be determined from the LQCD results for He most precisely by re-expressing it in terms of thequantities in Table I as α , G ( He) = 13 (cid:104) x (cid:105) ( He) u − d (cid:104) x (cid:105) ( p ) u − d − (cid:104) x (cid:105) ( p ) u − d . (S21)Since renormalization effects cancel in the ratio, it is more precisely determined than the individual momentumfractions themselves and, with a naive error propagation, computing α , G ( 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:1. Determine the counterterm α , G ( He) at m π = 806 MeV via Eq. (S21), using the LQCD calculations of theratio of (cid:104) x (cid:105) ( He) u − d / (cid:104) x (cid:105) ( p ) u − d , and (cid:104) x (cid:105) ( p ) u − d itself.2. While the momentum fractions themselves have nonanalytic quark mass dependence, the counterterms α , andnuclear factors G ( N, Z ) are expected to have only mild quark mass dependence and so their their combinationis extrapolated to the physical quark masses by assuming the same central value and increasing the uncertaintyby 50%.3. The extrapolated value of α , G ( He) is combined with the value of (cid:104) x (cid:105) ( p ) u − d from phenomenology to produce aphysical-point value of (cid:104) x (cid:105) ( He) u − dd / (cid:104) x (cid:105) ( p ) u − dd