# The axial charge of the triton from lattice QCD

Assumpta Parreño, Phiala E. Shanahan, Michael L. Wagman, Frank Winter, Emmanuel Chang, William Detmold, Marc Illa

IICCUB-21-001, FERMILAB-PUB-21-026-T, MIT-CTP/5274

The axial charge of the triton from lattice QCD

Assumpta Parre˜no, Phiala E. Shanahan,

2, 3

Michael L. Wagman,

2, 4

Frank Winter, Emmanuel Chang, William Detmold,

2, 3 and Marc Illa (NPLQCD collaboration) 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 Center for Theoretical Physics, MassachusettsInstitute of Technology, Cambridge, MA 02139, USA The NSF AI Institute for Artiﬁcial Intelligence and Fundamental Interactions Fermi National Accelerator Laboratory, Batavia, IL 60510, USA Jeﬀerson Laboratory, 12000 Jeﬀerson Avenue, Newport News, VA 23606, USA (Dated: February 9, 2021)The axial charge of the triton is investigated using lattice quantum chromody-namics (QCD). Extending previous work at heavier quark masses, calculations areperformed using three ensembles of gauge ﬁeld conﬁgurations generated with quarkmasses corresponding to a pion mass of 450 MeV. Finite-volume energy levels for thetriton, as well as for the deuteron and diproton systems, are extracted from analysisof correlation functions computed on these ensembles, and the corresponding ener-gies are extrapolated to inﬁnite volume using ﬁnite-volume pionless eﬀective ﬁeldtheory (FVEFT). It is found with high likelihood that there is a compact boundstate with the quantum numbers of the triton at these quark masses. The axial cur-rent matrix elements are computed using background ﬁeld techniques on one of theensembles and FVEFT is again used to determine the axial charge of the proton andtriton. A simple quark mass extrapolation of these results and earlier calculationsat heavier quark masses leads to a value of the ratio of the triton to proton axialcharges at the physical quark masses of g H A /g pA = 0 . +0 . − . . This result is consistentwith the ratio determined from experiment and prefers values less than unity (inwhich case the triton axial charge would be unmodiﬁed from that of the proton),thereby demonstrating that QCD can explain the modiﬁcation of the axial charge ofthe triton. a r X i v : . [ h e p - l a t ] F e b I. INTRODUCTION

The triton ( H) is the simplest nucleus that undergoes weak decay, and as such is animportant system with which to investigate the Standard Model (SM) origins of nuclearphysics. While at the quark level the weak decay of the triton is mediated by the weakinteractions, it is the strong interactions, described by Quantum Chromodynamics (QCD),that dictate the embedding of quarks inside nuclei and thus are key to the nuclear decayrate. A notable feature of the β -decay of the triton and other nuclei is a reduction of theGamow-Teller (isovector axial-vector) transition rate as compared to that of the neutron;this reduction scales approximately with the atomic number of the system, A , for mediummass nuclei [1–4] and can be described phenomenologically by an ad hoc reduction of thein-nucleus axial charge of the proton or by the introduction of two-nucleon interactionswith the weak current [5]. For the triton, this manifests as a suppression of the isovectoraxial matrix element; analysis of precision tritium decay experiments ﬁnds that the ratioof the axial charge of the triton to that of the proton is g H A /g pA = 0 . β -decay is a very promising avenue through which to improveconstraints on neutrino masses [7–10]. Careful comparison of decay measurements withtheoretical predictions may also reveal physics beyond that of the SM [11, 12]. Moreover,nuclear eﬀects in Gamow-Teller transitions are important for understanding neutrinoful andneutrinoless double- β decay matrix elements and thereby for sensitive searches for leptonnumber violation [13].The half-life of tritium, t / , is related to the Fermi ( (cid:104) F (cid:105) ) and Gamow-Teller ( (cid:104) GT (cid:105) )reduced matrix elements through t / = K/G V (1 + δ R ) f V (cid:104) F (cid:105) + f A /f V g A (cid:104) GT (cid:105) , (1)where the constants K , G V , and δ R are known precisely from theory or experiment [14], f A,V denote known Fermi functions [14], and g A ≡ g pA is the axial charge of the nucleon. Fromthe Ademollo-Gatto theorem [15], (cid:104) F (cid:105) ∼ A ak = ¯ qγ k γ τ a q, (2)where q = ( u, d ) T denotes a doublet of quark ﬁelds, γ µ are Dirac matrices, and τ a are Paulimatrices in ﬂavor space. In particular, the Gamow-Teller matrix element is deﬁned from thezero-momentum projected current ˜ A ak as (cid:68) He , s (cid:12)(cid:12)(cid:12) ˜ A + k (cid:12)(cid:12)(cid:12) H , s (cid:48) (cid:69) = ¯ U s γ k γ τ + U s (cid:48) g A (cid:104) GT (cid:105) , (3)where g H A ≡ g A (cid:104) GT (cid:105) , U s is a relativistic spinor for the nucleus spin component s ∈ { + , −} ,and τ + = τ + iτ . Eq. (3) is valid for zero electron mass and vanishing electron momentum.Because of the low energy scale of the β -decay process, determinations of the axial chargesof hadronic and nuclear systems probe QCD in the non-perturbative regime; theoretical Label

L/a T /a L [fm] T [fm] m π L m π T N cfg N meas E24 24 64 2.80 7.47 6.390 17.04 2124 5 . × E32 32 96 3.73 11.2 8.514 25.54 2850 2 . × E48 48 96 5.60 11.2 12.78 25.49 980 4 . × TABLE I. Details of the ensembles of gauge ﬁeld conﬁgurations used in this calculation, and ofthe measurements performed on them. In all cases, the gauge coupling is β = 6 .

1, the tadpole-improved Sheikholeslami-Wohlert parameter is c SW = 1 . am (bare) ud = − . am (bare) s = − . N cfg conﬁgurations on which a total of N meas measurements were made (slightlydiﬀerent numbers of sources were used on each conﬁguration). determinations must therefore be undertaken using lattice QCD (LQCD) which is the onlyknown systematically improvable tool for calculations in this regime. The axial charge of theproton has been calculated using LQCD for many years following the ﬁrst studies in Ref. [16],see Ref. [17] for a recent summary of results. A ﬁrst calculation of the axial charge of thetriton was presented in Ref. [18], albeit at an SU (3)-symmetric point with quark massescorresponding to a pion mass m π = 806 MeV. This work extends Ref. [18] with calculationsat quark masses that are considerably closer to their physical values, corresponding to m π =450 MeV and a kaon mass m K = 595 MeV [19]. At these quark masses, the inﬁnite volumeextrapolated Gamow-Teller matrix element is determined to be (cid:104) GT (cid:105) L →∞ = 0 . (cid:104) GT (cid:105) = g H A /g pA = 0 . +0 . − . . These are the main results of this work and showthat the phenomenological modiﬁcation of the axial charge of the triton can be reproduceddirectly from QCD. II. LATTICE QCD DETAILS

