First-principles calculation of electroweak box diagrams from lattice QCD
Xu Feng, Mikhail Gorchtein, Lu-Chang Jin, Peng-Xiang Ma, Chien-Yeah Seng
MMITP/20-018
First-principles calculation of electroweak box diagrams from lattice QCD
Xu Feng,
1, 2, 3, 4
Mikhail Gorchtein,
5, 6, 7
Lu-Chang Jin,
8, 9
Peng-Xiang Ma, and Chien-Yeah Seng School of Physics, Peking University, Beijing 100871, China Collaborative Innovation Center of Quantum Matter, Beijing 100871, China Center for High Energy Physics, Peking University, Beijing 100871, China State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, China Helmholtz Institute Mainz, Mainz, Germany GSI Helmholtzzentrum f¨ur Schwerionenforschung, Darmstadt, Germany Johannes Gutenberg University, Mainz, Germany Department of Physics, University of Connecticut, Storrs, CT 06269, USA RIKEN-BNL Research Center, Brookhaven National Laboratory, Building 510, Upton, NY 11973 Helmholtz-Institut f¨ur Strahlen- und Kernphysik and Bethe Center for Theoretical Physics,Universit¨at Bonn, 53115 Bonn, Germany (Dated: May 4, 2020)We present the first realistic lattice QCD calculation of the γW -box diagrams relevant for betadecays. The nonperturbative low-momentum integral of the γW loop is calculated using a latticeQCD simulation, complemented by the perturbative QCD result at high momenta. Using the pionsemileptonic decay as an example, we demonstrate the feasibility of the method. By using domainwall fermions at the physical pion mass with multiple lattice spacings and volumes, we obtain theaxial γW -box correction to the semileptonic pion decay, ◻ V AγW ∣ π = . ( ) stat ( ) sys × − , withthe total uncertainty controlled at the level of ∼ γW -box correction to the neutron decay, which plays a decisive role in thedetermination of ∣ V ud ∣ . Introduction – The precise determination of theCabibbo-Kobayashi-Maskawa (CKM) matrix elements,which are fundamental parameters of the StandardModel, is one of the central themes in modern par-ticle physics. In the CKM matrix, V ud is the mostaccurately-determined element from the study of superal-lowed 0 + → + nuclear beta decays [1] which are pure vec-tor transitions at tree level and are theoretically clean dueto the protection of the conserved vector current. Goingbeyond tree level, the electroweak radiative correctionsinvolving the axial-vector current become important andultimately dominate the theoretical uncertainties.Among various electroweak radiative corrections, theaxial γW -boson box contribution ◻ V AγW contains a signif-icant sensitivity to low-energy hadronic effects, and is adominant source of the total theoretical uncertainty [2].The recent dispersive analysis [3, 4] reduced this uncer-tainty by a factor of 2 comparing to the previous study byMarciano and Sirlin [5], and the updated result of ∣ V ud ∣ raised a 4 standard-deviation tension with the first-rowCKM unitarity (barring possibly underestimated nucleareffects: see Ref. [4, 6]). The main difference betweenthose works is the use of inclusive neutrino and antineu-trino scattering data that Refs. [3, 4] used to estimatethe contribution of the intermediate momenta inside the γW loop integral, 0 . ≲ Q ≲ , prone to non-perturbative hadronic effects. To further improve thedetermination of ∣ V ud ∣ , it requires either better-qualityexperimental input or the direct, precise lattice QCD cal-culations of the γW -box contribution.Lattice QCD has played an important role in the de-termination of the nonperturbative hadronic matrix el- ements needed to constrain the CKM unitarity. Recentlattice results are averaged and summarized by the FLAGreport 2019 [7]. With lattice QCD simulations havingreached an impressive level of precision for tree-levelparameters of the electroweak interaction, it becomestimely and important to study higher-order electroweakcorrections. The examples of such lattice applications in-clude the QED corrections to hadron masses [8–15] andleptonic decay rates [16–19] and a series of higher-orderelectroweak effects, such as K L - K S mass difference [20–22], (cid:15) K [23], rare kaon decays [24–29] and double betadecays [30–35]. As for the γW -box contribution, whichis a QED correction to semileptonic decays, it still re-mains a new horizon for lattice QCD.It has been proposed to use the Feynman-Hellmanntheorem to calculate the γW -box contribution [36, 37].In this work, we opt for a more straightforward way toperform the lattice calculation. To demonstrate the fea-sibility of the method, we carry out the exploratory studyfor the case of the pion semileptonic decays. The calcula-tion is performed at the physical pion mass with variouslattice spacings and volumes, which allows us to controlthe systematic effects in the lattice results. Combiningthe results from lattice QCD together with the pertur-bative QCD, we obtain the axial γW -box correction topion decay amplitude with a relative ∼
1% uncertainty.
The γW -box contribution – In the theoreti-cal analysis of the superallowed nuclear beta decayrates, the dominant uncertainty arises from the nucleus-independent electroweak radiative correction, ∆ VR , whichis universal for both nuclear and free neutron beta de-cay [1]. Among various contributions to ∆ VR , Sirlin es- a r X i v : . [ h e p - l a t ] M a y tablished [2] that only the axial γW -box contribution issensitive to hadronic scales; see Fig. 1 for the γW dia-grams. The relevant hadronic tensor T V Aµν is defined as T V Aµν = ∫ d x e iqx ⟨ H f ( p )∣ T [ J emµ ( x ) J W,Aν ( )] ∣ H i ( p )⟩ , (1)for a semileptonic decay process H i → H f e ¯ ν e . Above, H i / f are given by neutron and proton for the neutronbeta decay, and by π − and π for the pion semileptonicdecay, respectively. Furthermore, J emµ = ¯ uγ µ u − ¯ dγ µ d − ¯ sγ µ s is the electromagnetic quark current, and J W,Aν = ¯ uγ ν γ d is the axial part of the weak charged current. Figure 1. The γW -box diagrams for the semileptonic decayprocess H i → H f e ¯ ν e . The spin-independent part of T V Aµν has only one term, T V Aµν = i(cid:15) µναβ q α p β T + . . . , where T is a scalar function.For the neutron beta decay, the spin-dependent contri-butions, denoted by the ellipses here, are absorbed intothe definition of the nucleon axial charge g A , which canbe measured directly from experiments. According tocurrent algebra [2], it is this spin-independent term thatgives rise to the hadron structure-dependent contributionand dominates the uncertainty in the theoretical predic-tion. Using T as input, the axial γW -box correction tothe tree-level amplitude is given as [3] ◻ V AγW ∣ H = F H + α e π ∫ ∞ dQ m W m W + Q × ∫ √ Q −√ Q dQ π ( Q − Q ) ( Q ) T ( Q , Q ) . (2)Here Q = − q > F H + arises from the lo-cal matrix element ⟨ H f ( p ′ )∣ J W,Vµ ∣ H i ( p )⟩ = ( p + p ′ ) µ F H + ,with F H + = √ Methodology – In the framework of lattice QCD, thehadronic tensor T V Aµν in Euclidean spacetime is given by T V Aµν = ∫ dt e − iQ t ∫ d x e − i ⃗ Q ⋅⃗ x H V Aµν ( t, ⃗ x ) (3)with H V Aµν ( t, ⃗ x ) defined as H V Aµν ( t, ⃗ x ) ≡ ⟨ H f ( P )∣ T [ J emµ ( t, ⃗ x ) J W,Aν ( )] ∣ H i ( P )⟩ . (4)Here the Euclidean momenta P and Q are chosen as P = ( im H , ⃗ ) , Q = ( Q , ⃗ Q ) (5) with m H the hadron mass.By multiplying (cid:15) µναβ Q α P β to T V Aµν , we can extract thefunction T ( Q , Q ) through T ( Q , Q ) = − I m H ∣ ⃗ Q ∣ , I = (cid:15) µναβ Q α P β T V Aµν . (6)Here I can be written in terms of H V Aµν as I = i (cid:15) µνα Q α m H ∫ dt e − iQ t ∫ d ⃗ x e − i ⃗ Q ⋅⃗ x H V Aµν = m H ∫ dt e − iQ t ∫ d ⃗ x e − i ⃗ Q ⋅⃗ x (cid:15) µνα ∂ H V Aµν ∂x α . (7)We can average over the spatial directions for ⃗ Q and have I = m H ∫ dt e − iQ t ∫ d ⃗ x j (∣ ⃗ Q ∣∣⃗ x ∣) (cid:15) µνα ∂ H V Aµν ∂x α = m H ∫ dt e − iQ t ∫ d ⃗ x ∣ ⃗ Q ∣∣⃗ x ∣ j (∣ ⃗ Q ∣∣⃗ x ∣) (cid:15) µνα x α H V Aµν , (8)where j n ( x ) are the spherical Bessel functions. A keyingredient in this approach is that once the Lorentz scalarfunction (cid:15) µνα x α H V Aµν is prepared, e.g. from a latticeQCD calculation, one can determine T ( Q , Q ) directly.Putting Eqs. (8) and (6) into Eq. (2) and changingthe variables as ∣ ⃗ Q ∣ = √ Q cos θ and Q = √ Q sin θ , weobtain the master formula ◻ V AγW ∣ H = α e π ∫ dQ Q m W m W + Q M H ( Q ) (9)with M H ( Q ) = −
16 1 F H + √ Q m H ∫ d x ω ( t, ⃗ x ) (cid:15) µνα x α H V Aµν ( t, ⃗ x ) ,ω ( t, ⃗ x ) = ∫ π − π cos θ dθπ j (√ Q ∣⃗ x ∣ cos θ )∣⃗ x ∣ cos (√ Q t sin θ ) . (10)For small Q , lattice QCD can determine the function M H ( Q ) with lattice discretization errors under control.For large Q , we utilize the operator product expansion12 ∫ d xe − iQx T [ J emµ ( x ) J W,Aν ( )]= i Q { C a ( Q ) δ µν Q α − C b ( Q ) δ µα Q ν − C c ( Q ) δ να Q µ } J W,Aα ( )+ Q C d ( Q ) (cid:15) µναβ Q α J W,Vβ ( ) + ⋯ . (11)There are only four possible local operators at leadingtwist. (For the pion decay, the hadronic matrix ele-ments for the first three operators vanish.) Multiplying (cid:15) µναβ Q α P β to the relation (11) we obtain T ( Q , Q ) = C d ( Q ) Q F H + + ⋯ ,M H ( Q ) = C d ( Q ) + ⋯ , (12)where the ellipses remind us that the higher-twist con-tributions are not included yet. The Wilson coefficient C d ( Q ) is calculated to four-loop accuracy [38, 39] C d ( Q ) = ∑ n c n a ns , a s = α s ( Q ) π , (13)with coefficients c n given in Eq. (12) of Ref. [39]. Here α s is the strong coupling constant.We introduce a momentum-squared scale Q thatseparates the two regimes, and split the integral in Eq. (9)into two parts ◻ V AγW ∣ H = ◻ V A, ≤ γW ∣ H + ◻ V A, > γW ∣ H = (∫ Q dQ Q + ∫ ∞ Q dQ Q ) m W m W + Q M H ( Q ) . (14)With Eq. (10) we use the lattice data to determine theintegral for Q ≤ Q , while with Eq. (12) we use per-turbation theory to determine the integral for Q > Q . π − π J emµ J W,Aν (A) π − π J emµ J W,Aν (B) π − π J emµ J W,Aν (C) π − π J emµ J W,Aν (D)
Figure 2. Four types of quark contractions for pion γW -boxdiagrams. Lattice setup – We use five lattice QCD gauge en-sembles at the physical pion mass, generated by RBCand UKQCD Collaborations using 2 + ⟨ φ π ( t f ) J emµ ( x ) J W,Aν ( y ) φ † π − ( t i )⟩ with t f = max { t x , t y } + ∆ t and t i = min { t x , t y } − ∆ t . We use the wall-source pioninterpolating operators φ π and φ † π − , which have a good overlap with the π ground state, and find the ground-state saturation for ∆ t ≳ t are chosen conservatively as shown in Table I. Foreach ensemble we use the gauge configurations with suffi-ciently long separation, i.e. each separated by at least 10trajectories. The number of configurations used is listedin Table I. Ensemble m π [MeV] L T a − [GeV] N conf N r ∆ t / a
24D 141.2(4) 24 64 1 .
015 46 1024 832D 141.4(3) 32 64 1 .
015 32 2048 832D-fine 143.0(3) 32 64 1 .
378 71 1024 1048I 135.5(4) 48 96 1 .
730 28 1024 1264I 135.3(2) 64 128 2 .
359 62 1024 18Table I. Ensembles used in this work. For each ensemblewe list the pion mass m π , the spatial and temporal extents, L and T , the inverse of lattice spacing a − , the number of config-urations used, N conf , the number of point-source light-quarkpropagator generated for each configuration, N r , and the timeseparation, ∆ t , used for the π ground-state saturation. There are four types of contractions for γW -box dia-grams as shown in Fig. 2. We produce wall-source quarkpropagators on all time slices. Using the techniques de-scribed in Ref. [35] type (A) and (B) diagrams can be cal-culated with high precision by performing the spacetime-translation average over L × T measurements. Under the γ Hermitian conjugation of the Euclidean quark propa-gators, one can confirm that type (B) does not contributeto the axial γW -box diagrams. Type (C) diagram is cal-culated by treating one current as the source and theother as the sink. We calculate point-source propagatorsat N r random spacetime locations. The values of N r areshown in Table I. These point-source propagators can beplaced at either electromagnetic current or weak current.We thus average the type (C) correlation functions over2 N r measurements. This is similar with the treatmenttaken by Ref. [41]. We neglect the disconnected contri-bution (D), which vanishes in the flavor SU ( ) limit. Numerical results – In practice, the integral inEq. (10) can be performed within a range of √ t + ⃗ x ≤ R .Taking the ensemble 64I as an example, M π ( Q ) as afunction of the integral range R is shown in Fig. 3. Wefind that for all the momenta Q ∈ [ , ] GeV , the in-tegral is saturated at large R . We choose the truncationrange R ≃ √ t + ⃗ x > R is negligible, indicating thatthe finite-volume effects are well under control in ourcalculation. We can further verify this conclusion by adirect comparison using the 24D and 32D data.The lattice results of M π ( Q ) as a function of Q areshown in Fig. 4 together with the perturbative results.Ensemble 24D and 32D have the same pion mass and lat-tice spacing but different volumes. The good agreementbetween these two ensembles indicates that the finite- M π ( Q ) Q = 0.05 GeV Q = 0.10 GeV Q = 0.20 GeV Q = 0.40 GeV Q = 0.60 GeV Q = 1.00 GeV Q = 4.00 GeV Figure 3. For ensemble 64I, lattice results of M π ( Q ) as afunction of the integral range R . [GeV ]00.020.040.060.08 M π ( Q ) [GeV ] Cont. Limit, IwasakiCont. Limit, DSDRPT (n f =4 match with n f =3)PT (n f =4 down to 1 GeV)Q = 2 GeV Figure 4. M π ( Q ) as a function of Q . In the left panel,the lattice results for ensembles 64I, 48I, 32D-fine, 32D and24D are represented by turquoise, indigo, dark green, red andblack bands, respectively. Taking 64I and 24D as examples,the results for type (A) diagram are also plotted. In the rightpanel, Iwasaki and DSDR results at the continuum limit areshown by the gray and brown bands. The orange curve showsthe results from perturbation theory by decoupling the charmquark at 1.6 GeV while the magenta one is compiled using the4-flavor perturbation theory continuously down to 1 GeV. volume effects are smaller than the statistical errors. At Q ≳ , the lattice discretization effects dominatethe uncertainties. In the left panel of Fig. 4 an obviousdiscrepancy is observed at large Q for the lattice resultswith different lattice spacings.For the perturbation theory, the Wilson coefficient C d ( Q ) is determined using the RunDec package [42],where α s is calculated to four-loop accuracy. At low Q the results still contain large systematic uncertainties dueto the lack of higher-loop and higher-twist contributions.In Fig. 4 we show two curves from perturbation theory.One is compiled using 4-flavor theory down to 1 GeV, while the other decouples the charm quark at 1.6 GeVand uses 3-flavor theory for ( ) ≤ Q ≤ ( ) .The discrepancy between the two curves suggests anO(14%) systematic effect in the perturbative determina-tion of M H ( Q ) at Q ≈ . Estimate of systematic effects – For ◻ V A, ≤ γW thelargest uncertainties arise from the lattice discretizationeffects. Since Iwasaki and DSDR ensembles have differ-ent lattice discretizations, we treat them separately. Af-ter the linear extrapolation in a , the Iwasaki and DSDRresults at continuum limit are shown in the right panelof Fig. 4. Using Q = we obtain ◻ V A, ≤ γW ∣ π = ⎧⎪⎪⎨⎪⎪⎩ . ( ) × − for Iwasaki0 . ( ) × − for DSDR . (15)We take the Iwasaki result as the central value and es-timate the residual O ( a ) lattice artifacts using the dis-crepancy between Iwasaki and DSDR.For ◻ V A, > γW the largest uncertainties arise from thehigher-loop and higher-twist truncation effects. We es-timate the former by comparing the 4-loop and 3-loopresults from perturbation theory. For the latter, unfor-tunately the complete information is not available. Con-sidering the fact that type (A) diagram, which has twocurrents located at different quarks lines, only containsthe higher-twist contributions, we use it to estimate thesize of higher twist. At Q = we have ◻ V A, > γW ∣ π = . ( ) HL ( ) HT × − , (16)where the central value is compiled using the 4-flavor the-ory (see the magenta curve in the right panel of Fig. 4).The first error indicates the higher-loop effects. The sec-ond one stands for the higher-twist effects, which arecompiled from the integral of Q > Q using the type(A) data as input. Summary of results – After combining the resultsof ◻ V A, ≤ γW from lattice QCD and ◻ V A, > γW from perturbationtheory, we obtain the total contribution of ◻ V AγW ◻ V AγW ∣ π = . ( ) stat ( ) PT ( ) a ( ) FV × − = . ( ) stat ( ) sys × − , (17)where the first uncertainty is statistical, and the re-maining errors account for perturbative truncationand higher-twist effects, lattice discretization effects,and lattice finite-volume effects by comparing the24D and 32D results. We add these systematic er-rors in quadrature to obtain the final systematic er-ror. For comparison, we also calculate ◻ V AγW ∣ π = . ( ) stat ( ) PT ( ) a ( ) FV × − at Q = and ◻ V AγW ∣ π = . ( ) stat ( ) PT ( ) a ( ) FV × − at Q = . Both results are consistent with Eq. (17).For the pion semileptonic decay, the PIBETA experi-ment [43] has improved the measurement of the branch-ing ratio to 0.6%. The Standard Model prediction of thedecay rate is given by [2, 43]Γ π(cid:96) = G F ∣ V ud ∣ m π ∣ f π + ( )∣ π ( + δ ) I π (18)with G F = . ( ) × − GeV − the Fermi’s con-stant measured from the muon decay, m π the chargedpion mass, f π + ( ) = I π = . ( ) × − a known kinematic fac-tor. Numerically, Γ πl = . ( ) s − after taking intoaccount the updated value of π + → e + ν ( γ ) branchingratio as an overall normalization [44]. The effects ofradiative corrections are contained in δ . The existinganalysis from chiral perturbation theory (ChPT) yields δ = . ( ) LEC ( ) HO [44–47] with an overall theoret-ical uncertainty of Γ π(cid:96) at a level of 0 . ∣ V ud ∣ = . ( ) exp ( ) th with a 0.3% uncertainty.We now show how our calculation reduces the uncer-tainty in δ . We adopt Sirlin’s parameterization [2] withslight modifications: δ = α e π [ ¯ g + m Z m p + ln m Z m W + ˜ a g ] + δ QEDHO + ◻ V AγW . (19)By separating the axial γW -box part into ◻ V AγW , the re-maining contributions are model independent and aregiven as follows. • Sirlin’s function ¯ g arises from a structure-indepenent, UV-finite one-loop integral. It ac-counts for the infrared contributions involving thevector γW -box and the bremsstrahlung correc-tions. It contains a 3 ln m p term that cancels the m p -dependence in 3 ln m Z m p . Here m p is the protonmass that appears just as a matter of convention.Numerically, one has α e π ¯ g = . × − [2, 49]. • ˜ a g represents the O ( α s ) QCD correction to all one-loop diagrams except for the axial γW box. Theintegral in ˜ a g is dominated by the high-energy scale Q ≃ m W , where α s is small. As a consequence α e π ˜ a g ≈ − . × − is a small contribution [2, 50]. • δ QEDHO = . ( ) summarizes the leading-loghigher-order QED effects which can be accountedfor through the running of α e . The uncertainty as-signment follows Ref. [48].Although the detailed uncertainties for ¯ g and ˜ a g are notgiven, by power counting the intrinsic precision for theterms in the square brackets (multiplied by α e π ) is of theorder G F m p ∼ − .Combining the ◻ V AγW in Eq. (17), we now obtain δ = . ( ) γW ( ) HO , (20) which corresponds to an almost complete removal ofthe dominant LEC uncertainties in the ChPT expres-sion, and a reduction of the total uncertainty by a fac-tor of 3. Therefore, any theoretical improvement in thefuture will unavoidably require a complete electroweaktwo-loop analysis. Consequently, the ∣ V ud ∣ determinedfrom the pion semileptonic decay now reads: ∣ V ud ∣ = . ( ) exp ( ) th . Conclusion – In this work we perform the first real-istic lattice QCD calculation of the γW -box correctionto the pion semileptonic decay, ◻ V AγW ∣ π . The final resultcombines the lattice data at low momentum and pertur-bative calculation at high momentum. We use multiplelattice spacings and volumes at the physical pion massto control the continuum and infinite-volume limits andobtain ◻ V AγW ∣ π with a total error of ∼ σ tension persists.The combined experimental measurement of 14 nuclearsuperallowed beta decays [1] is 10 times more accuratethan the current pion semileptonic decay experiment. Onthe other hand, the free neutron decay [52, 53] leads toa 4.5 times better precision. In these two cases, the non-perturbative, structure-dependent γW -box contributionplays a decisive role. The technique presented in thiswork can be straightforwardly generalized to a lattice cal-culation of the nucleon γW -box corrections, which areuniversal for both free and bound neutron decay. Thelatter is the key to a precise determination of ∣ V ud ∣ and astringent test of CKM unitarity. Acknowledgments – X.F. and L.C.J. gratefully ac-knowledge many helpful discussions with our colleaguesfrom the RBC-UKQCD Collaborations. We thank Yan-Qing Ma and Guido Martinelli for inspiring discussions.X.F. and P.X.M. were supported in part by NSFCof China under Grant No. 11775002. M.G. is sup-ported by EU Horizon 2020 research and innovation pro-gramme, STRONG-2020 project, under grant agreementNo 824093 and by the German-Mexican research collab-oration Grant No. 278017 (CONACyT) and No. SP778/4-1 (DFG). L.C.J. acknowledges support by DOEgrant DE-SC0010339. The work of C.Y.S. is supportedin part by the DFG (Grant No. TRR110) and the NSFC(Grant No. 11621131001) through the funds provided tothe Sino-German CRC 110 Symmetries and the Emer-gence of Structure in QCD, and also by the Alexandervon Humboldt Foundation through the Humboldt Re-search Fellowship. The computation is performed underthe ALCC Program of the US DOE on the Blue Gene/Q(BG/Q) Mira computer at the Argonne Leadership ClassFacility, a DOE Office of Science Facility supported underContract DE-AC02-06CH11357. Computations for thiswork were carried out in part on facilities of the USQCDCollaboration, which are funded by the Office of Scienceof the U.S. Department of Energy. The calculation is alsocarried out on Tianhe 3 prototype at Chinese NationalSupercomputer Center in Tianjin. [1] J. C. Hardy and I. S. Towner, Phys. Rev.
C91 , 025501(2015), arXiv:1411.5987 [nucl-ex].[2] A. Sirlin, Rev. Mod. Phys. , 573 (1978), [Erratum:Rev. Mod. Phys.50,905(1978)].[3] C.-Y. Seng, M. Gorchtein, H. H. Patel, and M. J.Ramsey-Musolf, Phys. Rev. Lett. , 241804 (2018),arXiv:1807.10197 [hep-ph].[4] C. Y. Seng, M. Gorchtein, and M. J. Ramsey-Musolf,Phys. Rev. D100 , 013001 (2019), arXiv:1812.03352 [nucl-th].[5] W. J. Marciano and A. Sirlin, Phys. Rev. Lett. , 032002(2006), arXiv:hep-ph/0510099 [hep-ph].[6] M. Gorchtein, Phys. Rev. Lett. , 042503 (2019),arXiv:1812.04229 [nucl-th].[7] S. Aoki et al. (Flavour Lattice Averaging Group), Eur.Phys. J. C80 , 113 (2020), arXiv:1902.08191 [hep-lat].[8] A. Duncan, E. Eichten, and H. Thacker, Phys. Rev. Lett. , 3894 (1996), arXiv:hep-lat/9602005 [hep-lat].[9] A. Duncan, E. Eichten, and H. Thacker, Phys. Lett. B409 , 387 (1997), arXiv:hep-lat/9607032 [hep-lat].[10] T. Blum, T. Doi, M. Hayakawa, T. Izubuchi, and N. Ya-mada, Phys. Rev.
D76 , 114508 (2007), arXiv:0708.0484[hep-lat].[11] T. Blum, R. Zhou, T. Doi, M. Hayakawa, T. Izubuchi,S. Uno, and N. Yamada, Phys. Rev.
D82 , 094508 (2010),arXiv:1006.1311 [hep-lat].[12] T. Ishikawa, T. Blum, M. Hayakawa, T. Izubuchi,C. Jung, and R. Zhou, Phys. Rev. Lett. , 072002(2012), arXiv:1202.6018 [hep-lat].[13] S. Borsanyi et al. , Science , 1452 (2015),arXiv:1406.4088 [hep-lat].[14] P. Boyle, V. G¨ulpers, J. Harrison, A. J¨uttner, C. Lehner,A. Portelli, and C. T. Sachrajda, JHEP , 153 (2017),arXiv:1706.05293 [hep-lat].[15] X. Feng and L. Jin, Phys. Rev. D100 , 094509 (2019),arXiv:1812.09817 [hep-lat].[16] N. Carrasco, V. Lubicz, G. Martinelli, C. T. Sachrajda,N. Tantalo, C. Tarantino, and M. Testa, Phys. Rev.
D91 , 074506 (2015), arXiv:1502.00257 [hep-lat].[17] V. Lubicz, G. Martinelli, C. T. Sachrajda, F. Sanfilippo,S. Simula, and N. Tantalo, Phys. Rev.
D95 , 034504(2017), arXiv:1611.08497 [hep-lat].[18] D. Giusti, V. Lubicz, G. Martinelli, C. T. Sachrajda,F. Sanfilippo, S. Simula, N. Tantalo, and C. Tarantino,Phys. Rev. Lett. , 072001 (2018), arXiv:1711.06537[hep-lat].[19] M. Di Carlo, D. Giusti, V. Lubicz, G. Martinelli, C. T.Sachrajda, F. Sanfilippo, S. Simula, and N. Tantalo,Phys. Rev.
D100 , 034514 (2019), arXiv:1904.08731 [hep-lat].[20] N. H. Christ, T. Izubuchi, C. T. Sachrajda, A. Soni, andJ. Yu (RBC, UKQCD), Phys. Rev.
D88 , 014508 (2013),arXiv:1212.5931 [hep-lat].[21] Z. Bai, N. H. Christ, T. Izubuchi, C. T. Sachrajda,A. Soni, and J. Yu, Phys. Rev. Lett. , 112003 (2014), arXiv:1406.0916 [hep-lat].[22] B. Wang, (2020), arXiv:2001.06374 [hep-lat].[23] Z. Bai,
Proceedings, 34th International Symposium onLattice Field Theory (Lattice 2016): Southampton, UK,July 24-30, 2016 , PoS
LATTICE2016 , 309 (2017),arXiv:1611.06601 [hep-lat].[24] N. H. Christ, X. Feng, A. Portelli, and C. T. Sachra-jda (RBC, UKQCD), Phys. Rev.
D92 , 094512 (2015),arXiv:1507.03094 [hep-lat].[25] N. H. Christ, X. Feng, A. Portelli, and C. T. Sachra-jda (RBC, UKQCD), Phys. Rev.
D93 , 114517 (2016),arXiv:1605.04442 [hep-lat].[26] N. H. Christ, X. Feng, A. J¨uttner, A. Lawson, A. Portelli,and C. T. Sachrajda, Phys. Rev.
D94 , 114516 (2016),arXiv:1608.07585 [hep-lat].[27] Z. Bai, N. H. Christ, X. Feng, A. Lawson, A. Portelli, andC. T. Sachrajda, Phys. Rev. Lett. , 252001 (2017),arXiv:1701.02858 [hep-lat].[28] Z. Bai, N. H. Christ, X. Feng, A. Lawson, A. Portelli,and C. T. Sachrajda, Phys. Rev.
D98 , 074509 (2018),arXiv:1806.11520 [hep-lat].[29] N. H. Christ, X. Feng, A. Portelli, and C. T. Sachra-jda (RBC, UKQCD), Phys. Rev.
D100 , 114506 (2019),arXiv:1910.10644 [hep-lat].[30] P. E. Shanahan, B. C. Tiburzi, M. L. Wagman, F. Win-ter, E. Chang, Z. Davoudi, W. Detmold, K. Orginos,and M. J. Savage, Phys. Rev. Lett. , 062003 (2017),arXiv:1701.03456 [hep-lat].[31] B. C. Tiburzi, M. L. Wagman, F. Winter, E. Chang,Z. Davoudi, W. Detmold, K. Orginos, M. J. Savage,and P. E. Shanahan, Phys. Rev.
D96 , 054505 (2017),arXiv:1702.02929 [hep-lat].[32] A. Nicholson et al. , Phys. Rev. Lett. , 172501 (2018),arXiv:1805.02634 [nucl-th].[33] X. Feng, L.-C. Jin, X.-Y. Tuo, and S.-C. Xia, Phys. Rev.Lett. , 022001 (2019), arXiv:1809.10511 [hep-lat].[34] W. Detmold and D. Murphy,
Proceedings, 36th Interna-tional Symposium on Lattice Field Theory (Lattice 2018):East Lansing, MI, United States, July 22-28, 2018 , PoS
LATTICE2018 , 262 (2019), arXiv:1811.05554 [hep-lat].[35] X.-Y. Tuo, X. Feng, and L.-C. Jin, Phys. Rev.
D100 ,094511 (2019), arXiv:1909.13525 [hep-lat].[36] C. Bouchard, C. C. Chang, T. Kurth, K. Orginos,and A. Walker-Loud, Phys. Rev.
D96 , 014504 (2017),arXiv:1612.06963 [hep-lat].[37] C.-Y. Seng and U.-G. Meißner, Phys. Rev. Lett. ,211802 (2019), arXiv:1903.07969 [hep-ph].[38] S. A. Larin and J. A. M. Vermaseren, Phys. Lett.
B259 ,345 (1991).[39] P. A. Baikov, K. G. Chetyrkin, and J. H. Kuhn, Phys.Rev. Lett. , 132004 (2010), arXiv:1001.3606 [hep-ph].[40] T. Blum et al. (RBC, UKQCD), Phys. Rev.
D93 , 074505(2016), arXiv:1411.7017 [hep-lat].[41] X. Feng, Y. Fu, and L.-C. Jin, (2019), arXiv:1911.04064[hep-lat].[42] K. G. Chetyrkin, J. H. Kuhn, and M. Steinhauser,Comput. Phys. Commun. , 43 (2000), arXiv:hep-ph/0004189 [hep-ph].[43] D. Pocanic et al. , Phys. Rev. Lett. , 181803 (2004),arXiv:hep-ex/0312030 [hep-ex].[44] A. Czarnecki, W. J. Marciano, and A. Sirlin, (2019),arXiv:1911.04685 [hep-ph].[45] W. Jaus, Phys. Rev. D63 , 053009 (2001), arXiv:hep-ph/0003015 [hep-ph]. [46] V. Cirigliano, M. Knecht, H. Neufeld, and H. Pichl, Eur.Phys. J.
C27 , 255 (2003), arXiv:hep-ph/0209226 [hep-ph].[47] M. Passera, K. Philippides, and A. Sirlin, Phys. Rev.
D84 , 094030 (2011), arXiv:1109.1069 [hep-ph].[48] J. Erler, Rev. Mex. Fis. , 200 (2004), arXiv:hep-ph/0211345 [hep-ph].[49] D. H. Wilkinson and B. E. F. Macefield, Nucl. Phys. A158 , 110 (1970).[50] C.-Y. Seng, D. Galviz, and U.-G. Meißner, JHEP ,069 (2020), arXiv:1910.13208 [hep-ph].[51] C.-Y. Seng, X. Feng, M. Gorchtein, and L.-C. Jin,(2020), arXiv:2003.11264 [hep-ph].[52] M. Tanabashi et al. (Particle Data Group), Phys. Rev. D98 , 030001 (2018).[53] B. Markisch et al. , Phys. Rev. Lett.122