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 Artificial Intelligence and Fundamental Interactions Fermi National Accelerator Laboratory, Batavia, IL 60510, USA Jefferson Laboratory, 12000 Jefferson 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 field configurations 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 infinite volume using finite-volume pionless effective fieldtheory (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 field 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 unmodified from that of the proton),thereby demonstrating that QCD can explain the modification 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 finds 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 effects 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 fields, γ µ are Dirac matrices, and τ a are Paulimatrices in flavor space. In particular, the Gamow-Teller matrix element is defined 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 field configurations 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 configurations on which a total of N meas measurements were made (slightlydifferent numbers of sources were used on each configuration). 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 first studies in Ref. [16],see Ref. [17] for a recent summary of results. A first 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 infinite 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 modification 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-figurations generated with a L¨uscher-Weisz [20] gauge action and N f = 2 + 1 flavors 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 efficiently 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 significantlymore 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 gaugefield configurations, as discussed in Refs. [19, 22].The ground-state triton energy and its difference 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 fitting 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 effects For the L = { , , } ensembles, smearing parameters of N smear = { , , } steps of smearing withGaussian smearing parameters ρ = { . , . , . } were used. arising from the finite 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 fit parameters; while correlation functions for different propagator smearingshave different overlap factors, the energies are common and are thus fit simultaneously. Toquantify the systematic uncertainties arising from the choice of source-sink separations, t ,included in the fit, and of the truncation of the sum in Eq. (5), 200 fit windows are sampledat random from the set of all choices of contiguous time-ranges longer than t plat = 6 and withfinal times less than t max (defined by the point at which the noise-to-signal ratio exceeds tol noise = 0 .
25 for the given correlation function). For each fit range, the Akaike InformationCriterion [30] (AIC) is used to perform model selection ( i.e. , fix 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 fit. Ineach case, combined correlated fits 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 fits with a χ /N dof < tol χ = 2 are included in a set of‘accepted’ fits (accepted fits must also pass tests that the fit 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 fits, and c) agree within a prescribed toleranceof tol median = 2 σ to the median bootstrap result, as in Refs. [22, 29]). The final value anduncertainties for the energy are then computed from the results of all N success accepted fitsusing 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 fit result from a given fit labelled by f , and w f is the associatedweight. For each fit range, the statistical uncertainties δE f are computed using 67% confi-dence intervals determined from the results of the N boot resampled correlation function fitsdescribed above; the total statistical uncertainty is defined as the statistical uncertainty ofthe fit 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 fit provides a statistically unbiased estimate of the ground-state energy, the rela-tive weights w f of each fit 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 fit f assuming χ -distributedgoodness-of-fit parameters.The resulting fitted masses are summarised for each volume in Table II and the fits areshown graphically for the proton in Fig. 1 and for the triton in Fig. 2. In each figure, the □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ FIG. 1. Left column: effective mass functions (Eq. (8)) for the proton correlation function on theE24 (top), E32 (center), and E48 (bottom) ensembles. The orange circles and blue squares showeffective masses constructed from SP and SS correlation functions respectively, and the correspond-ing orange and blue curves display the highest-weight fit contributing to the weighted average ofEq. (6); the light (dark) grey bands show the mass extracted from the weighted average (singlehighest-weight fit). Right column: masses extracted from each successful fit to the correlators.The opacity of each point is set by the contribution of the fit to the weighted average. The hor-izontal line and grey band in each figure correspond to the final result for the mass and its totaluncertainty. effective 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 fits to the correlationfunctions and the final result for E ( h )0 . The difference between the triton mass and three □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □ □ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ □ □ □ □ □ □ □ □ □ □ □ □○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ ○ FIG. 2. Effective mass functions and fit summaries for the fits 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 differencesof the fitted energies computed under bootstrap. The results for this quantity are alsopresented in Table II.In the limit in which m π L is sufficiently large and the pion mass is sufficiently 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 modified 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 effectively 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 fit to the infinite volumeproton mass from the masses determined on the three volumes. This fit, displayed in Figure3, results in an infinite volume mass of M ∞ p = 1 . Figure 3 also shows the difference 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 finite-volume effective field 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-finite 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 finite-volume three-particle energies to the infinite-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 fit 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 finite-volume dependence inpionless EFT fitted to the LQCD data. The blue point at the right indicates the infinite 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 effort is made to takeinto account the differences 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 figure, 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 fits 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 flavor-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 infinite-volume matrix element.To compute the finite-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 field that couples toup and down quarks separately. Extended quark propagators are defined 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 flavor f ∈ { u, d } and λ f is the strength of theapplied field for the given flavor. These quantities are calculated for five values of the exter-nal field 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-field 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 flavor f inthe correlation function. With computations at n f different field 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 effective isovector axial charge function, which asymptotesto the desired bare axial charge in the corresponding hadronic or nuclear state at largetemporal separations, can be defined 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 first 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 differ at the source and sink, and the bare(transition) charge g h ( m,n )(0) A is defined 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 fits of two-point functions to Eq. (5). Notably, the O ( a ) term is absent for a single-state correlation The correlation function is defined 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 effective matrix element functions, g pA ( t )and g H A ( t ), respectively, calculated on the E32 ensemble. In each figure, the orange (blue) datacorrespond to the SP (SS) calculations and the corresponding shaded bands illustrate the highestweight fit. The light (dark) shaded gray bands denote the extracted values of the matrix elementsarising from the combined analysis (single highest-weight fit). The right panels show the valuesobtained for all the successful fits to different time ranges, with the opacity determined by thecontribution of the fit in the final weighted average. The horizontal line and gray band in the rightpanels show the final central value and uncertainty. function model. Multi-state fits 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 thefit. For both the triton and proton, this AIC test prefers the fit without O ( a ) terms for all fitrange choices. Combined fits 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 fitting 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 fit are obtained using bootstrapconfidence intervals, and a weighted average performed analogously to the two-point functioncase described above is used to determine the final ground-state matrix element values andstatistical plus fitting systematic uncertainties. Results for the ground-state matrix elements g h (0) A obtained using this fitting 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 effective 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 fit 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 fitting systematic uncertainties are combined in quadrature. The fitting systematicuncertainties account for variation in the choice of fit ranges and multi-state fit 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 fits 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 fits 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 effective field theory in the finite 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 defined 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 infinite 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 different volumes, includingthe infinite volume wave function that allows the infinite-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 effects are small. The resulting infinite-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 field 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 effective field 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 fits 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 infinite volume and the physical quark masses, these results aredetermined at a single lattice spacing and QED and isospin breaking effects are absent andthe uncertainties from these systematic effects can as yet only be estimated. Lattice spacingeffects are expected to contribute to the matrix elements at O ( a Λ QCD ), however it is likelythat there will be partial cancellations in these effects between the proton and triton axialcharges and a full evaluation of this uncertainty will require further calculations. The leadingQED effects cancel in the ratio of triton to proton axial charges and isospin breaking effectsin g pA have been estimated as (cid:46) e − m π L ∼ − so these effects are expected to be negligible.The axial charge of the triton is thus extracted from lattice QCD for the first 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 Office 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 Office of Science of the U.S. Department of Energy, and theresources of the National Energy Research Scientific Computing Center (NERSC), a U.S. De-partment of Energy Office 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 finite vol-ume axial matrix element to its infinite 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 defined 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 figure showthe FVEFT two-body energy shifts for L ∈ { , , } (using the same color scheme) for theEFT cutoff 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 finite volumeenergies to the LQCD results, with the allowed range of LECs determined by a fit to theconstraints from the three volumes. With the LECs fixed, volume-extrapolated energies areobtained using wavefunctions optimized at infinite 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 cutoffparameter Λ = √ /r with r ∈ { . , . } fm result in indistinguishable results from thosewith r = 0 . D and leading to the infinite-volume extrapolationpresented in Fig. 3.Using the optimized wavefunctions, matrix elements of the isovector axial-vector currentcan be computed for the same finite volume as the LQCD calculations to fix the corre-sponding LECs, and the extrapolation to infinite-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, specifically 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 finite-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 infinite volumelimit. Data points denote lattice results (from Ref. [19]), displaced slightly horizontally for clarity,while the shaded bands and infinite-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 figure show the corresponding FVEFT energy shifts forthe EFT cutoff 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 Effective 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) flavor-symmetric point with m ud (cid:39) m phys s : a first 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 identification,” 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. Buchoff, 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 effects 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 finite 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 usingEffective 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 flavor 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 flavor 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-flavor 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 effective fieldtheory in a finite 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 effective fieldtheory,” (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 effective field 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 field 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