The lattice QCD calculations presented in this work make use of isotropic gluon con-ﬁgurations generated with a L¨uscher-Weisz [20] gauge action and N f = 2 + 1 ﬂavors offermions implemented using a tadpole-improved clover quark action [21]. All ensembles aregenerated using a gauge coupling of β = 6 . m π = 450 MeV and a strange quark mass that correspondsto a kaon mass of m K = 595 MeV. Performing scale setting using upsilon spectroscopyextrapolated to the physical quark masses results in a lattice spacing of a = 0 . III. SPECTROSCOPY AND INFINITE VOLUME EXTRAPOLATION AT m π ∼ MEV

In order to determine the ground-state energy of the triton and He, which are degeneratefor the isospin-symmetric quark masses used in this calculation, two-point correlation func-tions projected to zero three-momentum are constructed using the methodology introducedin Refs. [23, 24]. The correlation functions are C h,s ( t ) = (cid:88) x Γ ( s ) βα (cid:68) χ ( h ) (cid:48) α ( x , t ) χ ( h ) † β ( , (cid:69) for s ∈ { + , −} , (4)where Γ ( ± ) = (1 + γ )(1 ± iγ γ ) is a projector onto the given positive energy spin compo-nent and χ ( h ) (cid:48) α and χ ( h ) α are interpolating operators with the quantum numbers of the hadron h ∈ { p, H } . For the triton, the interpolating operator is built from three color-singlet nu-cleons that are independently projected to zero three-momentum in Eq. (4). The quarkpropagators used to construct the correlation functions are computed using APE smeared[25] sources and point or APE smeared sinks; the resulting correlation functions are referredto as ‘SP’ and ‘SS’ respectively. An advantage of using local multi-hadron sources in thismanner is that they can be eﬃciently computed for light nuclei such as the triton usingbaryon block algorithms whose cost scales linearly with the spatial volume [23]. A disadvan-tage of this approach is that since source and sink operators are distinct, it is not possibleto build a variational basis of operators that would explicitly account for excited-state con-tamination from unbound multi-nucleon states that have small energy separations to theground state for large volumes. Similar issues arise in the two-nucleon sector, and Ref. [22]discusses a number of consistency checks that have been applied to two-nucleon results us-ing multiple interpolating operators on the same ensemble as used here. Given tensionsbetween two-nucleon ground-state energy results with m π ∼

800 MeV obtained using prod-ucts of zero-momentum baryon sources [26, 27] compared with results obtained using localor approximately local two-baryon sources [24, 28], it will be important to pursue variationalstudies of both the two-nucleon and three-nucleon sectors, including operators overlappingwith bound and unbound states, in the future. However, multi-nucleon variational studiesrequire a large set of interpolating operators whose correlation functions are signiﬁcantlymore costly to compute than those calculated here and are deferred to future work. Whilehere the focus is on the triton and the nucleon, other single hadron and two-baryon systemshave been studied using the same approach as applied here, on the same ensembles of gaugeﬁeld conﬁgurations, as discussed in Refs. [19, 22].The ground-state triton energy and its diﬀerence from the mass of three non-interactingnucleons are extracted in each volume from analysis of the time dependence of C H , ± ( t ) and C p, ± ( t ) using the same ﬁtting methodology as applied and detailed in Refs. [22, 29]. Forcompleteness, the approach is summarized here. Provided that the temporal separation ofthe source and sink, t , is larger than the extent of the lattice action, and small comparedwith the temporal extent of the lattice geometry ( t (cid:28) T ), the correlation functions are givenby a sum of exponentials whose exponents are determined by the energies of the eigenstatesof the given quantum numbers: C SS (cid:48) h, ± ( t ) = N ex (cid:88) n =0 Z ( h ) S n Z ( h ) S (cid:48) ∗ n exp( − E ( h ) n t ) for {S , S (cid:48) } ∈ { S , P } , (5)where N ex excited states contribute to the sum, Z ( h )S/P n denotes the overlap factor of thecorresponding interpolating operators onto the n th energy eigenstate, and thermal eﬀects For the L = { , , } ensembles, smearing parameters of N smear = { , , } steps of smearing withGaussian smearing parameters ρ = { . , . , . } were used. arising from the ﬁnite temporal extent of the lattice geometry have been neglected. Fitsof the correlation functions to Eq. (5) determine the ground state energies M h = E ( h )0 among the ﬁt parameters; while correlation functions for diﬀerent propagator smearingshave diﬀerent overlap factors, the energies are common and are thus ﬁt simultaneously. Toquantify the systematic uncertainties arising from the choice of source-sink separations, t ,included in the ﬁt, and of the truncation of the sum in Eq. (5), 200 ﬁt windows are sampledat random from the set of all choices of contiguous time-ranges longer than t plat = 6 and withﬁnal times less than t max (deﬁned by the point at which the noise-to-signal ratio exceeds tol noise = 0 .

25 for the given correlation function). For each ﬁt range, the Akaike InformationCriterion [30] (AIC) is used to perform model selection ( i.e. , ﬁx the number of exponentialterms in the sum above). The truncation is set as the largest N ex for which the change inAIC is ∆AIC < − . N dof , where N dof is the number of degrees of freedom of the ﬁt. Ineach case, combined correlated ﬁts are performed to average correlation functions as well asto N boot = 200 bootstrap resamplings of the correlation functions using covariance matricescomputed using optimal shrinkage [31–33] and using variable projection (VarPro) [34, 35]to eliminate overlap factors. All ﬁts with a χ /N dof < tol χ = 2 are included in a set of‘accepted’ ﬁts (accepted ﬁts must also pass tests that the ﬁt results are a) invariant to thechoice of the minimizer that is used to within tol solver = 10 − , b) agree within a prescribedtolerance of tol corr = 5 σ with uncorrelated ﬁts, and c) agree within a prescribed toleranceof tol median = 2 σ to the median bootstrap result, as in Refs. [22, 29]). The ﬁnal value anduncertainties for the energy are then computed from the results of all N success accepted ﬁtsusing a weighted average. The central value and statistical and systematic uncertainties arecomputed as E = N success (cid:88) f =1 w f E f , δ stat E = δE f max , δ sys E = N success (cid:88) f =1 w f (cid:16) E f − E (cid:17) , (6)where E f denotes the ﬁt result from a given ﬁt labelled by f , and w f is the associatedweight. For each ﬁt range, the statistical uncertainties δE f are computed using 67% conﬁ-dence intervals determined from the results of the N boot resampled correlation function ﬁtsdescribed above; the total statistical uncertainty is deﬁned as the statistical uncertainty ofthe ﬁt f max with maximum weight w f max . The statistical and systematic uncertainties arecombined in quadrature to give a total uncertainty δE = (cid:113) δ stat E + δ sys E . Since eachsuccessful ﬁt provides a statistically unbiased estimate of the ground-state energy, the rela-tive weights w f of each ﬁt in the weighted average can be chosen arbitrarily (in the limit oflarge statistics) [33]. Here, as in Refs. [22, 29, 33], the weights are set as w f = p f (cid:16) δE f (cid:17) − (cid:80) N success f (cid:48) =1 p f (cid:48) (cid:16) δE f (cid:48) (cid:17) − , (7)where p f = Γ( N dof / , χ f / / Γ( N dof /

2) is the p -value of ﬁt f assuming χ -distributedgoodness-of-ﬁt parameters.The resulting ﬁtted masses are summarised for each volume in Table II and the ﬁts areshown graphically for the proton in Fig. 1 and for the triton in Fig. 2. In each ﬁgure, the □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ FIG. 1. Left column: eﬀective mass functions (Eq. (8)) for the proton correlation function on theE24 (top), E32 (center), and E48 (bottom) ensembles. The orange circles and blue squares showeﬀective masses constructed from SP and SS correlation functions respectively, and the correspond-ing orange and blue curves display the highest-weight ﬁt contributing to the weighted average ofEq. (6); the light (dark) grey bands show the mass extracted from the weighted average (singlehighest-weight ﬁt). Right column: masses extracted from each successful ﬁt to the correlators.The opacity of each point is set by the contribution of the ﬁt to the weighted average. The hor-izontal line and grey band in each ﬁgure correspond to the ﬁnal result for the mass and its totaluncertainty. eﬀective mass functions aM SS (cid:48) h ( t ) = ln (cid:32) (cid:80) s = ± C SS (cid:48) h,s ( t ) (cid:80) s = ± C SS (cid:48) h,s ( t + a ) (cid:33) , (8)are shown for the SS and SP data along with the functional form of best ﬁts to the correlationfunctions and the ﬁnal result for E ( h )0 . The diﬀerence between the triton mass and three □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ FIG. 2. Eﬀective mass functions and ﬁt summaries for the ﬁts of the triton correlation functionon the E24 (top), E32 (center), and E48 (bottom) ensembles. All features are as in Fig. 1.Ensemble aE ( p )0 aE ( H)0 a ∆ E H E24 0 . . − . . . − . . . − . E H = E ( H)0 − E ( p )0 , determined in lattice units on each ensemble. times the proton mass, ∆ E H = E ( H)0 − E ( p )0 , is determined from the correlated diﬀerencesof the ﬁtted energies computed under bootstrap. The results for this quantity are alsopresented in Table II.In the limit in which m π L is suﬃciently large and the pion mass is suﬃciently smallthat the volume dependence of the proton mass is described by p -regime chiral perturbationtheory ( χ PT), that dependence takes the form [36, 37] M p ( L ) − M ∞ p = − g A π f π K (0 , L ) − g N π f π K (∆ , L ) , (9)where f π is the chiral limit pion decay constant and K (∆ , L ) ≡ (cid:90) ∞ dλβ ∆ (cid:88) (cid:126)n (cid:54) = (cid:126) (cid:2) ( L | (cid:126)n | ) − K ( Lβ ∆ | (cid:126)n | ) − β ∆ K ( Lβ ∆ | (cid:126)n | ) (cid:3) . (10)Here, L denotes the spatial extent of the lattice geometry, K n are modiﬁed Bessel functions, β ≡ λ +2 λ ∆+ m π , ∆ denotes the mass splitting between the nucleon and the ∆ baryon, and g ∆ N is the nucleon-∆ transition axial coupling. The sum in Eq. (10) is over integer tripletsexcluding (cid:126)n = (0 , , ∼ e − m π L /L .While m π = 450 MeV is not in the regime where chiral perturbation theory shows rapidconvergence for baryon properties, it has been found to eﬀectively describe the volume-dependence of baryon masses in this mass regime [38]. With the physical values of g A = 1 . g ∆ N = 1 . f π = 132 MeV and ∆ = 300 MeV, Eq. (9) is used to ﬁt to the inﬁnite volumeproton mass from the masses determined on the three volumes. This ﬁt, displayed in Figure3, results in an inﬁnite volume mass of M ∞ p = 1 . Figure 3 also shows the diﬀerence between the triton energy and the three-nucleon thresh-old for each of the three volumes. Unlike for the proton, the form of the volume dependenceof the triton energy is not known analytically, however a numerical description can be foundby matching to ﬁnite-volume eﬀective ﬁeld theory (FVEFT). The procedure of matching pi-onless EFT to LQCD results for the binding energies of light nuclei using the stochastic vari-ational method has been demonstrated in Ref. [39] using LQCD results with m π = 806 MeV.The same procedure as detailed in that work is followed for the data presented here, andfurther details of the FVEFT approach used here are provided in Appendix A. The in-ﬁnite volume extrapolation leads to an energy shift ∆ E ∞ H = − . .

2) and 11 . .

6) MeV from the FVEFT extrapolation respectively (illustrated in Fig. 9in Appendix A; note also that these results are consistent, within 1 standard deviation, withthose obtained from this data via L¨uscher’s method in Ref. [19, 22]). While the latter casesof 2+1 body systems can not be ruled out, there is a strong preference (80% likelihood, usingthe most conservative two-body binding energies) that the state is a compact three-bodysystem. In what follows, this is assumed to be the case and the state is referred to as the‘triton’.In Fig. 4, the resulting binding energy of the triton is compared to the results of othercalculations, including that of Ref. [24] using the same action but heavier quark massescorresponding to m π = 806 MeV. The extracted binding energy is compatible with those This is in agreement with the value M ∞ p = 1 . As well as the EFT approach, recent generalizations of the quantization conditions derived by L¨uscherfor two particles relate ﬁnite-volume three-particle energies to the inﬁnite-volume three-particle scatteringamplitude and constrain bound states when present, see Ref. [40] for a recent review. □ □ □ □□ (a) □□ □□□ □ □◦ ◦ ◦ ◦ ◦ ◦◦◦ □□◦ ◦ ◦ ◦ ◦ ◦ ◦◦ ◦ ◦ ◦ ◦ ◦ ◦ - - - (b) FIG. 3. (a) Mass of the proton as a function of the spatial extent of the lattice geometry, alongwith the ﬁt to this dependence via Eq. (9). (b) The binding energy of the triton as a functionof the spatial extent of the lattice geometry. Grey squares correspond to the binding energiesdetermined in the LQCD calculations, while the blue curve shows the ﬁnite-volume dependence inpionless EFT ﬁtted to the LQCD data. The blue point at the right indicates the inﬁnite volumeextrapolation of the binding energy. For comparison, the binding energies of the pp (purple) anddeuteron (green) systems obtained via an analogous FVEFT study are also shown. of other calculations at nearby quark masses [41, 42], although no eﬀort is made to takeinto account the diﬀerences between the lattice actions or scale-setting schemes that areused. Naive extrapolations of the current result and that from Ref. [24] that are linear orquadratic in m π are consistent with the experimental value for the binding energy of thetriton, albeit with large uncertainties. The strong evidence for binding at the other massesshown in the ﬁgure, and the assumption of smooth behavior under variation in the pionmass provides additional support for the conclusion that the triton is a compact three-bodysystem at these quark masses.0 □ □□ □ ** - - FIG. 4. Binding energy of the triton as a function of the pion mass. The current result isshown in blue, along with the result obtained in a previous calculation [24] using the same actionat m π = 806 MeV. The results of Refs. [41, 42] are shown in orange. The physical value of thebinding is shown as a red star at the physical pion mass, as indicated by the vertical line. Linearand quadratic ﬁts in m π to the blue data are shown as the green and blue bands respectively. IV. GAMOW-TELLER MATRIX ELEMENT FOR TRITIUM β -DECAY To extract the axial charge of the triton, and hence the Gamow-Teller contribution totritium β -decay, LQCD calculations of matrix elements of the isovector axial current in thetriton and proton are performed using the E32 ensemble. Only a single ensemble is useddue to computational cost and for technical simplicity, the ﬂavor-diagonal matrix element (cid:104) H | ˜ A | H (cid:105) is studied and is related by isospin symmetry to the tritium β -decay matrixelement. The resulting values are matched to FVEFT to determine the relevant low energyconstants (LECs), which are then used to predict the inﬁnite-volume matrix element.To compute the ﬁnite-volume matrix elements, the extended propagator technique dis-cussed in detail in Refs. [18, 43–46] is used. This requires extraction of hadronic and nuclearcorrelation functions at a range of values of an applied constant axial ﬁeld that couples toup and down quarks separately. Extended quark propagators are deﬁned as S ( f ) λ f ( x, y ) = S ( f ) ( x, y ) + λ f (cid:90) d zS ( f ) ( x, z ) γ γ S ( f ) ( z, y ) , (11)where S ( f ) ( x, y ) is a quark propagator of ﬂavor f ∈ { u, d } and λ f is the strength of theapplied ﬁeld for the given ﬂavor. These quantities are calculated for ﬁve values of the exter-nal ﬁeld strength λ u,d ∈ { , ± . , ± . } in lattice units. Two-point correlation functions C ( λ u ,λ d ) h,s ( t ) are constructed from these extended propagators using the same contractionmethods as for the zero-ﬁeld correlation functions discussed in the previous section. Forclarity, the smearing labels SS/SP are suppressed in this section. The two-point correlationfunctions are polynomials in λ f of order n f , the number of quark propagators of ﬂavor f inthe correlation function. With computations at n f diﬀerent ﬁeld strength values, the termsin the polynomial can thus be extracted uniquely [44] and are labelled as C ( λ f , h,s (cid:12)(cid:12)(cid:12) O ( λ kf ) and1 C (0 ,λ f ) h,s (cid:12)(cid:12)(cid:12) O ( λ kf ) for k ∈ { , . . . , n f } . It is straightforward to show [18, 44] through insertions ofcomplete sets of states that an eﬀective isovector axial charge function, which asymptotesto the desired bare axial charge in the corresponding hadronic or nuclear state at largetemporal separations, can be deﬁned as g hA ( t ) = R h ( t + a ) − R h ( t ) −→ t →∞ g h (0) A + O ( e − δt ) , (12)where the (0) superscript denotes that the charge is not renormalized, δ is the energy dif-ference between the ground state and the ﬁrst excited state with the quantum numbers ofthe state h , and R h ( t ) = (cid:88) s = ± η s C ( λ u , h,s ( t ) (cid:12)(cid:12)(cid:12) O ( λ u ) − C (0 ,λ d ) h,s ( t ) (cid:12)(cid:12)(cid:12) O ( λ d ) a C (0 , h,s ( t ) (13)for η ± = ∓ O ( λ ) correlation functionsappearing in Eq. (13) can be constructed as C ( λ u , h,s ( t ) (cid:12)(cid:12)(cid:12) O ( λ u ) − C (0 ,λ d ) h,s ( t ) (cid:12)(cid:12)(cid:12) O ( λ d ) = aη s t/a (cid:88) τ/a =0 (cid:88) n,m Z ( h ) n Z ( h ) (cid:48)∗ m e − E ( h ) n ( t − τ ) e − E ( h ) m τ g h ( m,n )(0) A , (14)where the prime on the second Z -factor is included to denote that, although smearings aresuppressed in this section, the overlap factors diﬀer at the source and sink, and the bare(transition) charge g h ( m,n )(0) A is deﬁned from the corresponding matrix element as (cid:68) h, n, s (cid:12)(cid:12)(cid:12) ˜ A (cid:12)(cid:12)(cid:12) h, m, s (cid:48) (cid:69) = ¯ U n,s γ γ τ U m,s (cid:48) g h ( m,n )(0) A , (15)where | h, n, s (cid:105) denotes states with the quantum numbers of h and U m,s denotes the spinorfor state m with spin s . This can be used to derive a spectral representation for g hA ( t ): g hA ( t ) = (cid:88) n Z ( h ) n Z ( h ) (cid:48)∗ n g h ( n,n )(0) A (cid:32) ( t/a + 1) e − E ( h ) n ( t + a ) (cid:80) k Z ( h ) k Z ( h ) (cid:48)∗ k e − E ( h ) k ( t + a ) − ( t/a ) e − E ( h ) n t (cid:80) k Z ( h ) k Z ( h ) (cid:48)∗ k e − E ( h ) k t (cid:33) + O ( a ) , (16)where the O ( a ) contribution is detailed in Ref. [47] and depends on excited-state transitionmatrix elements as well as combinations of overlap factors not determined from ﬁts of two-point functions to Eq. (5). Notably, the O ( a ) term is absent for a single-state correlation The correlation function is deﬁned in Euclidean spacetime, and the sum over τ extends only over thetemporal range between the source and the sink because of the scalar isoscalar nature of the vacuum(exponentially small contributions that are suppressed by the mass of the lightest axial-vector meson areignored). □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □□ □□ □□ □□ □□ □□ □□ □□ □□ □□ □□ FIG. 5. The left panels show the proton and triton bare eﬀective matrix element functions, g pA ( t )and g H A ( t ), respectively, calculated on the E32 ensemble. In each ﬁgure, the orange (blue) datacorrespond to the SP (SS) calculations and the corresponding shaded bands illustrate the highestweight ﬁt. The light (dark) shaded gray bands denote the extracted values of the matrix elementsarising from the combined analysis (single highest-weight ﬁt). The right panels show the valuesobtained for all the successful ﬁts to diﬀerent time ranges, with the opacity determined by thecontribution of the ﬁt in the ﬁnal weighted average. The horizontal line and gray band in the rightpanels show the ﬁnal central value and uncertainty. function model. Multi-state ﬁts have been performed both with and without these O ( a )terms and the AIC is used to determine whether the O ( a ) terms should be included in theﬁt. For both the triton and proton, this AIC test prefers the ﬁt without O ( a ) terms for all ﬁtrange choices. Combined ﬁts of (cid:80) s = ± C h,s ( t ) to Eq. (5) and g hA ( t ) to Eq. (16) without O ( a )terms using both SP/SS interpolating operator combinations are therefore used to determine E ( h ) n , the product Z ( h ) n Z ( h ) (cid:48)∗ n , and g h ( n,n ) A . For each ﬁtting interval, ln( E ( h ) n ) and g h (0) A = g h ( n,n )(0) A are used as nonlinear optimizer parameters with Z ( h ) n obtained from C h,s ( t ) usingVarPro as above and g h ( n,n )(0) A subsequently obtained from g hA ( t ) using VarPro. Statisticaluncertainties on the ground-state matrix elements for each ﬁt are obtained using bootstrapconﬁdence intervals, and a weighted average performed analogously to the two-point functioncase described above is used to determine the ﬁnal ground-state matrix element values andstatistical plus ﬁtting systematic uncertainties. Results for the ground-state matrix elements g h (0) A obtained using this ﬁtting procedure for both one- and three-nucleon systems are shownin Table III.The quantities g pA ( t ) and g H A ( t ) constructed on ensemble E32 for both SS and SP source–3 □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ FIG. 6. The eﬀective ratio function g H A ( t ) /g pA ( t ), constructed through Eq. (12), which asymptotesto g H A /g pA and is independent of axial current renormalization. The orange (blue) data show theSP (SS) correlation functions on the E32 ensemble and the shaded band indicates the extractedvalue of the matrix element (not a ﬁt to the displayed data).h g h (0) A g hA g hA /g pA p . . H 1 . . . g h (0) A , and renormalized, g hA , axial charges of the proton andtriton as well as the corresponding ratio to the axial charge of the proton. The statistical uncer-tainties and ﬁtting systematic uncertainties are combined in quadrature. The ﬁtting systematicuncertainties account for variation in the choice of ﬁt ranges and multi-state ﬁt models used toobtain g h (0) A and the resulting correlated ratio of ground-state matrix elements g hA /g pA as describedin the main text. sink pairs are shown in Fig. 5. Also shown are the values of the bare matrix elementsdetermined from ﬁts to the time-dependence of these functions as discussed above. Ta-ble III displays the extracted bare couplings g h (0) A as well as renormalized couplings g hA = Z A g h (0) A obtained by multiplying by the appropriate axial current renormalization constant, Z A = 0 . g H A ( t ) /g pA ( t ), which at large timesasymptotes to the GT reduced matrix element (cid:104) GT (cid:105) = g H A /g pA , and is independent of therenormalization of the axial current, is shown in Fig. 6 along with the value of (cid:104) GT (cid:105) thatis extracted from the ﬁts to the individual axial couplings.The FV three-body matrix element in Table III can be used to constrain the leadingtwo-nucleon axial current counterterm of pionless eﬀective ﬁeld theory in the ﬁnite volume.To do so, the approach developed in Ref. [49] is followed, whereby EFT wavefunctions,determined variationally and matched to the LQCD spectrum computed in the E32 volume,are used to evaluate the FVEFT matrix elements of the EFT current: A i,a = g A N † τ a σ i N + L ,A (cid:104)(cid:0) N T P i N (cid:1) † (cid:0) N T ¯ P a N (cid:1) + h.c. (cid:105) + . . . , (17)4 - - - (a) ◦ ◦ ◦ ◦ ◦ ◦□□ □□ (b) FIG. 7. (a) Determination of the LEC ratio (cid:101) L ,A by matching the matrix element deﬁned inEq. (19), computed in the optimized FVEFT wavefunction for the E32 volume, to the LQCDresult in the E32 volume, which is depicted as the horizontal band. (b) The FVEFT extrapolationof the ratio of triton to proton axial charges to inﬁnite volume. where g A is the single-nucleon axial coupling, and the projectors P i ≡ √ σ σ i τ , P a ≡ √ σ τ τ a , (18)form spin-triplet and spin-singlet two-nucleon states, respectively. The two-body couterterm L ,A is regulator-dependent, and in this work a Gaussian regulator scheme is used with ascale Λ as discussed in Appendix A.The ratio of the triton to proton matrix elements in FVEFT is given by2 g A (cid:104) Ψ H , s = + |A , | Ψ H , s = + (cid:105)(cid:104) Ψ H , s = + | Ψ H , s = + (cid:105) = (cid:18) L ,A g A h H (Λ , L ) (cid:19) , (19)where | Ψ H , s (cid:105) is the wavefunction for the spin- s H state and h H (Λ , L ) is the spatial ex-pectation value of a regulated form of A , in the triton spatial wavefunction as detailed inAppendix A. The ratio of LECs (cid:101) L ,A ≡ L ,A /g A is determined by demanding that Eq. (19)for the E32 volume reproduces the LQCD ratio of axial charges for H and for the protoncomputed on that volume. This value of (cid:101) L ,A is then used to compute the axial current ma-trix element in variationally-optimised triton wavefunctions for diﬀerent volumes, includingthe inﬁnite volume wave function that allows the inﬁnite-volume matrix element to be deter-mined (more details on this procedure are given in Appendix A). While the counterterm L ,A is scheme and scale-dependent, the triton axial charge for any volume is scale-independent.Figure 7 shows the result of this matching procedure and the volume dependence of theratio of triton to proton axial matrix elements. As expected from the deep binding of thissystem, the volume eﬀects are small. The resulting inﬁnite-volume GT matrix element is (cid:104) GT (cid:105) L = ∞ = g H A /g pA = 0 . V. DISCUSSION

In order to connect to the physical quark masses, the result for the ratio of triton andproton axial charges described in the last section can be combined with the previous de-termination of g H A /g pA = 0 . SU (3)-symmetric point (where m π = m K = 806 MeV) using the same action and scale setting procedure. Tritium β -decayhas been investigated in χ PT in Ref. [50] (see Refs. [51, 52] for related work in pionless ef-fective ﬁeld theory) and so the mass dependence of the ratio is in principle known. However,the quark masses in the calculation of Ref. [18], and potentially in the current work, arebeyond the regime of applicability of χ PT. Additionally, for the three-body system of thetriton, the eﬀective ﬁeld theory results are determined numerically rather than analytically.At present, the above discussion, and the paucity of LQCD data, motivates extrapolationof the axial current matrix element ratio with the phenomenological forms of linear andquadratic dependence on the pion mass. The calculated GT matrix elements and bothof these extrapolations are shown in Fig. 8; the two ﬁts result in values of 0.90(8) and0.92(6), respectively. Given the model dependence of the forms used in this extrapolation,the envelope of the extrapolated uncertainties is taken as the extrapolated result, leading to g H A /g pA = 0 . +0 . − . .While extrapolated to inﬁnite volume and the physical quark masses, these results aredetermined at a single lattice spacing and QED and isospin breaking eﬀects are absent andthe uncertainties from these systematic eﬀects can as yet only be estimated. Lattice spacingeﬀects are expected to contribute to the matrix elements at O ( a Λ QCD ), however it is likelythat there will be partial cancellations in these eﬀects between the proton and triton axialcharges and a full evaluation of this uncertainty will require further calculations. The leadingQED eﬀects cancel in the ratio of triton to proton axial charges and isospin breaking eﬀectsin g pA have been estimated as (cid:46) e − m π L ∼ − so these eﬀects are expected to be negligible.The axial charge of the triton is thus extracted from lattice QCD for the ﬁrst time inthis work. The extrapolated coupling ratio is in agreement with the phenomenological value g H A /g pA = 0 . β -decay matrix element. In the future, the two-body pionless EFTcurrents that are determined using these methods will also allow for calculations of the GTmatrix elements in larger nuclei and a more comprehensive investigation of the QCD originof the phenomenological quenching of the axial charge. ACKNOWLEDGMENTS

We are grateful to S. R. Beane, Z. Davoudi, K. Orginos, M. J. Savage, and B. Tiburzifor extensive discussions and collaboration in the early stages of this work. This researchused resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge NationalLaboratory, which is supported by the Oﬃce of Science of the U.S. Department of Energyunder Contract number DE-AC05-00OR22725, as well as facilities of the USQCD Collabo-ration, which are funded by the Oﬃce of Science of the U.S. Department of Energy, and theresources of the National Energy Research Scientiﬁc Computing Center (NERSC), a U.S. De-partment of Energy Oﬃce of Science User Facility operated under Contract No. DE-AC02-6 □ □ ** FIG. 8. Ratio of the axial charge of tritium to that of the single nucleon as a function of thepion mass. The result from this work and that of Ref. [18] are shown as the blue points while thephysical value [6] is shown in red at the physical pion mass (indicated by the vertical line).

Appendix A: FVEFT

Following Ref. [49], the stochastic variational method is used to connect the ﬁnite vol-ume axial matrix element to its inﬁnite volume limit in pionless EFT. The pionless EFTLagrangian relevant to the interaction of few-nucleon systems is L = L + L + L , (A1)where the strong interactions between nucleons arise from L = N † (cid:18) iD + D M N (cid:19) N + . . . , (A2) L = − C S (cid:0) N T P i N (cid:1) † (cid:0) N T P i N (cid:1) − C T (cid:0) N T P a N (cid:1) † (cid:0) N T P a N (cid:1) + . . . , (A3) L = − D N † N ) + . . . , (A4)where P i and P a are the projectors deﬁned in Eq. (18), M N is the nucleon mass, and C S , C T , and D are the relevant two and three-nucleon low-energy constants (LECs). Thetwo-nucleon couplings can also be expressed in terms of alternative LECs C , through therelations C T = C − C and C S = C + C . (A5)The n -particle non-relativistic Hamiltonian corresponding to Eq. (A1) is H = − M N (cid:88) i ∇ i + (cid:88) i

50 0 - - - - - - - - -

50 0 - - - - - FIG. 9. Determination of the EFT two-body couplings, C S,T from the deuteron (left) and pp (right) energies. The horizontal bands correspond to the energy shifts from the LQCD calculationson the E24 (blue), E32 (orange), and E48 (green) ensembles. The curves in each ﬁgure showthe FVEFT two-body energy shifts for L ∈ { , , } (using the same color scheme) for theEFT cutoﬀ parameter r = 0 . wavefunctions optimized in the same volumes as the LQCD calculations, the LECs C S , C T ,and D of the pionless EFT Lagrangian can be constrained by matching the ﬁnite volumeenergies to the LQCD results, with the allowed range of LECs determined by a ﬁt to theconstraints from the three volumes. With the LECs ﬁxed, volume-extrapolated energies areobtained using wavefunctions optimized at inﬁnite volume.Figure 9 shows the determination of the two-body LECs from the LQCD results forthe deuteron and dineutron energy shifts in the three lattice volumes. The couplings areregulator-scale–dependent but the resulting energy shifts are not; calculations with cutoﬀparameter Λ = √ /r with r ∈ { . , . } fm result in indistinguishable results from thosewith r = 0 . D and leading to the inﬁnite-volume extrapolationpresented in Fig. 3.Using the optimized wavefunctions, matrix elements of the isovector axial-vector currentcan be computed for the same ﬁnite volume as the LQCD calculations to ﬁx the corre-sponding LECs, and the extrapolation to inﬁnite-volume can then be undertaken in thesame manner as that used for the binding energies themselves. The isovector axial-vectorcurrent is expressed in pionless EFT as given in Eq. (17) in the main text. In position space,the two-nucleon contribution proportional to L ,A is implemented using the same gaussianregulator as for the potential, speciﬁcally L ,A (cid:104)(cid:0) N T ( r i ) P i N ( r i ) (cid:1) † (cid:0) N ( r j ) T ¯ P a N ( r j ) (cid:1) + h.c. (cid:105) g Λ ( r ij ) . (A9)The optimised triton wavefunction corresponding to volume of the E32 ensemble is usedto compute the ﬁnite-volume axial matrix element in Eq. (19). Approximating the tritonwavefunction as a tensor product | Ψ H , s (cid:105) = | χ s (cid:105) ⊗ | ψ ( r , r , r ) (cid:105) , (A10)9 □ □ □□ □ □ □□□□◦ ◦ ◦ ◦ ◦ ◦ ◦◦ ◦ ◦ ◦ ◦ ◦ ◦ - - - - - FIG. 10. Extrapolation of the deuteron (green) and pp (purple) energies to the inﬁnite volumelimit. Data points denote lattice results (from Ref. [19]), displaced slightly horizontally for clarity,while the shaded bands and inﬁnite-volume points show the result of the FVEFT extrapolation. - - - FIG. 11. Determination of the EFT three-body couplings D . As in Fig. 9, the horizontal bandscorrespond to the energy shifts from the LQCD calculations on the E24 (blue), E32 (orange), andE48 (green) ensembles. The curves in each ﬁgure show the corresponding FVEFT energy shifts forthe EFT cutoﬀ parameter r = 0 . the computation of the matrix elements separates into the spin-isospin and spatial parts.The simplest spin-isospin wavefunction for the spin-up component is given by | χ s =+ (cid:105) = 1 √ (cid:2)(cid:12)(cid:12) n ↑ p ↑ n ↓ (cid:11) − (cid:12)(cid:12) n ↓ p ↑ n ↑ (cid:11) − (cid:12)(cid:12) p ↑ n ↑ n ↓ (cid:11) + (cid:12)(cid:12) p ↑ n ↓ n ↑ (cid:11) − (cid:12)(cid:12) n ↑ n ↓ p ↑ (cid:11) + (cid:12)(cid:12) n ↓ n ↑ p ↑ (cid:11)(cid:3) , (A11)with an analogous expression for the spin-down wavefunction. The spatial part of the matrix0element is determined from the variationally-optimized triton wavefunction and is given by h H (Λ , L ) = (cid:82) (cid:81) k d r k (cid:80) i

450 MeV from lattice QCD,” Phys. Rev. D , 114512 (2015), [Erratum:Phys. Rev. D , 039903 (2020)], arXiv:1508.07583 [hep-lat].[20] M. L¨uscher and P. Weisz, “On-Shell Improved Lattice Gauge Theories,” Commun. Math.Phys. , 59 (1985), [Erratum: Commun. Math. Phys.98,433(1985)].[21] B. Sheikholeslami and R. Wohlert, “Improved Continuum Limit Lattice Action for QCD withWilson Fermions,” Nucl. Phys. B259 , 572 (1985).[22] Marc Illa et al. , “Low-energy Scattering and Eﬀective Interactions of Two Baryons at m π ∼

450 MeV from Lattice Quantum Chromodynamics,” (2020), arXiv:2009.12357 [hep-lat].[23] William Detmold and Kostas Orginos, “Nuclear correlation functions in lattice QCD,” Phys.Rev.

D87 , 114512 (2013), arXiv:1207.1452 [hep-lat].[24] S. R. Beane, E. Chang, S. D. Cohen, William Detmold, H. W. Lin, T. C. Luu, K. Orginos,A. Parre˜no, M. J. Savage, and A. Walker-Loud (NPLQCD), “Light Nuclei and Hypernucleifrom Quantum Chromodynamics in the Limit of SU(3) Flavor Symmetry,” Phys. Rev.

D87 ,034506 (2013), arXiv:1206.5219 [hep-lat].[25] M. Falcioni, M.L. Paciello, G. Parisi, and B. Taglienti, “Again on SU(3) glueball mass,” Nucl.Phys. B , 624–632 (1985).[26] A. Francis, J.R. Green, P.M. Junnarkar, Ch. Miao, T.D. Rae, and H. Wittig, “Lattice QCDstudy of the H dibaryon using hexaquark and two-baryon interpolators,” Phys. Rev. D ,074505 (2019), arXiv:1805.03966 [hep-lat].[27] Ben H¨orz et al. , “Two-nucleon S-wave interactions at the SU (3) ﬂavor-symmetric point with m ud (cid:39) m phys s : a ﬁrst lattice QCD calculation with the stochastic Laplacian Heaviside method,”(2020), arXiv:2009.11825 [hep-lat].[28] Michael L. Wagman, Frank Winter, Emmanuel Chang, Zohreh Davoudi, William Detmold,Kostas Orginos, Martin J. Savage, and Phiala E. Shanahan, “Baryon-Baryon Interactions andSpin-Flavor Symmetry from Lattice Quantum Chromodynamics,” Phys. Rev. D , 114510(2017), arXiv:1706.06550 [hep-lat].[29] S.R. Beane et al. , “Charged multi-hadron systems in lattice QCD+QED,” (2020),arXiv:2003.12130 [hep-lat].[30] H. Akaike, “A new look at the statistical model identiﬁcation,” IEEE Transactions on Auto-matic Control , 716–723 (1974).[31] Charles Stein, “Inadmissibility of the usual estimator for the mean of a multivariate normaldistribution,” in Proceedings of the Third Berkeley Symposium on Mathematical Statisticsand Probability, Volume 1: Contributions to the Theory of Statistics (University of California Press, Berkeley, Calif., 1956) pp. 197–206.[32] Olivier Ledoit and Michael Wolf, “A well-conditioned estimator for large-dimensional covari-ance matrices,” Journal of Multivariate Analysis , 365 – 411 (2004).[33] Enrico Rinaldi, Sergey Syritsyn, Michael L. Wagman, Michael I. Buchoﬀ, Chris Schroeder,and Joseph Wasem, “Lattice QCD determination of neutron-antineutron matrix elements withphysical quark masses,” Phys. Rev. D99 , 074510 (2019), arXiv:1901.07519 [hep-lat].[34] Gene Golub and Victor Pereyra, “Separable nonlinear least squares: the variable projectionmethod and its applications,” Inverse Problems , R1 (2003).[35] Dianne P. O’Leary and Bert W. Rust, “Variable projection for nonlinear least squares prob-lems,” Computational Optimization and Applications , 579–593 (2013).[36] A. Ali Khan et al. (QCDSF-UKQCD), “The Nucleon mass in N(f) = 2 lattice QCD: Finitesize eﬀects from chiral perturbation theory,” Nucl. Phys. B , 175–194 (2004), arXiv:hep-lat/0312030.[37] Silas R. Beane, “Nucleon masses and magnetic moments in a ﬁnite volume,” Phys. Rev. D , 034507 (2004), arXiv:hep-lat/0403015.[38] S.R. Beane, E. Chang, W. Detmold, H.W. Lin, T.C. Luu, K. Orginos, A. Parre˜no, M.J.Savage, A. Torok, and A. Walker-Loud, “High Statistics Analysis using Anisotropic CloverLattices: (IV) Volume Dependence of Light Hadron Masses,” Phys. Rev. D , 014507 (2011),arXiv:1104.4101 [hep-lat].[39] Moti Eliyahu, Betzalel Bazak, and Nir Barnea, “Extrapolating Lattice QCD Results usingEﬀective Field Theory,” Phys. Rev. C , 044003 (2020), arXiv:1912.07017 [nucl-th].[40] Raul A. Brice˜no, Jozef J. Dudek, and Ross D. Young, “Scattering processes and resonancesfrom lattice QCD,” Rev. Mod. Phys. , 025001 (2018), arXiv:1706.06223 [hep-lat].[41] Takeshi Yamazaki, Ken-ichi Ishikawa, Yoshinobu Kuramashi, and Akira Ukawa, “Heliumnuclei, deuteron and dineutron in 2+1 ﬂavor lattice QCD,” Phys. Rev. D , 074514 (2012),arXiv:1207.4277 [hep-lat].[42] Takeshi Yamazaki, Ken-ichi Ishikawa, Yoshinobu Kuramashi, and Akira Ukawa, “Study ofquark mass dependence of binding energy for light nuclei in 2+1 ﬂavor lattice QCD,” Phys.Rev. D , 014501 (2015), arXiv:1502.04182 [hep-lat].[43] Phiala E. Shanahan, Brian C. Tiburzi, Michael L. Wagman, Frank Winter, Emmanuel Chang,Zohreh Davoudi, William Detmold, Kostas Orginos, and Martin J. Savage, “Isotensor AxialPolarizability and Lattice QCD Input for Nuclear Double- β Decay Phenomenology,” Phys.Rev. Lett. , 062003 (2017), arXiv:1701.03456 [hep-lat].[44] Brian C. Tiburzi, Michael L. Wagman, Frank Winter, Emmanuel Chang, Zohreh Davoudi,William Detmold, Kostas Orginos, Martin J. Savage, and Phiala E. Shanahan, “Double- β Decay Matrix Elements from Lattice Quantum Chromodynamics,” Phys. Rev.

D96 , 054505(2017), arXiv:1702.02929 [hep-lat].[45] Chris Bouchard, Chia Cheng Chang, Thorsten Kurth, Kostas Orginos, and Andre Walker-Loud, “On the Feynman-Hellmann Theorem in Quantum Field Theory and the Calculationof Matrix Elements,” Phys. Rev.

D96 , 014504 (2017), arXiv:1612.06963 [hep-lat].[46] Zohreh Davoudi, William Detmold, Kostas Orginos, Assumpta Parre˜no, Martin J. Savage,Phiala Shanahan, and Michael L. Wagman, “Nuclear matrix elements from lattice QCD forelectroweak and beyond-Standard-Model processes,” (2020), arXiv:2008.11160 [hep-lat].[47] William Detmold, Marc Illa, David J. Murphy, Patrick Oare, Kostas Orginos, Phiala E.Shanahan, Michael L. Wagman, and Frank Winter, “Lattice QCD constraints on the partondistribution functions of He,” (2020), arXiv:2009.05522 [hep-lat]. [48] Boram Yoon et al. , “Isovector charges of the nucleon from 2+1-ﬂavor QCD with cloverfermions,” Phys. Rev. D95 , 074508 (2017), arXiv:1611.07452 [hep-lat].[49] W. Detmold and P. E. Shanahan, “Few-nucleon matrix elements in pionless eﬀective ﬁeldtheory in a ﬁnite volume,” (2021).[50] A. Baroni et al. , “Local chiral interactions, the tritium Gamow-Teller matrix element, and thethree-nucleon contact term,” Phys. Rev.

C98 , 044003 (2018), arXiv:1806.10245 [nucl-th].[51] Hilla De-Leon, Lucas Platter, and Doron Gazit, “Tritium β -decay in pionless eﬀective ﬁeldtheory,” (2016), arXiv:1611.10004 [nucl-th].[52] Hilla De-Leon, Lucas Platter, and Doron Gazit, “Calculation of an A = 3 bound-state matrixelement in pionless eﬀective ﬁeld theory,” (2019), arXiv:1902.07677 [nucl-th].[53] Xue-min Jin, “Isospin breaking in the nucleon isovector axial charge from QCD sum rules,”(1996), arXiv:hep-ph/9602298.[54] Norbert Kaiser, “Isospin breaking in neutron beta decay and SU(3) violation in semileptonichyperon decays,” Phys. Rev. C , 028201 (2001), arXiv:nucl-th/0105043.[55] Robert G. Edwards and Balint Jo´o (SciDAC, LHPC, UKQCD), “The Chroma software systemfor lattice QCD,” Lattice ﬁeld theory. Proceedings, 22nd International Symposium, Lattice2004, Batavia, USA, June 21-26, 2004 , Nucl. Phys. Proc. Suppl. , 832 (2005), [,832(2004)],arXiv:hep-lat/0409003 [hep-lat].[56] A. Pochinsky, “Qlua, https://usqcd.lns.mit.edu/qlua.”.[57] M.A. Clark, R. Babich, K. Barros, R.C. Brower, and C. Rebbi, “Solving Lattice QCD systemsof equations using mixed precision solvers on GPUs,” Comput. Phys. Commun. , 1517–1528 (2010), arXiv:0911.3191 [hep-lat].[58] R. Babich, M.A. Clark, B. Joo, G. Shi, R.C. Brower, and S. Gottlieb, “Scaling Lattice QCDbeyond 100 GPUs,” in