F K / F π from Möbius domain-wall fermions solved on gradient-flowed HISQ ensembles
Nolan Miller, Henry Monge-Camacho, Chia Cheng Chang, Ben Hörz, Enrico Rinaldi, Dean Howarth, Evan Berkowitz, David A. Brantley, Arjun Singh Gambhir, Christopher Körber, Christopher J. Monahan, M. A. Clark, Bálint Joó, Thorsten Kurth, Amy Nicholson, Kostas Orginos, Pavlos Vranas, André Walker-Loud
LLLNL-JRNL-809712, RIKEN-iTHEMS-Report-20, JLAB-THY-20-3192 F K /F π from M¨obius Domain-Wall fermions solved on gradient-flowed HISQ ensembles Nolan Miller, Henry Monge-Camacho, Chia Cheng Chang ( 張 家 丞 ),
2, 3, 4
Ben H¨orz, Enrico Rinaldi,
5, 2
Dean Howarth,
6, 3
Evan Berkowitz,
7, 8
David A. Brantley, Arjun Singh Gambhir,
9, 3
Christopher K¨orber,
4, 3
Christopher J. Monahan,
10, 11
M.A. Clark, B´alint Jo´o, Thorsten Kurth, Amy Nicholson, Kostas Orginos,
10, 11
Pavlos Vranas,
6, 3 and Andr´e Walker-Loud
3, 6, 4 Department of Physics and Astronomy, University of North Carolina, Chapel Hill, NC 27516-3255, USA Interdisciplinary Theoretical and Mathematical Sciences Program(iTHEMS), RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA Department of Physics, University of California, Berkeley, CA 94720, USA Arithmer Inc., R&D Headquarters, Minato, Tokyo 106-6040, Japan Physics Division, Lawrence Livermore National Laboratory, Livermore, CA 94550, USA Department of Physics, University of Maryland, College Park, MD 20742, USA Institut f¨ur Kernphysik and Institute for Advanced Simulation, Forschungszentrum J¨ulich, 54245 J¨ulich Germany Design Physics Division, Lawrence Livermore National Laboratory, Livermore, CA 94550, USA Department of Physics, The College of William & Mary, Williamsburg, VA 23187, USA Theory Center, Thomas Jefferson National Accelerator Facility, Newport News, VA 23606, USA NVIDIA Corporation, 2701 San Tomas Expressway, Santa Clara, CA 95050, USA Scientific Computing Group, Thomas Jefferson National Accelerator Facility, Newport News, VA 23606, USA (Dated: May 19, 2020 - 1:17)We report the results of a lattice QCD calculation of F K /F π using M¨obius Domain-Wall fermionscomputed on gradient-flowed N f = 2 + 1 + 1 HISQ ensembles. The calculation is performedwith five values of the pion mass ranging from 130 (cid:46) m π (cid:46)
400 MeV, four lattice spacingsof a ∼ . , . , .
09 and 0 .
06 fm and multiple values of the lattice volume. The interpola-tion/extrapolation to the physical pion and kaon mass point, the continuum, and infinite volumelimits are performed with a variety of different extrapolation functions utilizing both the relevantmixed-action effective field theory expressions as well as discretization-enhanced continuum chiralperturbation theory formulas. We find that the a ∼ .
06 fm ensemble is helpful, but not necessaryto achieve a sub-percent determination of F K /F π . We also include an estimate of the strong isospinbreaking corrections and arrive at a final result of F ˆ K + /F ˆ π + = 1 . | V us | / | V ud | = 0 . CONTENTS
I. Introduction 1II. Details of the lattice calculation 2A. MDWF on gradient-flowed HISQ 2B. Correlation function construction andanalysis 41. Analysis strategy 4III. Extrapolation Functions 4A. N LO χ PT 5B. Discretization Corrections 8C. Finite Volume Corrections 9D. N LO Corrections 10IV. Extrapolation Details and Uncertainty Analysis 10A. Prior widths for LECs 10B. Physical Point 11C. Model Averaging Procedure 11D. Full analysis and uncertainty breakdown 121. Impact of a06m310L ensemble 122. Convergence of the chiral expansion 13E. QCD isospin breaking corrections 14 V. Summary and Discussion 17Acknowledgments 19A. Models included in final analysis 19B. NLO Mixed Action Formulas 19C. HMC for new ensembles 20References 23
I. INTRODUCTION
Leptonic decays of the charged pions and kaons pro-vide a means for probing flavor-changing interactions ofthe Standard Model (SM). In particular, the SM pre-dicts that the Cabibbo-Kobayashi-Maskawa (CKM) ma-trix is unitary, providing strict constraints on varioussums of the matrix elements. Thus, a violation of theseconstraints is indicative of new, beyond the SM (BSM)physics. There is a substantial flavor physics programdedicated to searching indirectly for potential violations. a r X i v : . [ h e p - l a t ] M a y CKM matrix elements may be determined through acombination of experimental leptonic decay widths andtheoretical determinations of the meson decay constants.For example, the ratio of the kaon and pion decay con-stants, F K , F π , respectively, may be related to the ratioof light and strange CKM matrix elements | V us | , | V ud | via [1, 2],Γ( K → l ¯ ν l )Γ( π → l ¯ ν l ) = | V us | | V ud | F K F π m K m π (cid:16) − m l m K (cid:17) (cid:16) − m l m π (cid:17) × (cid:2) δ EM + δ SU (2) (cid:3) . (1.1)In this expression, l = e, µ , the one-loop radiative QEDcorrection is δ EM [3, 4] and δ SU (2) is the strong isospinbreaking correction that relates F K /F π in the isospinlimit to F K + /F π + that includes m d − m u corrections [5] F K + F π + = F K F π (cid:2) δ SU (2) (cid:3) . Using lattice QCD calculations of the ratio of decayconstants in the above expression yields one of the mostprecise determinations of | V us | / | V ud | [6]. Combining theresults obtained through lattice QCD with independentdeterminations of the CKM matrix elements, such assemileptonic meson decays, provides a means for testingthe unitarity of the CKM matrix and obtaining signalsof new physics. F K /F π is a so-called gold-plated quantity [7] for calcu-lating within lattice QCD. This dimensionless ratio skirtsthe issue of determining a physical scale for the lattices,and gives precise results due to the correlated statisti-cal fluctuations between numerator and denominator, aswell as the lack of signal-to-noise issues associated withcalculations involving, for instance, nucleons. LatticeQCD calculations of F K /F π are now a mature endeavor,with state-of-the-art calculations determining this quan-tity consistently with sub-percent precision. The mostrecent review by the Flavour Lattice Averaging Group(FLAG), which performs global averages of quantitiesthat have been calculated and extrapolated to the phys-ical point by multiple groups, quotes a value of F ˆ K + F π + = 1 . N f = 2 + 1 + 1 dynamical quark flavors, includingstrong-isospin breaking corrections [8].This average includes calculations derived from twodifferent lattice actions, one [9] with twisted-massfermions [10, 11] and the other two [12, 13] with thehighly improved staggered quark (HISQ) action [14]. Theresults obtained using the HISQ action are approximatelyseven times more precise than those from twisted massand so the universality of the continuum limit for F K /F π from N f = 2 + 1 + 1 results has not been tested withprecision yet: in the continuum limit, all lattice actions should reduce to a single universal limit, that of SMQCD, provided all systematics are properly accountedfor. Thus, in addition to lending more confidence toits global average, the calculation of a gold-plated quan-tity also allows for precise testing of new lattice actions,and the demonstration of control over systematic un-certainties for a given action. FLAG also reports av-erages for N F = 2 + 1, F K ± /F π ± = 1 . N f = 2, F K ± /F π ± = 1 . N f = 2 + 1 + 1 results just for simplicity.In this work, we report a new determination of F K /F π calculated with M¨obius Domain-Wall fermions computedon gradient-flowed N f = 2 + 1 + 1 HISQ ensembles [22].Our final result in the isospin symmetric limit, Sec. IV D,including a breakdown in terms of statistical ( s ), pion andkaon mass extrapolation ( χ ), continuum limit ( a ), infi-nite volume limit ( V ), physical point (phys) and modelselection ( M ) uncertainties, is F K F π = 1 . s (12) χ (20) a (01) V (15) phys (12) M = 1 . . (1.3)With our estimated strong isospin breaking corrections,Sec. IV E, our result including m d − m u effects is F ˆ K + F ˆ π + = 1 . iso = 1 . , (1.4)where the first uncertainty in the first line is the combi-nation of those in Eq. (1.3).In the following sections we will discuss details of ourlattice calculation, including a brief synopsis of the actionand ensembles used, as well as our strategy for extractingthe relevant quantities from correlation functions. Wewill then detail our procedure for extrapolating to thephysical point via combined continuum, infinite volume,and physical pion and kaon mass limits and the resultinguncertainty breakdown. We discuss the impact of the a ∼ .
06 fm ensemble on our analysis, the convergence ofthe SU (3)-flavor chiral expansion, and the estimate of thestrong isospin breaking corrections. We conclude with anestimate of the impact our result has on improving theextraction of | V us | / | V ud | and an outlook. II. DETAILS OF THE LATTICE CALCULATIONA. MDWF on gradient-flowed HISQ
There are many choices for discretizing QCD, witheach choice being commonly referred to as a lattice ac-tion. These actions correspond to different UV theo-ries that share a common low-energy theory, QCD. Suf-ficiently close to the continuum limit, the discrete lat-tice actions can be expanded as a series of local oper-ators known as the Symanzik Expansion [26, 27], the
TABLE I. Input parameters for our lattice action. The abbreviated ensemble name [23] indicates the approximate latticespacing in fm and pion mass in MeV. The S, L, XL which come after an ensemble name denote a relatively small, large andextra-large volume with respect to m π L = 4. ensemble β N cfg volume am l am s am c L /a aM b , c am val l am res l × am val s am res s × σ N N src a15m400 a ×
48 0.0217 0.065 0.838 12 1.3 1.50, 0.50 0.0278 9.365(87) 0.0902 6.937(63) 3.0 30 8a15m350 a ×
48 0.0166 0.065 0.838 12 1.3 1.50, 0.50 0.0206 9.416(90) 0.0902 6.688(62) 3.0 30 16a15m310 5.80 1000 16 ×
48 0.013 0.065 0.838 12 1.3 1.50, 0.50 0.0158 9.563(67) 0.0902 6.640(44) 4.2 45 24a15m220 5.80 1000 24 ×
48 0.0064 0.064 0.828 16 1.3 1.75, 0.75 0.00712 5.736(38) 0.0902 3.890(25) 4.5 60 16a15m135XL a ×
64 0.002426 0.06730 0.8447 24 1.3 2.25, 1.25 0.00237 2.706(08) 0.0945 1.860(09) 3.0 30 32a12m400 a ×
64 0.0170 0.0509 0.635 8 1.2 1.25, 0.25 0.0219 7.337(50) 0.0693 5.129(35) 3.0 30 8a12m350 a ×
64 0.0130 0.0509 0.635 8 1.2 1.25, 0.25 0.0166 7.579(52) 0.0693 5.062(34) 3.0 30 8a12m310 6.00 1053 24 ×
64 0.0102 0.0509 0.635 8 1.2 1.25, 0.25 0.0126 7.702(52) 0.0693 4.950(35) 3.0 30 8a12m220S 6.00 1000 24 ×
64 0.00507 0.0507 0.628 12 1.2 1.50, 0.50 0.00600 3.990(42) 0.0693 2.390(24) 6.0 90 4a12m220 6.00 1000 32 ×
64 0.00507 0.0507 0.628 12 1.2 1.50, 0.50 0.00600 4.050(20) 0.0693 2.364(15) 6.0 90 4a12m220L 6.00 1000 40 ×
64 0.00507 0.0507 0.628 12 1.2 1.50, 0.50 0.00600 4.040(26) 0.0693 2.361(19) 6.0 90 4a12m130 6.00 1000 48 ×
64 0.00184 0.0507 0.628 20 1.2 2.00, 1.00 0.00195 1.642(09) 0.0693 0.945(08) 3.0 30 32a09m400 a ×
64 0.0124 0.037 0.44 6 1.1 1.25, 0.25 0.0160 2.532(23) 0.0491 1.957(17) 3.5 45 8a09m350 a ×
64 0.00945 0.037 0.44 6 1.1 1.25, 0.25 0.0121 2.560(24) 0.0491 1.899(16) 3.5 45 8a09m310 6.30 780 32 ×
96 0.0074 0.037 0.44 6 1.1 1.25, 0.25 0.00951 2.694(26) 0.0491 1.912(15) 6.7 167 8a09m220 6.30 1001 48 ×
96 0.00363 0.0363 0.43 8 1.1 1.25, 0.25 0.00449 1.659(13) 0.0491 0.834(07) 8.0 150 6a09m135 a ×
96 0.001326 0.03636 0.4313 12 1.1 1.50, 0.50 0.00152 0.938(06) 0.04735 0.418(04) 3.5 45 16a06m310L a ×
96 0.0048 0.024 0.286 6 1.0 1.25, 0.25 0.00617 0.225(03) 0.0309 0.165(02) 3.5 45 8 a Additional ensembles generated by CalLat using the MILC code. The m350 and m400 ensembles were made on the Vulcansupercomputer at LLNL while the a15m135XL, a09m135, and a06m310L ensembles were made on the Sierra and Lassensupercomputers at LLNL and the Summit supercomputer at OLCF using
QUDA [24, 25]. These configurations are available to anyintersted party upon request, and will be available for easy anonymous downloading—hopefully soon. low-energy Effective Field Theory (EFT) for the discretelattice action. The Symanzik EFT contains a series ofoperators having higher dimension than those in QCD,multiplied by appropriate powers of the lattice spacing, a . For all lattice actions, the only operators of mass-dimension ≤ a →
0, in that all lattice actions,if calculated using sufficiently small lattice spacing, willrecover the target theory of QCD, provided there are nosurprises from non-perturbative effects.Performing LQCD calculations with different actions istherefore valuable to test this universality, to help ensurea given action is not accidentally in a different phase ofQCD, and to protect against unknown systematic uncer-tainties arising from a particular calculation with a par-ticular action. In this work, we use a mixed-action [28] inwhich the discretization scheme for the valence-quarks isthe M¨obius Domain-Wall Fermion (MDWF) action [29–31] while the discretization scheme for the sea-quarks isthe Highly Improved Staggered Quark action [14]. Beforesolving the MDWF propagators, we apply a gradient-flow [32–34] smoothing algorithm [35, 36] to the glu-ons to dampen UV fluctuations, which also significantlyimproves the chiral symmetry properties of the MDWFaction [22] (for example, the residual chiral symmetrybreaking scale of domain-wall fermions m res is held toless than 10% of m l for reasonable values of L and M , see Tab. I). Our motivation to perform this calcu-lation is to improve our understanding of F K /F π and totest the MDWF on gradient-flowed HISQ action we haveused to compute the π − → π + neutrinoless double beta decay matrix elements arising from prospective higher-dimension lepton-number-violating physics [37], and theaxial coupling of the nucleon g A [38, 39]. As there isan otherwise straightforward path to determining g A tosub-percent precision with pre-exascale computing suchas Summit at OLCF and Lassen at LLNL [40], it is im-portant to ensure this action is consistent with knownresults at this level of precision.There are several motivations for chosing this mixed-action (MA) scheme [28, 41]. The MILC Collabora-tion provides their gauge configurations to any inter-ested party and we have made heavy use of them. Theyhave generated the configurations covering a large pa-rameter space allowing one to fully control the physicalpion mass, infinite volume and continuum limit extrap-olations [42, 43]. The good chiral symmetry propertiesof the Domain Wall (DW) action [44–46] significantlysuppress sources of chiral symmetry breaking from anysea-quark action, motivating the use of this mixed-actionsetup. While this action is not unitary at finite latticespacing, we have tuned the valence quark masses suchthat the valence pion mass matches the taste-5 HISQ pionmass within a few percent, so as the continuum limit istaken, we recover a unitary theory.EFT can be used to understand the salient features ofsuch MALQCD calculations. Chiral Perturbation The-ory ( χ PT) [47–49] can be extended to incorporate dis-cretization effects into the analytic formula describingthe quark-mass dependence of various hadronic quan-tities [50]. The MAEFT [51] for DW valence fermionson dynamical rooted staggered fermions is well devel-oped [52–59]. The use of valence fermions which re-spect chiral symmetry leads to a universal form of theMAEFT extrapolation formulae at NLO (next-to-leadingorder) in the joint quark mass and lattice spacing expan-sions [55, 58], which follows from the suppression of chiralsymmetry breaking discretization effects.
B. Correlation function construction and analysis
The correlation function construction and analysis fol-lows closely the strategy of Ref. [22] and [38, 39]. Herewe summarize the relevant details for this work.The pesudoscalar decay constants F can be obtainedfrom standard two-point correlation functions by makinguse of the 5D Ward-Takahashi Identity [60, 61] F q q = z q q p m q + m q res + m q + m q res ( E q q ) / , (2.1)where q and q denote the quark content of the mesonwith lattice input masses m q and m q respectively. Thepoint-sink ground-state overlap-factor z p and ground-state energy E are extracted from a two-point correla-tion function analysis with the model C q q ( ss ) ps ( t ) = (cid:88) n z q q n ( s ) p z q q † ns (cid:16) e − E q q n t + e − E q q n ( T − t ) (cid:17) , (2.2)where n encompasses in general an infinite tower ofstates, t is the source-sink time separation, T is thetemporal box size and we have both smeared ( s ) andpoint ( p ) correlation functions which both come fromsmeared sources. From Ref. [22], we show that gradient-flow smearing leads to the suppression of the domain-wall fermion oscillating mode (which also decouples as M →
1, at least in free-field [62]), and therefore thismode is not included in the correlator fit model. Finally,the residual chiral symmetry breaking m res is calculatedby the ratio of two-point correlation functions evaluatedat the midpoint of the fifth dimension L / m res ( t ) = (cid:80) x (cid:104) π ( t, x , L / π (0 , , (cid:105) (cid:80) x (cid:104) π ( t, x , π (0 , , (cid:105) , (2.3)where π ( t, x , s ) ≡ ¯ q ( t, x , w ) γ q ( t, x , w ) is the pseu-doscalar interpolating operator at time t , space x andfifth dimension s . We extract m res by fitting Eq. (2.3) toa constant.
1. Analysis strategy
For all two-point correlation function parameters(MDWF and mixed MDWF-HISQ), we infer posteriorparameter distributions in a Bayesian framework us-ing a 4-state model which simultaneously describes thesmeared- and point-sink two-point correlation functions(the source is always smeared). The joint posterior dis-tribution is approximated by a multivariate normal dis-tribution (we later refer to this procedure as fitting ). The two-point correlation functions are folded in time to dou-ble the statistics. The analysis of the pion, kaon, ¯ sγ s ,and mixed MDWF-HISQ mesons are performed indepen-dently, with correlations accounted for under bootstrapresampling.We analyze data of source-sink time separations be-tween 0.72 fm to 3.6 fm for all 0.09 fm and 0.12 fm latticespacing two-point correlation functions, and 0.75 fm to3.6 fm for all 0.15 fm lattice spacing two-point correlationfunctions.We choose normally distributed priors for the ground-state energy and all overlap factors, and log-normal dis-tributions for excited-state energy priors. The ground-state energy and overlap factors are motivated by theplateau values of the effective masses and scaled corre-lation function, and a prior width of 10% of the centralvalue. The excited-state energy splittings are set to thevalue of two pion masses with a width allowing for fluc-tuations down to one pion mass within one standard de-viation. The excited-state overlap factors are set to zero,with a width set to the mean value of the ground-stateoverlap factor.Additionally, we fit a constant to the correlation func-tions in Eq. (2.3). For the 0.09 fm to 0.12 fm ensembles,we analyze source-sink separations that are greater than0.72 fm. For the 0.12 fm ensemble, the minimum source-sink separation is 0.75 fm. The prior distribution for theresidual chiral symmetry breaking parameter is set to theobserved value per ensemble, with a width that is 100%of the central value. The uncertainty is propagated withbootstrap resampling.We emphasize that all input fit parameters ( i.e. num-ber of states, fit region, priors) are chosen to have thesame values in physical units for all observables, to theextent that a discretized lattice allows. Additionally, wenote that the extracted ground-state observables fromthese correlation functions are insensitive to variationsaround the chosen set of input fit parameters. III. EXTRAPOLATION FUNCTIONS
We now turn to the extrapolation/interpolation to thephysical point. We have three ensembles at the physi-cal pion mass with relatively high statistics, a15m135XL,a12m130, and a09m135 with precise determinations of F K /F π , see Tab. II, such that the physical quark massextrapolation is an interpolation. Nevertheless, we ex-plore how the ensembles with heavier pion masses impactthe physical point prediction and we use our dataset toexplore uncertainty arising in the SU (3)-flavor chiral ex-pansion.We begin by assuming a canonical power countingscheme for our MA LQCD action [52] in which O( m π ) ∼ O( m K ) ∼ O( a Λ ) are all treated as small scales. Forthe quark mass expansion, the dimensionless small pa-rameters ( m P / πF ) naturally emerge from χ PT where P ∈ { π, K, η } . For the discretization corrections, while TABLE II. Extracted masses and decay constants from correlation functions. An HDF5 file is provided with this publicationwhich includes the resulting bootstrap samples which can be used to construct the correlated uncertainties. The small param-eters in this table are defined as (cid:15) π,K = m π,K / (4 πF π ), (cid:15) a = a/ (2 w ). The normalization of (cid:15) a is chosen such that as a smallparameter, it spans the range of (cid:15) π (cid:46) (cid:15) a (cid:46) (cid:15) K .ensemble am π am K (cid:15) π (cid:15) K m π L (cid:15) a α S aF π aF K F K /F π a15m400 0.30281(31) 0.42723(27) 0.09216(33) 0.18344(62) 4.85 0.19378(13) 0.58801 0.07938(12) 0.08504(09) 1.0713(09)a15m350 0.26473(30) 0.41369(28) 0.07505(28) 0.18326(60) 4.24 0.19378(13) 0.58801 0.07690(11) 0.08370(09) 1.0884(09)a15m310 0.23601(29) 0.40457(25) 0.06223(17) 0.18285(48) 3.78 0.19378(13) 0.58801 0.07529(09) 0.08293(09) 1.1015(13)a15m220 0.16533(19) 0.38690(21) 0.03269(11) 0.17901(48) 3.97 0.19378(13) 0.58801 0.07277(08) 0.08196(10) 1.1263(15)a15m135XL 0.10293(07) 0.38755(14) 0.01319(05) 0.18704(59) 4.94 0.19378(13) 0.58801 0.07131(11) 0.08276(10) 1.1606(18)a12m400 0.24347(16) 0.34341(14) 0.08889(30) 0.17685(63) 5.84 0.12376(18) 0.53796 0.06498(11) 0.06979(07) 1.0739(17)a12m350 0.21397(20) 0.33306(16) 0.07307(37) 0.17704(83) 5.14 0.12376(18) 0.53796 0.06299(14) 0.06851(07) 1.0876(27)a12m310 0.18870(17) 0.32414(21) 0.05984(25) 0.17657(69) 4.53 0.12376(18) 0.53796 0.06138(11) 0.06773(10) 1.1033(21)a12m220S 0.13557(32) 0.31043(22) 0.03384(19) 0.1774(10) 3.25 0.12376(18) 0.53796 0.05865(16) 0.06673(11) 1.1378(27)a12m220L 0.13402(15) 0.31021(19) 0.03289(15) 0.17621(79) 5.36 0.12376(18) 0.53796 0.05881(13) 0.06631(17) 1.1276(29)a12m220 0.13428(17) 0.31001(17) 0.03314(15) 0.17666(81) 4.30 0.12376(18) 0.53796 0.05870(13) 0.06636(11) 1.1306(22)a12m130 0.08126(16) 0.30215(11) 0.01287(08) 0.17788(71) 3.90 0.12376(18) 0.53796 0.05701(11) 0.06624(08) 1.1619(21)a09m400 0.18116(15) 0.25523(13) 0.08883(32) 0.17633(59) 5.80 0.06515(08) 0.43356 0.04837(08) 0.05229(07) 1.0810(09)a09m350 0.15785(20) 0.24696(12) 0.07256(32) 0.17761(68) 5.05 0.06515(08) 0.43356 0.04663(08) 0.05127(07) 1.0994(10)a09m310 0.14072(12) 0.24106(14) 0.06051(22) 0.17757(59) 4.50 0.06515(08) 0.43356 0.04552(07) 0.05053(08) 1.1101(16)a09m220 0.09790(06) 0.22870(09) 0.03307(14) 0.18045(70) 4.70 0.06515(08) 0.43356 0.04284(08) 0.04899(07) 1.1434(18)a09m135 0.05946(06) 0.21850(08) 0.01346(08) 0.18175(91) 3.81 0.06515(08) 0.43356 0.04079(10) 0.04804(06) 1.1778(22)a06m310L 0.09456(06) 0.16205(07) 0.06141(35) 0.1803(10) 6.81 0.02726(03) 0.29985 0.03037(08) 0.03403(07) 1.1205(17) a Λ is often used to estimate the relative size of cor-rections compared to typical hadronic mass scales, it is abit unatural to use this in a low-energy EFT as Λ QCD isa QCD scale that does not emerge in χ PT.We chose to use another hadronic scale to form a di-mensionless parameter with a , that being the gradientflow scale w ∼ .
17 fm [63]. This quantity is easy tocompute, has mild quark mass dependence, and the valueis roughly w − (cid:39) πF π . We then define the dimension-less small parameters for controlling the expansion to be (cid:15) P = (cid:18) m P Λ χ (cid:19) = (cid:16) m P πF (cid:17) , (cid:15) a = (cid:18) a w (cid:19) . (3.1)We leave F ambiguous, as we will explore taking F = F π , F = F K and F = F π F K in our definition of Λ χ . Thisparticular choice of (cid:15) a is chosen such that the rangeof values of this small parameter roughly correspond to (cid:15) π (cid:46) (cid:15) a (cid:46) (cid:15) K as the lattice spacing is varied, similar tothe variation of (cid:15) π itself over the range of pion massesused, see Tab. II. As we will discuss in Sec. IV, thischoice of (cid:15) a seems natural as determined by the size ofthe discretization LECs which are found in the analysis.Note, this differs from the choice used in our analysis of g A [38, 39].With this power counting scheme, the different ordersin the expansion are defined to beNLO: O( (cid:15) P ) ∼ O( (cid:15) a ) , N LO: O( (cid:15) P ) ∼ O( (cid:15) P (cid:15) a ) ∼ O( (cid:15) a ) , N LO: O( (cid:15) P ) ∼ O( (cid:15) P (cid:15) a ) ∼ O( (cid:15) P (cid:15) a ) ∼ O( (cid:15) a ) . (3.2)Even at finite lattice spacing, F K = F π in the SU (3)flavor symmetry limit, also known as the SU (3) vector limit SU (3) V , and so there can not be a pure O( (cid:15) a ) cor-rection as it must accompany terms which vanish in the SU (3) V limit, such as (cid:15) K − (cid:15) π . Therefore, at NLO, therecan not be any counterterms proportional to (cid:15) a and theonly discretization effects that can appear at NLO comethrough modification of the various meson masses thatappear in the MA EFT.We find that the precision of our results requires in-cluding terms higher than NLO, and we have to workat a hybrid N LO order to obtain a good description ofour data. Therefore, we will begin with a discussion ofthe full N LO χ PT theory expression for F K /F π in thecontinuum limit [64–67]. A. N LO χ PT The analytic expression for F K /F π up to N LO is [67] F K F π = 1 + 58 (cid:96) π − (cid:96) k − (cid:96) η + 4 ¯ L ( (cid:15) K − (cid:15) π )+ (cid:15) K F F (cid:18) m π m K (cid:19) + ˆ K r λ π + ˆ K r λ π λ K + ˆ K r λ π λ η + ˆ K r λ K + ˆ K r λ K λ η + ˆ K r λ η + ˆ C r λ π + ˆ C r λ K + ˆ C r λ η + ˆ C r . (3.3)The first line is the LO (1) plus NLO terms, while thenext three lines are the N LO terms. Several non-uniquechoices were made to arrive at this formula. Prior todiscussing these choices, we first define the parametersappearing in Eq. (3.3). First, the small parameters wereall defined as (cid:15) P = (cid:18) m P πF π ( m P ) (cid:19) , (3.4)where F π ( m P ) is the “on-shell” pion decay constant atthe masses m P . The quantities (cid:96) P are defined as (cid:96) P = (cid:15) P ln (cid:18) m P µ (cid:19) , (3.5)where µ is a renormalization scale. The coefficient ¯ L =(4 π ) L r ( µ ) is one of the regulated Gasser-LeutwylerLECs [68] which has a renormalization scale dependencethat exactly cancels against the dependence arising fromthe logarithms appearing at the same order. In the fol-lowing, we define all of the Gasser-Leutwyler LECs withthe extra (4 π ) for convenience:¯ L i ≡ (4 π ) L ri ( µ ) . (3.6)The η mass has been defined through the Gell-Mann–Okubo (GMO) relation m η ≡ m K − m π , (3.7)with the corrections to this relation being propagatedinto Eq. (3.3) for consistency at N LO. The logs are λ P ≡ ln (cid:18) m P µ (cid:19) . (3.8)The ln terms are encapsulated in the F F ( x ) function,defined in Eqs. (8-17) of Ref. [67], and the ˆ K ri λ P λ P (cid:48) terms whose coefficients are given by ˆ K r = 1124 (cid:15) π (cid:15) K − (cid:15) π , ˆ K r = − (cid:15) π (cid:15) K − (cid:15) π , ˆ K r = 1324 (cid:15) π (cid:15) K + 5996 (cid:15) π , ˆ K r = 1736 (cid:15) K + 7144 (cid:15) π (cid:15) K , ˆ K r = − (cid:15) K − (cid:15) π (cid:15) K + 332 (cid:15) π , ˆ K r = 241288 (cid:15) K − (cid:15) π (cid:15) K − (cid:15) π . (3.9)The single log coefficients ˆ C r − are combinations of theNLO Gasser-Leutwyler coefficientsˆ C ri = c ππi (cid:15) π + c Kπi (cid:15) K (cid:15) π + c KKi (cid:15) K , (3.10) They also provide an approximate formula which is easy to imple-ment, but our numerical results are sufficiently precise to requirethe exact expression. To implement this function in our analysis,we have modified an interface
C++ file provided by J. Bijnens to
CHIRON [69], the package for two loop χ PT functions. We haveprovided a
Python interface as well so that the function can becalled from our main analysis code, which is provided with thisarticle. We correct a typographical error in the K r term presented inRef. [67]: a simple power-counting reveals the ξ K = (cid:15) K accom-panying this term should not be there. where c ππ = − − L + 5 ¯ L ) −
132 ¯ L + 212 ¯ L ,c Kπ = − −
112 ¯ L ,c KK = c ππ = 0 ,c Kπ = 209144 + 3 ¯ L ,c KK = 5396 + 2(2 ¯ L + 5 ¯ L ) + 5 ¯ L − L ,c ππ = 19288 + 16 ¯ L + 116 ¯ L − L + ¯ L ) ,c Kπ = − −
43 ¯ L −
256 ¯ L + 16(2 ¯ L + ¯ L ) ,c KK = 1318 + 83 ¯ L −
23 ¯ L − L + ¯ L ) . (3.11)Finally, ˆ C r is a combination of these L ri coefficientsas well as counterterms appearing at N LO. At N LO,only two counterterm structures can appear due to the SU (3) V constraints:ˆ C r = ( (cid:15) K − (cid:15) π ) (cid:2) ( A K + L K ) (cid:15) K + ( A π + L π ) (cid:15) π (cid:3) (3.12)which are linear combinations of the N LO counterterms A K = 16(4 π ) ( C r + C r ) ,A π = 8(4 π ) ( C r + 2 C r ) , (3.13)and contributions from the Gasser-Leutwyler LECs(Eq. (7) of Ref. [67]) L K = 8 ¯ L (8( ¯ L − L ) + 3 ¯ L − L ) − L − ¯ L −
118 ¯ L + 43 ¯ L − L + ¯ L ) ,L π = 8 ¯ L (4( ¯ L − L ) + 5 ¯ L − L ) − L − ¯ L −
518 ¯ L −
43 ¯ L + 8(2 ¯ L + ¯ L ) . (3.14)There were several non-unique choices that went intothe determination of Eq. (3.3). When working with thefull N LO χ PT expression, the different choices one canmake result in different N LO or higher corrections andexploring these different choices in the analysis will ex-pose sensitivity to higher-order contributions that are notexplicitly included. The first choice we discuss is the Tay-lor expansion of the ratio of F K /F π F K F π = 1 + δF NLO K + δF N LO K + · · · δF NLO π + δF N LO π + · · · = 1 + δF NLO K − π + δF N LO K − π + (cid:0) δF NLO π (cid:1) − δF NLO π δF NLO K + · · · , (3.15)Where the · · · represent higher order terms in the ex-pansion and δF N LO K − π = δF N LO K − δF N LO π . Eq. (3.3) hasbeen derived from this standard Taylor-expanded formwith the choices mentioned above: the use of the onshellrenormalized value of F → F π and the definition of the η mass through the GMO relation. The NLO expressionsare the standard ones [68] δF NLO K = − (cid:96) π − (cid:96) K − (cid:96) η + 4 ¯ L (cid:15) K + 4 ¯ L ( (cid:15) π + 2 (cid:15) K ) ,δF NLO π = − (cid:96) π − (cid:96) K + 4 ¯ L (cid:15) π + 4 ¯ L ( (cid:15) π + 2 (cid:15) K ) ,δF NLO K − π = 58 (cid:96) π − (cid:96) K − (cid:96) η + 4 ¯ L ( (cid:15) K − (cid:15) π ) . (3.16)The δF N LO P terms have been determined in Ref. [64] andcast into analytic forms in Refs. [65, 66]. The NLO termsare of O(20%) and so Taylor expanding this ratio leads tosizeable corrections from the (cid:0) δF NLO π (cid:1) − δF NLO π δF NLO K contributions. Utilizing the full ratio expression could inprinciple lead to a noticeable difference in the analysis(a different determination of the values of the LECs forexample). Rather than implementing the full δF N LO P expressions for kaon and pion, we explore this conver-gence by instead just resumming the NLO terms whichwill dominate the potential differences in higher ordercorrections. A consistent expression at N LO is F K F π [(3.3)] = 1 + δF NLO K δF NLO π + δF N LO + δ N LOratio + · · · , (3.17)where δF N LO is the full N LO expression in Eq. (3.15) δF N LO = δF N LO K − π + (cid:0) δF NLO π (cid:1) − δF NLO π δF NLO K , (3.18)and the ratio correction is given by δ N LOratio = δF NLO π δF NLO K − (cid:0) δF NLO π (cid:1) . (3.19)Another choice we explore is the use of F → F π in thedefinition of the small parameters. Such a choice is veryconvenient as it allows one to express the small parame-ters entirely in terms of observables one can determinein the lattice calculation (unlike the bare parameters,such as χ PT’s F and Bm q , which must be determinedthrough extrapolation analysis). Equally valid, one couldhave chosen F → F K or F → F π F K . Each choice in-duces explicit corrections one must account for at N LOto have a consistent expression at this order. The NLOcorrections in Eq. (3.3) are proportional to (cid:15) P = m P (4 πF π ) = m P (4 π ) F π F K F K F π = m P (4 π ) F π F K (cid:0) δF NLO K − π (cid:1) = m P (4 πF K ) F K F π = m P (4 πF K ) (cid:0) δF NLO K − π (cid:1) , (3.20)plus higher order corrections. Related to this choice, Eq. (3.3) is implicitly definedat the standard renormalization scale [67] µ ρ = m ρ = 770 MeV . (3.21)While F K /F π of course does not depend upon this choice,the numerical values of the LECs do. Further, a scale set-ting would be required to utilize this or any fixed value of µ . Instead, as was first advocated in Ref. [70] to the bestof our knowledge, it is more convenient to set the renor-malization scale on each ensemble with a lattice quantity.For example, Ref. [70] used µ = f latt π = √ F latt π where F latt π is the lattice-determined value of the pion decayconstant on a given ensemble. The advantage of thischoice is that the entire extrapolation can be expressedin terms of ratios of lattice quantities such that a scalesetting is not required to perform the extrapolation tothe physical point.At NLO in the expansion, one is free to make thischoice as the corrections appear at N LO. In the presentwork, we must account for these corrections for a consis-tent expression at this order, which is still defined at afixed renormalization scale. To understand these correc-tions, we take as our fixed scale µ = 4 πF , (3.22)where F is the decay constant in the SU (3) chiral limit.Define µ π = 4 πF π and consider the NLO expression F K F π = 1 + 58 (cid:96) µ π − (cid:96) µ K − (cid:96) µ η + 4( (cid:15) K − (cid:15) π ) ¯ L ( µ )= + 58 (cid:15) π ln (cid:18) (cid:15) π µ π µ (cid:19) − (cid:15) K ln (cid:18) (cid:15) K µ π µ (cid:19) − (cid:15) η ln (cid:18) (cid:15) η µ π µ (cid:19) + 4( (cid:15) K − (cid:15) π ) ¯ L ( µ )= 1 + 58 (cid:96) µ π π − (cid:96) µ π K − (cid:96) µ π η + 4( (cid:15) K − (cid:15) π ) ¯ L ( µ )+ ln (cid:18) µ π µ (cid:19) (cid:20) (cid:15) π − (cid:15) K − (cid:15) η (cid:21) , (3.23)where we have introduced the notation (cid:96) µP = (cid:15) P ln (cid:18) (cid:15) P µ (cid:19) . (3.24)If we chose the renormalization scale µ π and add thesecond term of the last equality, then this expression isequivalent to working with the scale µ through N LO.The convenience of this choice becomes clear as µ π /µ has a familiar expansion µ π µ = 1 + δF NLO π + · · · . (3.25)Using the GMO relation Eq. (3.7) and expanding ln(1+ x )for small x , this expression becomes F K F π = 1 + 58 (cid:96) µ π π − (cid:96) µ π K − (cid:96) µ π η + 4( (cid:15) K − (cid:15) π ) ¯ L ( µ ) −
32 ( (cid:15) K − (cid:15) π ) δF NLO π . (3.26)Similar expressions can be derived for the choices µ πK =4 πF πK (where F πK = √ F π F K ) and µ K = 4 πF K whichare made more convenient if one also makes the re-placements F π → { F π F K , F K } in the definition of thesmall parameters plus the corresponding N LO correc-tions that accompany these choices.If we temporarily expose the implicit dependence ofthe expression for F K /F π on the choices of F and µ ,such that Eq. (3.3) is defined as F K F π [(3.3)] = F K F π ( F π , µ ρ ) , (3.27)then the following expressions are all equivalent at N LO F K F π ( F π , µ ) = F K F π ( F K , µ K ) + δ N LO F K = F K F π ( F πK , µ πK ) + δ N LO F πK = F K F π ( F π , µ π ) + δ N LO F π , (3.28)where δ N LO F K = −
32 ( (cid:15) K − (cid:15) π ) δF NLO K + 2( δF NLO K − π ) δ N LO F πK = −
34 ( (cid:15) K − (cid:15) π )( δF NLO K + δF NLO π ) + ( δF NLO K − π ) δ N LO F π = −
32 ( (cid:15) K − (cid:15) π ) δF NLO π (3.29)and the LECs in these expressions are related to thoseat the standard scale by evolving them from µ ρ → µ with their known scale dependence [68]. Implicit in theseexpressions is the normalization of the small parameters (cid:15) P = m P (4 πF π ) , for F → F πm P (4 π ) F π F K , for F → √ F π F Km P (4 πF K ) , for F → F K . (3.30)We have described several choices one can make in pa-rameterizing the χ PT formula for F K /F π . The key pointis that if the underlying chiral expansion is well behaved,the formula resulting from each choice are all equivalentthrough N LO in the SU (3) chiral expansion, with dif-ferences only appearing at N LO and beyond. Therefore,by studying the variance in the extrapolated answer uponthese choices, one is assessing some of the uncertaintyarising from the truncation of the chiral extrapolationformula.
B. Discretization Corrections
We now turn to the discretization corrections. We ex-plore two parametrizations for incorporating the correc-tions arising at finite lattice spacing. The simplest ap-proach is to use the continuum extrapolation formula and enhance it by adding contributions from all allowed pow-ers of (cid:15) P and (cid:15) a to a given order in the expansion. Thisis very similar to including only the contributions fromlocal counter-terms that appear at the given order. AtN LO, the set of discretization corrections is given by δ N LO a = A s (cid:15) a ( (cid:15) K − (cid:15) π ) + A α S α S (cid:15) a ( (cid:15) K − (cid:15) π ) , (3.31)where A s and A α S are the LECs and α S is the runningQCD coupling that emerges in the Symanzik expansion ofthe lattice expansion through loop corrections. Each con-tribution at this order must vanish in the SU (3) V limitbecause the discretization corrections are flavor-blind andso we have the limiting constraintlim m l → m s F K F π = 1 , (3.32)at any lattice spacing.From a purist EFT perspective, we should instead uti-lize the MA EFT expression. Unfortunately, the MAEFT expression is only known at NLO [52] and our re-sults require higher orders to provide good fits. Never-theless, we can explore the utility of the MA EFT byreplacing the NLO χ PT expression with the NLO MAEFT expression and using the continuum expression en-hanced with the local discretization corrections at higherorders, Eq. (3.31).Using the parameterization of the hairpin contribu-tions from Ref. [55], the NLO MA EFT expressions are δF MA π = − (cid:96) ju − (cid:96) ru L (cid:15) π + 4 ¯ L ( (cid:15) π + 2 (cid:15) K ) + (cid:15) a ¯ L a ,δF MA K = − (cid:96) ju (cid:96) π − (cid:96) ru − (cid:96) js − (cid:96) rs (cid:96) ss − (cid:96) X
8+ 4 ¯ L (cid:15) K + 4 ¯ L ( (cid:15) π + 2 (cid:15) K ) + (cid:15) a ¯ L a − δ ju d(cid:96) π − K πX ) − δ ju K (2 , πX + δ rs (cid:18) K ssX −
23 ( (cid:15) K − (cid:15) π ) K (2 , ssX (cid:19) + δ ju δ rs (cid:16) K (2 , ssX − K πssX (cid:17) . (3.33)In these expressions, we use the partially quenched flavornotation [71] in which π : valence-valence pion K : valence-valence kaon u : valence light flavor quark j : sea light flavor quark s : valence strange flavor quark r : sea strange flavor quark X : sea-sea eta meson , (3.34)and so, for example (cid:96) ju = m ju (4 πF π ) ln (cid:32) m ju (4 πF π ) (cid:33) , (3.35) TABLE III. Extracted masses of the mixed MDWF-HISQ mesons. We use the notation from Ref. [71] in which m π and m K denote the masses of the valence pion and kaon and j and r denote the light and strange flavors of the sea quarks while u and s denote the light and strange flavors of the valence quarks. Since we have tuned the valence MDWF pion and ¯ ss mesons to have the same mass as the HISQ sea pion and ¯ ss mesons within a few percent, the quantities m ju − m π andother splittings provide an estimate of the additive mixed-meson mass splitting due to discretization effects, a ∆ Mix [52] andadditional additive corrections [59]. At LO in MAEFT, these splittings are predicted to be quark mass independent, which wefind to be approximately true, with a notable decrease in the splitting as the valence-quark mass is increased as first observedin Ref. [56] as well as a milder decrease as the seq-quark mass is increased.ensemble am ju am js am ru am rs am ss w ∆ , ju w ∆ , js w ∆ , ru w ∆ , rs w a ∆ I a15m400 0.3597(17) 0.4586(24) 0.4717(19) 0.5537(11) 0.5219(02) 0.0486(15) 0.0359(28) 0.0516(23) 0.0440(16) 0.112(14)a15m350 0.3308(23) 0.4463(14) 0.4598(16) 0.5526(10) 0.5201(02) 0.0508(20) 0.0362(17) 0.0519(19) 0.0451(15) 0.112(14)a15m310 0.3060(17) 0.4345(16) 0.4508(14) 0.5490(12) 0.5188(02) 0.0489(13) 0.0324(18) 0.0511(17) 0.0416(16) 0.112(14)a15m220 0.2564(27) 0.4115(17) 0.4320(29) 0.5420(08) 0.5150(01) 0.0495(18) 0.0253(19) 0.0476(33) 0.0368(11) 0.112(14)a15m135XL 0.232(13) 0.4058(56) 0.4337(84) 0.5560(31) 0.5257(02) 0.0559(75) 0.0187(59) 0.0489(94) 0.0423(45) 0.112(14)a12m400 0.2678(06) 0.3560(08) 0.3624(07) 0.4333(06) 0.4207(01) 0.0251(07) 0.0177(12) 0.0271(10) 0.0217(11) 0.063(05)a12m350 0.2303(08) 0.3446(07) 0.3454(10) 0.4322(05) 0.4197(01) 0.0147(07) 0.0158(10) 0.0168(15) 0.0214(09) 0.063(05)a12m310 0.2189(09) 0.3344(10) 0.3439(09) 0.4305(05) 0.4180(02) 0.0248(08) 0.0136(14) 0.0266(13) 0.0213(09) 0.063(05)a12m220S 0.1774(14) 0.3187(12) 0.3323(17) 0.4286(10) 0.4158(02) 0.0264(10) 0.0105(16) 0.0283(24) 0.0219(18) 0.063(05)a12m220L 0.1774(14) 0.3187(12) 0.3323(17) 0.4286(10) 0.4156(02) 0.0273(10) 0.0107(16) 0.0286(23) 0.0222(18) 0.063(05)a12m220 0.1774(14) 0.3187(12) 0.3323(17) 0.4286(10) 0.4154(01) 0.0272(10) 0.0110(16) 0.0289(23) 0.0225(18) 0.063(05)a12m130 0.1491(20) 0.3080(15) 0.3240(26) 0.4271(08) 0.4141(01) 0.0316(12) 0.0073(19) 0.0276(34) 0.0220(14) 0.063(05)a09m400 0.1878(05) 0.2581(06) 0.2607(06) 0.3162(05) 0.3133(01) 0.0094(07) 0.0056(12) 0.0109(11) 0.0071(12) 0.020(02)a09m350 0.1654(06) 0.2498(05) 0.2526(06) 0.3159(04) 0.3124(01) 0.0093(07) 0.0054(10) 0.0108(12) 0.0083(11) 0.020(02)a09m310 0.1485(06) 0.2428(05) 0.2472(10) 0.3150(04) 0.3117(01) 0.0086(07) 0.0032(10) 0.0114(20) 0.0080(09) 0.020(02)a09m220 0.1090(09) 0.2303(06) 0.2334(07) 0.3115(03) 0.3094(01) 0.0088(07) 0.0028(10) 0.0083(12) 0.0051(08) 0.020(02)a09m135 0.0786(15) 0.2187(11) 0.2270(15) 0.3079(05) 0.3027(07) 0.0102(09) 0.0004(19) 0.0146(26) 0.0123(19) 0.020(02)a06m310L 0.0957(08) 0.1619(11) 0.1619(12) 0.2103(10) 0.2098(01) 0.0020(14) -0.0004(34) -0.0004(34) 0.0020(40) 0.004(00) where m ju is the mass of a mixed valence-sea pion.The partial quenching parameters δ ju and δ rs providea measure of the unitarity violation in the theory. Forour MDWF on HISQ action, at LO in MA EFT, theyare given by the splitting in the quark masses plus adiscretization correction arising from the taste-Identitysplitting δ ju = 2 B ( m j − m u ) + a ∆ I (4 πF π ) ,δ rs = 2 B ( m r − m s ) + a ∆ I (4 πF π ) . (3.36)For the tuning we have done, setting the valence-valencepion mass equal to the taste-5 sea-sea pion mass, theseparameters are given just by the discretization terms as m u = m j and m s = m r within 1-2%. The sea-sea etamass in this tuning is given at LO in MA EFT as m X = 43 m K − m π + a ∆ I . (3.37)These parameters, and the corresponding meson massesare provided in Tab. III. The expressions for d(cid:96) π , K φ φ , K (2 , φ φ and K φ φ φ are provided in Appendix B.At NLO in the MA EFT, the LECs which contributeto δF K and δF π are the same as in the continuum, L and L , plus a discretization LEC which we have denoted¯ L a . Just like the L contribution, the contribution from¯ L a exactly cancels in δF K − δF π . At N LO, beyond the
TABLE IV. Multiplicity factors for the finite volume correc-tions of the first 10 vector lengths, | n | . | n | √ √ √ √ √ √ √ √ √ c n continuum counterterm contributions, (3.12), there arethe two additional LECs contributions, (3.31). C. Finite Volume Corrections
We now discuss the corrections arising from the finitespatial volume. The leading FV corrections arise fromthe tadpole integrals which arise at NLO in both the χ PT and MA expressions. The well known modificationto the integral can be expressed as [72–74] (cid:96) µ π , FV P = (cid:96) µ π P + 4 (cid:15) P (cid:88) | n |(cid:54) =0 c n m P L | n | K ( m P L | n | ) , (3.38)where the sum runs over all non-zero integer three-vectors. Each value of | n | can be thought of as a windingof the meson P around the finite universe. The c n aremultiplicity factors counting all the ways to form a vec-tor of length | n | from triplets of integers, see Tab. IV forthe first few. K ( x ) is a modified Bessel Function of thesecond kind. In the asymptotically large volume limit,0 .
000 0 .
001 0 .
002 0 .
003 0 .
004 0 .
005 0 .
006 0 . e − m π L / ( m π L ) / . . . . . . . . F K / F π a m : δ N L O χ P T F V ( (cid:15) π , m π L ) m π L = 5 . m π L = 4 . m π L = 3 . FIG. 1. We compare the finite volume results on a12m220L,a12m220 and a12m220S to the predicted finite volume cor-rections from NLO χ PT. The uncertainty band is from thefull N LO χ PT extrapolation, plotted with fixed mesonsmasses ( (cid:15) P ) and fixed lattice spacing ( (cid:15) a ), determined fromthe a12m220L ensemble. At the one-sigma level, our data areconsistent with the leading FV corrections. the finite volume correction to these integrals is δ FV (cid:96) P ≡ (cid:96) FV P − (cid:96) P = (cid:15) P √ π e − m P L ( m P L ) / + (cid:15) P × O (cid:32) e − m P L √ ( m P L √ / , e − m P L ( m P L ) / (cid:33) . (3.39)The full finite volume corrections to the continuum for-mula are also known at N LO [75] as well as in the par-tially quenched χ PT [76]. In this work, we restrict thecorrections to those arising from the NLO corrections asour results are not sensitive to higher-order FV correc-tions. This is because, with the ensembles used in thiswork, all ensembles except a12m220S satisfy m π L (cid:38) χ PT. The uncertaintyband arises from an N LO fit using the full N LO con-tinuum χ PT formula enhanced with discretization LECsand N LO corrections arising from continuum and finitelattice spacing corrections. Even with one of the mostprecise fits, we see that the numerical results are consis-tent with the predicted NLO FV corrections.
D. N LO Corrections
The numerical data set in this work requires us to addN LO corrections to obtain a good fit quality. At thisorder, we only consider local counterterm contributions,of which there are three new continuum like correctionsand three discretization corrections. A non-unique, but complete parameterization is δ N LO = ( (cid:15) K − (cid:15) π ) (cid:26) (cid:15) a A s + (cid:15) a ( A s,K (cid:15) K + A s,π (cid:15) π )+ A Kπ (cid:15) K (cid:15) π + ( (cid:15) K − (cid:15) π )( A K (cid:15) K + A π (cid:15) π ) (cid:27) . (3.40)In principle, we could also add counterterms proportionalto higher powers of α S but with four lattice spacings, wewould not be able to resolve the difference between thecomplete set of operators including all possible additional α S corrections. The set of operators we do include issufficient to parametrize the approach to the continuumlimit. IV. EXTRAPOLATION DETAILS ANDUNCERTAINTY ANALYSIS
We now carry out the extrapolation/interpolation tothe physical point, which we perform in a Bayesian frame-work. To obtain a good fit, we must work to N LO in themixed chiral and continuum expansion. The results fromthe a06m310L ensemble drive this need, in particular,for higher-order discretization corrections to parametrizethe results from all the ensembles. We will explore theimpact of the a06m310L ensemble in more detail in thissection. First, we discuss the values of the priors we setand the definition of the physical point.
A. Prior widths for LECs
The number of additional LECs we need to determineat each order in the expansion isorder N L i N χ N a NLO 1 0 0N LO 7 2 2N LO 0 3 3Total 8 5 5 . N L i is the number of Gasser-Leutwyler coefficients, N χ the number of counterterms associated with the contin-uum χ PT expansion and N a is the number of coun-terterms associated with the discretization corrections.In total, there are 18 unknown LECs. While we uti-lize 18 ensembles in this analysis, the span of parameterspace is not sufficient to constrain all the LECs with-out prior knowledge. In particular, the introduction ofall 8 L i coefficients requires prior widths informed fromphenomenology.In the literature, the L i are typically quoted at therenormalization scale µ ρ = 770 MeV while in our work,we use the scale µ F = 4 πF . We can use the BE14values of the L i LECs from Ref. [77] and the known scaledependence [68] to convert them from µ ρ to µ F : L ri ( µ ) = L ri ( µ ) − Γ i (4 π ) ln (cid:18) µ µ (cid:19) , (4.1)1 TABLE V. The Γ i coefficients that appear in the scale dependence of the L i ( µ ). We evolve the L i ( µ ) from the typicalscale µ = 770 MeV, Eq. (3.21) to µ = 4 πF , beginning with the BE14 estimates from the review [77] (Table 3), using theirknown scale dependence [68], Eq. (4.1). We assign the following slightly more conservative uncertainty as a prior width in theminimization: If a value of L i is less than 0 . × − , we assign it a 100% uncertainty at the scale µ = 770 MeV; If the valueis larger than 0 . × − , we assign it the larger of 0.5 or 1/3 of the mean value. L i L L L L L L L L Γ i L i ( m ρ ) 0.53(50) 0.81(50) -3.1(1.0) 0.30(30) 1.01(50) 0.14(14) -0.34(34) 0.47(47)10 L i ( µ ) 0.37(50) 0.49(50) -3.1(1.0) 0.09(30) 0.38(50) 0.01(14) -0.34(34) 0.29(47) with the values of Γ i listed in Table V for convenience.We use F = 80 MeV, which is the value adopted byFLAG [8]. We set the central value of all the L i withthis procedure and the widths are set as described inTab. V.Next, we must determine priors for the N LO andN LO local counterterm coefficients, A nK,π,s . We set thecentral value of all these priors to 0 and then perform asimple grid search varying the widths to find preferredvalues of the width, as measured by the Bayes Factor.Our goal is not to optimize the width of each prior indi-vidually for each model used in the fit, but rather find aset of prior widths that is close to optimal for all models.To this end, we vary the width of the χ PT LECs to-gether at each order (N LO, N LO) and the discretiza-tion LECs together at each order (N LO, N LO) for afour-parameter search. We apply a very crude grid wherethe values of the widths are taken to be 2, 5, or 10.We find taking the width of all these A nK,π,s LECs equalto 2 results in good fits with near-optimal values. Thisprovides evidence the normalization of small parameterswe have chosen for (cid:15) P and (cid:15) a , Eq. (3.1), is “natural” andsupports the power counting we have assumed, Eq. (3.2).The N LO LECs mostly favor a width of 2 while theN LO discretization LECs prefer 5 and the N LO χ PTLECs vary from model to model with 5 a reasonable valuefor all. As a result of this search, we pick as our priors˜ A K,π = 0 ± , ˜ A s = 0 ± , ˜ A K,π = 0 ± , ˜ A s = 0 ± . (4.2) B. Physical Point
As our calculation is performed with isospin symmet-ric configurations and valence quarks, we must define aphysical point to quote our final result. We adopt the def-inition of the physical point from FLAG. FLAG[2017] [78]defines the isospin symmetric pion and kaon masses to be(Eq. (16)) ¯ M π = 134 . , ¯ M K = 494 . . (4.3)The values of F π + and F K + are taken from the N f = 2+1results from FLAG[2020] [8] (we divide the values by √ F phys π + = 92 . ,F phys K + = 110 . . (4.4)The isospin symmetric physical point is then defined byextrapolating our results to the values (for the choice F → F π ) ( (cid:15) phys π ) = (cid:32) ¯ M π πF phys π + (cid:33) , ( (cid:15) phys K ) = (cid:32) ¯ M K πF phys π + (cid:33) . (4.5) C. Model Averaging Procedure
Our model average is performed under a Bayesianframework following the procedure described in [39, 79].Suppose we are interested in estimating the posterior dis-tribution of Y = F K /F π , ie. P ( Y | D ) given our data D .To that end, we must marginalize over the different mod-els M k . P ( Y | D ) = (cid:88) k P ( Y | M k , D ) P ( M k | D ) (4.6)Here P ( Y | M k , D ) is the distribution of Y for a givenmodel M k and dataset D , while P ( M k | D ) is the posteriordistribution of M k given D . The latter can be written,per Bayes’ theorem, as P ( M k | D ) = P ( D | M k ) P ( M k ) (cid:80) l P ( D | M l ) P ( M l ) . (4.7)We can be more explicit with what the latter is in thecontext of our fits. First, mind that we are a priori ag-nostic in our choice of M k . We thus take the distribution P ( M k ) to be uniform over the different models. We cal-culate P ( D | M k ) by marginalizing over the parameters(LECs) in our fits: P ( D | M k ) = (cid:90) (cid:89) j d θ ( k ) j P ( D | θ ( k ) j , M k ) P ( θ ( k ) j | M k ) . (4.8)2After marginalization, P ( D | M k ) is just a number.Specifically, it is the Bayes Factor of M k : P ( D | M k ) =exp( logGBF ) M k , where logGBF is the log of the BayesFactor as reported by lsqfit [80]. Thus P ( M k | D ) = exp( logGBF ) M k (cid:80) Kl =1 exp( logGBF ) M l (4.9)with K the number of models included in our average.We emphasize that this model selection criterium notonly rates the quality of the decription of data, it alsopenalizes parameters which do not improve this descrip-tion. This helps rule out models which overparameterizedata.Now we can estimate the expectation value and vari-ance of Y .E[ Y ] = (cid:88) k E[ Y | M k ] P ( M k | D ) (4.10)Var[ Y ] = (cid:34)(cid:88) k Var[ Y | M k ] P ( M k | D ) (cid:35) (4.11)+ (cid:34)(cid:32)(cid:88) k E [ Y | M k ] P ( M k | D ) (cid:33) − E [ Y | D ] (cid:35) The variance V [ Y ] results from the total law of variance;the first term in brackets is known as the expected valueof the process variance (which we refer to as the modelaveraged variance ), while the latter is the variance of thehypothetical means (the root of which we refer to as the model uncertainty ). D. Full analysis and uncertainty breakdown
In total, we consider 216 different models of extrapo-lation/interpolation to the physical point. The variouschoices for building a χ PT or MA EFT model consist of × χ PT or MA EFT at NLO × F = { F π , F π F K , F K } in defining (cid:15) P × × LO, use full χ PT or just counterterms × α S term at N LO × × LO counterterms or not192 : total choices . We also consider pure Taylor expansion fits with onlycounterterms and no log corrections. For these fits, theset of models we explore is × LO or N LO × F = { F π , F π F K , F K } in defining (cid:15) P × α S term at N LO × . Based upon the quality of fit (gauged by the Bayesiananalogue to the p -value, Q , or the reduced chi square, χ ν ) and/or the weight determined as discussed in theprevious section, we can dramatically reduce the numberof models used in the final averaging procedure. First,any model which does not include the FV correction fromNLO is heavily penalized. This is not surprising given theobserved volume dependence on the a12m220 ensembles,Fig. 1. However, even if we remove the a12m220S en-semble from the analysis, the Taylor-expanded fits havea relative weight of e − or less compared to those thathave χ PT form at NLO.If we add FV corrections to the Taylor expansion fits(pure counterterm) and use all ensembles, F K F π = 1 + ¯ L ( (cid:15) K − (cid:15) π ) (cid:26) t FV (cid:88) | n |(cid:54) =0 c n m π L | n | K ( m π L | n | ) (cid:27) + · · · (4.12)they still have weights which are ∼ e − over the normal-ized model distribution and also contribute negligibly tothe model average.We observe that the fits which use the MA EFT atNLO are also penalized with a relative weight of ∼ e − ,and fits which only work to N LO have unfavorableweights by ∼ e − (and are also accompanied by poor χ ν values). Cutting all of these variations reduces ourfinal set of models to be N LO χ PT with the followingvariations × F = { F π , F π F K , F K } in defining (cid:15) P × × LO, use full χ PT or just counterterms × α S term at N LO24 : total choices which enter the model average . The final list of models, with their corresponding weightsand resulting extrapolated values to the isospin symmet-ric physical point, is given in Tab. VII in Appendix A.Our final result in the isospin symmetric limit, definedas in Eq. (4.5) and analogously for other choices of F ,including a breakdown in terms of statistical ( s ), pionmass extrapolation ( χ ), continuum limit ( a ), infinite vol-ume limit ( V ), physical point (phys) and model selection( M ) uncertainties, is as reported in Eq. (1.3) F K F π = 1 . s (12) χ (20) a (01) V (15) phys (12) M = 1 . . The finite volume uncertainty is assessed by removingthe a12m220S ensemble from the analysis, repeating themodel averaging procedure and taking the difference.The final probability distribution broken down into thethree choices of F is shown in Fig. 2.
1. Impact of a06m310L ensemble
Next, we turn to understanding the impact of thea06m310L ensemble on our analysis. The biggest dif-3 .
17 1 .
18 1 .
19 1 .
20 1 .
21 1 . F K /F π B a y e s M o d e l A v g P D F F → F π F → F π F K F → F K FIG. 2. Final probability distribution giving rise to Eq. (1.3),separated into the three choices of F = { F π , F π F K , F K } inthe definition of the small parameters, Eq. (3.1). The parent“gray” distribution is the final PDF normalized to 1 whenintegrated. ference upon removing the a06m310L ensemble is thatthe data are not able to constrain the various terms con-tributing to the continuum extrapolation as well, partic-ularly since there are up to three different types of scalingviolations: ( (cid:15) K − (cid:15) π ) × { (cid:15) a , α S (cid:15) a , (cid:15) a } , and thus, the statistical uncertainty of the results growsas well as the model variance, with a total uncertaintygrowth from ∼ . ∼ . LO fits become ac-ceptable, though they are still grossly outweighed by theN LO fits. Including both effects, the final model averageresult shifts from F K F π = 1 . → F K F π (cid:12)(cid:12)(cid:12)(cid:12) no a06 = 1 . . (4.13)In Fig. 3, we show the continuum extrapolation fromthree fits: • (Left): all ensembles, N LO χ PT with only countert-erms at N LO and N LO and F = F π ; • (Middle): no a06m310L, N LO χ PT with only coun-terterms at N LO and N LO and F = F π ; • (Right): no a06m310L, N LO χ PT with only countert-erms at N LO and F = F π .As can be seen from the middle plot, the a15, a12 anda09 ensembles prefer contributions from both (cid:15) a and (cid:15) a contributions and are perfectly consistent with the resulton the a06m310L ensemble. They are also consistent withan N LO fit (no (cid:15) a contributions) as can be seen in theright figure. However, the weight of the N LO fits is stillsignificantly greater than the N LO fits even without thea06m310L data. We conclude that the a06m310L ensemble is useful,but not necessary to obtain a sub-percent determinationof F K /F π with our lattice action. A more exhaustivecomparison can be performed with the analysis notebookprovided with this publication.In Fig. 4, we show the stability of our final result forvarious choices discussed in this section.
2. Convergence of the chiral expansion
While the numerical analysis favors a fit function inwhich only counterterms are used at N LO and higher,it is interesting to study the convergence of the chiralexpansion by studying the fits which use the full χ PTexpression at N LO.In Fig. 5, we show the resulting light quark mass de-pendence using the N LO extrapolation with the fullN LO χ PT formula. After the analysis is performed,the results from each ensemble are shifted to the phys-ical kaon mass point, leaving only dependence upon (cid:15) π and (cid:15) a as well as dependence upon the η mass defined bythe GMO relation. The magenta band represents the full68% confidence interval in the continuum, infinite volumelimit. The different colored curves are the mean values asa function of (cid:15) π at the four different lattice spacings. Wealso show the convergence of this fit in the lower panelplot. From this convergence plot, one sees that roughlythat at the physical pion mass (vertical gray line) theNLO contributions add a correction of ∼ .
16 comparedto 1 at LO, the N LO contributions add another ∼ . LO corrections are not detectible by eye. Theband at each order represents the sum of all terms upto that order determined from the full fit. The reduc-tion in the uncertainty as the order is increased is due onlarge part to the induced correlation between the LECsat different orders through the fitting procedure.In Fig. 2, we observe that the different choices of F areall consistent, indicating higher-order corrections (start-ing at N LO in the non-counterterm contributions) aresmaller than the uncertainty in our results. It is also in-teresting to note that choosing F πK or F K is penalizedby the analysis, indicating the numerical results preferlarger expansion parameters. In Tab. VI, we show the re-sulting χ PT LECs determined in this analysis for the twochoices F = { F π , F Kπ } , as well as whether the ratio formof the fit is used, Eq. (3.17). For the Gasser-LeutwylerLECs, we evolve the values back from µ → µ ρ for a sim-pler comparison with the values quoted in literature. Formost of the L i , we observe the numerical results have verylittle influence on the parameters as they mostly returnthe prior value (also listed in the table for convenience).The only LECs influenced by the fit are L , L , and L with L getting pulled about one sigma away from theprior value and L and L only shifting by a third orhalf of the prior width. What is most interesting fromthese results is that our fit prefers a value of L that isnoticeably smaller than the value obtained by MILC [16]4 .
000 0 .
025 0 .
050 0 .
075 0 .
100 0 .
125 0 .
150 0 .
175 0 . (cid:15) a = a / (2 w ) . . . . . . . . F K / F π xpt nnnlo FV ct PP a06 a09 a12 a15 .
000 0 .
025 0 .
050 0 .
075 0 .
100 0 .
125 0 .
150 0 .
175 0 . (cid:15) a = a / (2 w ) . . . . . . . . F K / F π xpt nnnlo FV PP a06 a09 a12 a15 .
000 0 .
025 0 .
050 0 .
075 0 .
100 0 .
125 0 .
150 0 .
175 0 . (cid:15) a = a / (2 w ) . . . . . . . . F K / F π xpt nnlo FV PP a06 a09 a12 a15 FIG. 3. LEFT: An N LO fit to all ensembles. MIDDLE: The same fit to all ensembles excluding a06m310L. RIGHT: Arepresentative N LO fit to all ensembles excluding a06m310L. In all plots, the results from each ensemble are shifted to thephysical values of (cid:15) π and (cid:15) K and the infinite volume limit, with only the (cid:15) a dependence remaining. The labels are are explainedin Appendix A; the data points at each spacing are slightly offset horizontally for visual clarity.TABLE VI. Resulting LECs from full N LO χ PT analy-sis (also including N LO counterterms). For the Gasser-Leutwyler LECs L i , we evolve them back to the standardscale µ = 770 MeV, while for the other LECs, we leave themat the scale µ = 4 πF (cid:39) LEC F = F π F = F π F K ratio ratio µ = 770 prior no yes no yes10 L . . . . . L . . . . . L − . . − . − . − . − . L . . . . . L . . . . . L . . . . . L − . − . − . − . − . L . . . . . µ = µ A K . .
42) 0 . .
41) 0 . .
6) 0 . . A π . .
2) 2 . .
2) 2 . .
3) 2 . . A Kπ . .
7) 2 . .
7) 1 . .
7) 2 . . A K . . . .
0) 0 . .
4) 0 . . A p . .
0) 2 . .
1) 2 . .
4) 2 . . and HPQCD [12] and also smaller than the BE14 resultfrom Ref. [77].In Fig. 6, we show the impact of using the fully ex-panded expression, Eq. (3.15) versus the expression inwhich the NLO terms are kept in a ratio, Eq. (3.17).To simplify the comparison we restrict it to the choice F = F π and the full N LO χ PT expression. We see thatfits without the ratio form are preferred, but the cen-tral value of the final result depends minimally upon thischoice.In Fig. 7, we show that the results strongly favor theuse of only counterterms at N LO as opposed to the full χ PT fit function at that order. We focus on the choice F = F π to simplify the comparison.Our results are not sufficient to understand why the fitfavors only counterterms at N LO and higher. While thelinear combination of LECs in Eq. (3.12) are redundant,the L i LECs also appear in the single-log coefficients,Eqs. (3.10) and (3.11) in different linear combinations.Neverthelss, we double check that the fit is not penalizedfor the counterterm redundancy, Eq. (3.12). Using the priors for L i from Tab. V, we find the contribution fromthe Gasser-Leutwyler LECs to these N LO counterterms,Eq. (3.14) are given by L K = 0 . . , L π = − . . (4.14)As the A P terms are priored at 0(2), it is sufficient tore-run the analysis by simply setting L P = 0. We findthis result marginally improves the Bayes Factors but notstatistically significantly, leaving us with the puzzle thatthe optimal fit is a hybrid NLO χ PT plus countertermsat higher orders.If the Taylor expansion fits (pure counterterm) weregood and favored over the χ PT fits, this could be asign that the SU (3) χ PT formula was failing to de-scribe the lattice results. However, we have to includethe NLO χ PT expression, including its predicted (coun-terterm free) volume dependence to describe the numer-ical results. It would be nice to have the full N LO MAEFT expression to understand why the hybrid MA EFTfits are so relatively disfavored in the analysis. Theremay be compensating discretization effects that cancelagainst those at NLO to some degree that might allowthe full N LO MA EFT to better describe the results. Itis unlikely, however, that this expression will be derivedas at two loops in χ PT, the universality of MA EFTexpressions [58] breaks down and they can no longer be“derived” from the corresponding partially quenched for-mula, which is known for F K /F π at two loops [76, 81–83]. E. QCD isospin breaking corrections
Finally, we discuss the correction to our result to ob-tain a direct determination of F K + /F π + including strongisospin breaking corrections, but excluding QED correc-tions. This is the standard value quoted in the FLAGreviews [8, 78]. Our calculation, like most, are performedin the isospin symmetric limit, and therefore, the strongisospin breaking correction must be estimated, ratherthan having a direct determination. The optimal ap-proach is to incorporate both QED and QCD isospinbreaking corrections into the calculations such that the5 .
19 1 . F K /F π model average ct N LO χ PTfull N LO χ PT F = F π F = F π F K F = F K w α S w/o α S Taylor expandedNLO ratio a & . fmN LON LOMA EFTwidths: N LO N LO .
00 0 .
25 0 .
50 0 .
75 1 . Weight Relative to Highest
FIG. 4. Stability plot of final result compared to various model choices. The black square at the top is our final answer forthe isospin symmetric determination of F K /F π . The vertical magenta band is the uncertainty of this fit to guide the eye. Thesolid magenta squares are various ways of decomposing the model selection that goes into the final average. The right panelshows the relative weight with respect to the maximum logGBF value, exp { logGBF i − logGBF max } . Below the set of modelsincluded in the average, we show sets of analyses that are not included in the average for comparison, which are indicatedwith open gray symbols. First, we show the impact of excluding the a06m310L ensemble. The logGBF can not be directlycompared between these fits and the main analysis as the number of data points are not the same, so the overall normalizationis different; their logGBF are shown as open diamonds. The relative logGBF between the N LO and N LO analysis can becompared which indicates a large preference for the N LO analysis. We also show the MA EFT analysis, which agrees wellwith the main analysis. Finally, we show the results if one were to change the widths of the N LO and N LO priors from thosechosen in Eq. (4.2). separation is not necessary, as was done in Ref [84] byincorporating both types of corrections through the per-turbative modification of the path integral and correla-tion functions [85, 86]. In this work, we have not per-formed these extensive computations and so we rely uponthe SU (3) χ PT prediction to estimate the correctiondue to strong isospin breaking. As we have observedin Sec. IV D, the SU (3) chiral expansion behaves andconverges nicely, so we expect this approximation to bereasonable.The NLO corrections to F K and F π including the strong isospin breaking corrections are given by δF NLO π ± = − (cid:96) π − (cid:96) π ± − (cid:96) K − (cid:96) K ±
4+ 4 ¯ L ( (cid:15) π ± + (cid:15) K ± + (cid:15) K ) + 4 ¯ L (cid:15) π ± ,δF NLO K ± = − (cid:96) π − (cid:96) π ± − (cid:96) η − (cid:96) K − (cid:96) K ±
2+ 14 ( (cid:15) K − (cid:15) K ± ) (cid:96) η − (cid:96) π (cid:15) η − (cid:15) π + 4 ¯ L ( (cid:15) π ± + (cid:15) K ± + (cid:15) K ) + 4 ¯ L (cid:15) K ± , (4.15)where we have kept explicit the contribution from eachflavor of meson propagating in the loop. There are threepoints to note in these expressions:1. At NLO in the SU (3) chiral expansion, there areno additional LECs that describe the isospin break-ing corrections beyond those that contribute to theisospin symmetric limit. Therefore, one can make6 .
00 0 .
02 0 .
04 0 .
06 0 . (cid:15) π = ( m π / πF π ) . . . . . . . . F K / F π xpt nnnlo FV PP a06a09a12a15 .
00 0 .
02 0 .
04 0 .
06 0 . (cid:15) π = ( m π / πF π ) . . . . . . . . F K / F π xpt nnnlo FV PP NLON LO N LO FIG. 5. Sample light quark mass dependence from a χ PTfit with F = F π (N LO χ PT + N LO counterterms). TOP:The curves are plotted with (cid:15) K = ( (cid:15) phys K ) , at fixed (cid:15) a for eachlattice spacing, as a function of (cid:15) π . The magenta band is thefull uncertainty in the continuum, infinite volume limit. Thedata points have all been shifted from the values of (cid:15) latt .K to (cid:15) phys K and to the infinite volume limit. BOTTOM: We showthe convergence of the resulting fit as a function of (cid:15) π . Eachband corresponds to all contributions up to that order withthe LECs determined from the full fit. The N LO band corre-sponds to the continuum extrapolated band in the top figure. .
17 1 .
18 1 .
19 1 .
20 1 .
21 1 . F K /F π P D F s f r o m r a t i oo r n o t w/o ratio fit ( F π ) w/ ratio fit ( F π ) FIG. 6. Comparison of fits with the fully expanded Eq. (3.15)and ratio Eq. (3.17) expressions, all with the choice F = F π .The PDFs are taken from the parent PDF, Fig. 2 withoutrenormalizing such that height in this figure reflects the rela-tive weight comopared to the total PDF. .
17 1 .
18 1 .
19 1 .
20 1 .
21 1 . F K /F π P D F s f r o m N L O c t o r f u ll X P T N2LO=ct ( F π ) N2LO=XPT ( F π ) FIG. 7. Comparison of N LO χ PT analysis with F = F π using the full N LO χ PT expression (smaller histogram) ver-sus only counterterms at N LO, Eq. (3.12) and (3.31). As inFig. 6, the PDFs are drawn from the parent PDF. a parameter-free prediction of the isospin breakingcorrections using lattice results from isospin symmet-ric calculations, with the only assumption being that SU (3) χ PT converges for this observable;2. If we expand these corrections about the isospin limit,they agree with the known results [5], and δF π ± is freeof isospin breaking corrections at this order;3. We have used the kaon mass splitting in place of thequark mass splitting, which is exact at LO in χ PT B ( m d − m u ) = ( ˆ M K − ˆ M K ± );The estimated shift of our isospin-symmetric result toincorporate strong isospin breaking is then δF iso K − π ≡ F ˆ K + F ˆ π + − F K F π = −
14 ( (cid:96) ˆ K + − (cid:96) ¯ K ) + 4 ¯ L ( (cid:15) K + − (cid:15) K )+ 14 ( (cid:15) K − (cid:15) K ± ) (cid:96) η − (cid:96) π (cid:15) η − (cid:15) π (4.16)Ref. [5] suggested replacing ¯ L with the NLO expressionequating it to the isospin symmetric F K /F π which yields δF iso (cid:48) K − π = − (cid:15) K − (cid:15) K ± (cid:15) η − (cid:15) π (cid:20) (cid:18) F K F π − (cid:19) + (cid:15) π ln (cid:18) (cid:15) K (cid:15) π (cid:19) − (cid:15) K + (cid:15) π (cid:21) . (4.17)In this expression, we have utilized the two relations (cid:96) ˆ K + − (cid:96) ¯ K = − (cid:15) K − (cid:15) K ± (cid:15) η − (cid:15) π ( (cid:15) K − (cid:15) π )(ln (cid:15) K + 1) (cid:15) η − (cid:15) π = 43 ( (cid:15) K − (cid:15) π ) (4.18)At this order, both of these expressions, Eqs. (4.16) and(4.17) are equivalent. However, they can result in shifts7that differ by more than one standard deviation. Fur-ther, the direct estimate of the strong isospin breakingcorrections [9] is larger in magnitude than either of them.Therefore, to estimate the strong isospin breaking cor-rections, we take the larger of the two corrections as themean and the larger uncertainty of the two, and thenadd an additional 25% uncertainty for SU (3) truncationerrors. In Sec. IV D 2 we observe the N LO correction is ∼
25% of the NLO correction (while NLO is ∼
16% ofLO).In order to evaluate these expressions, we have to de-fine the physical point with strong isospin breaking andwithout QED isospin breaking. We employ the valuesfrom FLAG[2017] [78] (except ˆ M π = 134 . M π = ˆ M π + = 134 . , ˆ M K = 497 . , ˆ M K + = 491 . . (4.19)With this definition of the physical point, we find (underthe same model average as Tab. VII) δF iso K − π = − . ,δF iso (cid:48) K − π = − . , (4.20)resulting in our estimated strong isospin breaking correc-tion F ˆ K + F ˆ π + − F K F π = − . F ˆ K + F ˆ π + = 1 . iso = 1 . , where the first uncertainty in the first line is the combi-nation of those in Eq. (1.3). V. SUMMARY AND DISCUSSION
The ratio F K /F π may be used, in combination withexperimental input for leptonic decay widths, to makea prediction for the ratio of CKM matrix elements, | V us | / | V ud | . Using the most recent data, Eq. (1.1) be-comes [6], | V us || V ud | F ˆ K + F ˆ π + = 0 . , (5.1)where strong isospin breaking effects must be included fordirect comparison with experimental data. Combiningthis expression with our final result, we find, | V us || V ud | = 0 . . (5.2) .
955 0 .
960 0 .
965 0 .
970 0 .
975 0 . | V ud | . . . . . . | V u s | FIG. 8. Result for the ratio of CKM matrix elements, | V us | / | V ud | , extracted from the ratio F K /F π reported in thiswork (red band). The global lattice value for | V us | extractedfrom a semileptonic decay form factor, f + (0) [8], is shown asa horizontal blue band, while the global experimental aver-age for | V ud | from nuclear beta decay [6] is given as a verticalgreen band. Note that the intersection between the red andgreen bands agrees well with the unitarity constraint for theCKM matrix, while the intersection between the red and bluebands shows ∼ σ tension. Utilizing the current global average, | V ud | =0 . | V us | = 0 . . (5.3)Finally, we may use our results, combined with the value | V ub | = (3 . × − , as a test of unitarity for theCKM matrix, which states that | V ud | + | V us | + | V ub | =1. From our calculation we find, | V ud | + | V us | + | V ub | = 0 . . (5.4)Alternatively, rather than using the experimental de-termination of | V ud | as input for our test of unitar-ity, we may instead use the global lattice average for | V us | = 0 . f + (0),the zero momentum transfer limit of a form factor rele-vant for the semileptonic decay K → π − lν . This leadsto, | V ud | + | V us | + | V ub | = 0 . , (5.5)leading to a roughly 2 σ tension with unitarity. Our re-sult, along with with the reported experimental resultsfor | V ud | and lattice results for | V us | , are shown in Fig. 8.One could also combine our results with the more preciseaverage in the FLAG review which would lead to a slightreduction their reported uncertainties, but we will leavethat to the FLAG Collaboration in their next update.Another motivation for this work was to precisely test(below 1%) whether the action we have used for our nu-cleon structure calculations [38–40] can be used to repro-duce an accepted value from other lattice calculationsthat are known at the sub-percent level. Our result pro-vides the first sub-percent cross-check of the universality8 .
00 0 .
05 0 .
10 0 .
15 0 . (cid:15) a = a / (2 w ) . . . . . . . . F K + / F π + FLAG[2020]: N F = 2 + 1 + 1 a06 a09 a12 a15 FIG. 9. A comparision of the continuum extrapolation from the MILC[2014] [87] result (LEFT) with the continuum extrap-olation in the present work with the MDWF on gradient-flowed HISQ action (right). In the right plot, we also include theFLAG[2020] average value from the N f = 2 + 1 + 1 calculations [8]. While the a06m310L ensemble is not necessary for us toextrapolate to a consistent value as this FLAG average (see Fig. 3), the overall size of our discretization effects are larger. Thisis not necessarily surprising as the HISQ action used by MILC has perturbatively removed all O( a ) corrections such that theleading scaling violations begin at O( α S a ), as implied by the x -axis of the left plot. of the continuum limit of this quantity with N f = 2+1+1dynamical flavors, albeit with the same sea-quark actionas used by MILC/FNAL and HPQCD [12, 13].Critical in obtaining a sub-percent determination ofany quantity is control over the continuum extrapolation.This is relevant to our pursuit of a sub-percent determi-nation of g A as another calculation, utilizing many of thesame HISQ ensembles but with a different valence ac-tion (clover fermions), obtains a result that is in tensionwith our own [88, 89]. While there has been speculationthat this discrepancy is due to the continuum extrapo-lations [89], new work suggests the original work under-estimated the systematic uncertainty in the correlationfunction analysis, and when accounted for, the tensionbetween our results goes away [90].In either case, to obtain a sub-percent determination of g A , which is relevant for trying to shed light on the neu-tron lifetime discrepancy [91], it is important to under-stand the scaling violations of our lattice action. While asmooth continuum extrapolation in one observable doesnot guarantee such a smooth extrapolation in another, itat least provides some reassurance of a well-behaved con-tinuum extrapolation. Furthermore, the determination of F K /F π involves the same axial current that is relevantfor the computation of the nucleon matrix element usedto compute g A .Fig. 3 shows the continuum extrapolation of F K /F π from our analysis. The size of the discretization effectsare noticeably larger than we observed in our calcula-tion of g A [39]. In Sec. IV D 1, we demonstrated that,while helpful, the a06m310L ensemble is not necessaryto achieve a sub-percent determination of F K /F π . Thisis in contrast to the determination by MILC which re-quires the a ∼ .
06 fm (or smaller) lattice spacings tocontrol the continuum extrapolation (though we note, theHPQCD calculation [12], also performed on the HISQ en-sembles, does not utilize the a ∼ .
06 fm ensembles butagrees with the MILC result). It should be noted, theMILC result does not rely on the heavier mass ensemblesexcept to adjust for the slight mistuning of the inputquark masses on their near-physical point ensembles. InFig. 9, we compare our continuum extrapolation to thatof MILC [87].In Ref. [87], they also utilize the same four lattice spac-ings as in this work (they have subsequently improvedtheir determination with an additional two finer lattice9spacings [13].) A strong competition between the O( a )and O( a ) corrections was observed in that work, suchthat the a ∼ .
06 fm ensemble is much more instrumentalfor a reliable continuum extrapolation than is the case inour setup. At the same time, the overall scale of their dis-cretization effects is much smaller than we observe in theMDWF on gradient-flowed HISQ action for this quan-tity. This is not entirely surprising as the HISQ actionhas been tuned to perturbatively remove all O( a ) cor-rections such that the leading corrections formally beginas O( α S a ). ACKNOWLEDGMENTS
We would like to thank V. Cirigliano, S. Simula, J.Simone, and T. Kaneko for helpful correspondence anddiscussion regarding the strong isospin breaking correc-tions to F K /F π . We would like to thank J. Bijnens forhelpful correspondence on χ PT and a
C++ interface to
CHIRON [69] that we used for the analysis presented inthis work. We thank the MILC Collaboration for provid-ing some of the HISQ configurations used in this work,and A. Bazavov, C. Detar and D. Toussaint for guidanceon using their code to generate the new HISQ ensemblesalso used in this work. We would like to thank P. Lep-age for enhancements to gvar [92] and lsqfit [80] thatenable the pickling of lsqfit.nonlinear fit objects.Computing time for this work was provided throughthe Innovative and Novel Computational Impact on The-ory and Experiment (INCITE) program and the LLNLMultiprogrammatic and Institutional Computing pro-gram for Grand Challenge allocations on the LLNL su- percomputers. This research utilized the NVIDIA GPU-accelerated Titan and Summit supercomputers at OakRidge Leadership Computing Facility at the Oak RidgeNational Laboratory, which is supported by the Office ofScience of the U.S. Department of Energy under Con-tract No. DE-AC05-00OR22725 as well as the Surface,RZHasGPU, Pascal, Lassen, and Sierra supercomputersat Lawrence Livermore National Laboratory.The computations were performed utilizing
LALIBE [93] which utilizes the
Chroma software suite [94]with
QUDA solvers [24, 25] and HDF5 [95] for I/O [96].They were efficiently managed with
METAQ [97, 98] andstatus of tasks logged with EspressoDB [99]. The HMCwas performed with the MILC Code [100], and for the en-sembles new in this work, running on GPUs using
QUDA .The final extrapolation analysis utilized gvar v11.2 [92]and lsqfit v11.5.1 [80] and
CHIRON v0.54 [69]. The anal-ysis and data for this work can be found at this git repo: https://github.com/callat-qcd/project_fkfpi .This work was supported by the NVIDIA Corpo-ration (MAC), the Alexander von Humboldt Founda-tion through a Feodor Lynen Research Fellowship (CK),the DFG and the NSFC Sino-German CRC110 (EB),the RIKEN Special Postdoctoral Researcher Program(ER), the U.S. Department of Energy, Office of Sci-ence, Office of Nuclear Physics under Award NumbersDE-AC02-05CH11231 (CCC, CK, BH, AWL), DE-AC52-07NA27344 (DAB, DH, ASG, PV), DE-FG02-93ER-40762 (EB), DE-AC05-06OR23177 (BJ, CM, KO), DE-FG02-04ER41302 (KO); the Office of Advanced ScientificComputing (BJ); the Nuclear Physics Double Beta De-cay Topical Collaboration (DAB, HMC, AN, AWL); andthe DOE Early Career Award Program (CCC, AWL).
Appendix A: Models included in final analysis
We list the models that have entered the final analysis as described in Sec. IV D and listed in Tab. VII. For example,the model xpt-ratio nnnlo FV alphaS PPindicates the model uses the continuum χ PT fit function through N LO with discretization corrections added as inEqs. (3.31) and (3.40). The NLO contributions are kept in a ratio form, Eq. (3.17) and we have included the corre-sponding N LO ratio correction δ N LOratio . The finite volume corrections have been included at NLO. The discretizationterms at N LO include the α S (cid:15) a ( (cid:15) K − (cid:15) π ) counterterm. The renormalizations scale appearing in the logs is µ = 4 πF π ,indicated by PP and we have included the corresponding N LO correction δ N LO F π , Eq. (3.29), to hold the actualrenormalization scale fixed at µ = 4 πF .If ct appears in the model name, it means that the only N LO terms that are added are from the local countertermsand all chiral log corrections are set to zero.
Appendix B: NLO Mixed Action Formulas
The expression for d(cid:96) π arises from the integral d(cid:96) π = (cid:90) R d d k (2 π ) d i ( k − m π ) = 1 + ln( m π /µ )(4 π ) , (B1)0 TABLE VII. List of models used in final result, as described in the text.model χ ν Q logGBF weight F K /F π xpt-ratio nnnlo FV ct PP 0.847 0.645 77.728 0.273 1.1968(40)xpt-ratio nnnlo FV alphaS ct PP 0.843 0.650 77.551 0.229 1.1962(46)xpt nnnlo FV ct PP 0.908 0.569 76.830 0.111 1.1974(40)xpt nnnlo FV alphaS ct PP 0.902 0.576 76.668 0.095 1.1966(46)xpt-ratio nnnlo FV ct PK 1.014 0.439 76.343 0.068 1.1952(37)xpt-ratio nnnlo FV alphaS ct PK 1.006 0.449 76.234 0.061 1.1944(42)xpt nnnlo FV PP 0.949 0.517 75.371 0.026 1.1989(40)xpt nnnlo FV alphaS PP 0.946 0.522 75.196 0.022 1.1983(46)xpt nnnlo FV ct PK 1.135 0.309 75.084 0.019 1.1950(36)xpt nnnlo FV alphaS ct PK 1.123 0.321 75.007 0.018 1.1941(41)xpt-ratio nnnlo FV PP 1.014 0.439 74.765 0.014 1.1987(40)xpt-ratio nnnlo FV alphaS PP 1.009 0.445 74.599 0.012 1.1980(46)xpt nnnlo FV PK 1.100 0.344 74.421 0.010 1.1969(37)xpt nnnlo FV alphaS PK 1.093 0.352 74.306 0.009 1.1962(42)xpt-ratio nnnlo FV ct KK 1.262 0.202 74.014 0.007 1.1920(36)xpt-ratio nnnlo FV alphaS ct KK 1.244 0.215 74.004 0.007 1.1912(39)xpt-ratio nnnlo FV PK 1.159 0.286 73.880 0.006 1.1967(37)xpt-ratio nnnlo FV alphaS PK 1.150 0.295 73.780 0.005 1.1959(41)xpt nnnlo FV KK 1.288 0.184 72.757 0.002 1.1938(36)xpt nnnlo FV alphaS KK 1.273 0.194 72.718 0.002 1.1930(40)xpt-ratio nnnlo FV KK 1.338 0.152 72.348 0.001 1.1938(36)xpt-ratio nnnlo FV alphaS KK 1.322 0.162 72.323 0.001 1.1929(39)xpt nnnlo FV alphaS ct KK 1.536 0.068 71.459 0.001 1.1900(38)xpt nnnlo FV ct KK 1.558 0.061 71.430 0.001 1.1909(35) Bayes Model Average which has been regulated and renormalized with the standard χ PT modified dimensional-regularization scheme [48].The finite volume corrections to δ(cid:96) π are given by δ FV d(cid:96) π = (cid:88) | n |(cid:54) =0 c n (4 π ) (cid:20) K ( mL | n | ) mL | n | − K ( mL | n | ) − K ( mL | n | ) (cid:21) (B2)The expression for K φ φ arises from the integral K φ φ = (4 π ) (cid:90) R d d k (2 π ) d i ( k − m φ )( k − m φ ) = (cid:96) φ − (cid:96) φ (cid:15) φ − (cid:15) φ . (B3)Similarly, K (2 , φ φ is given by K (2 , φ φ = (cid:90) R d d k (2 π ) d i (4 π ) (4 πF ) ( k − m φ ) ( k − m φ ) = (cid:96) φ − (cid:96) φ ( (cid:15) φ − (cid:15) φ ) − d(cid:96) φ (cid:15) φ − (cid:15) φ . (B4)Finally, K φ φ φ is given by K φ φ φ = (cid:90) R d d k (2 π ) d i (4 π ) (4 πF ) ( k − m φ )( k − m φ )( k − m φ )= (cid:96) φ ( (cid:15) φ − (cid:15) φ )( (cid:15) φ − (cid:15) φ ) + (cid:96) φ ( (cid:15) φ − (cid:15) φ )( (cid:15) φ − (cid:15) φ ) + (cid:96) φ ( (cid:15) φ − (cid:15) φ )( (cid:15) φ − (cid:15) φ ) . (B5)In each of these expressions, the corresponding expression including FV corrections are given by replacing (cid:96) φ → (cid:96) FV φ ,Eq. (3.38). Appendix C: HMC for new ensembles
We present various summary information for the three new ensembles used in this work, a06m310L, a15m135XL anda09m135. In Tab. VIII, we list the parameters of the HISQ ensembles used in the HMC. In Fig. 10, we show the MDTU1
TABLE VIII. Input parameters and measured acceptance rate for the new HISQ ensembles. In addition to the columnsstandardly reported by MILC (see Table IV of Ref. [43]), we list the abbreviated ensemble name, number of streams N stream and total number of configurations N cfg . For a given ensemble, each stream has an equal number of configurations. Thegauge coupling, light, strange, and charm quark masses on each ensembles are given as well as the tadpole factor, u and theNaik-term added to the charm quark action, (cid:15) N . The total length in molecular dynamics time units (MDTU) between eachsaved configuration is s while the length between accept/reject steps is Len. The microstep size (cid:15) used in the HMC is providedas Len./ N steps which was input with single precision. The average acceptance rate over all streams is listed as well as thenumber of streams.ensemble 10 /g am l am s am c u (cid:15) N s Len. (cid:15)
Acc. N stream N cfg a15m135XL 5.80 0.002426 0.06730 0.8447 0.85535 − − − MDTU − − ∆ S a15m135XL b acc. = 0.635c acc. = 0.631 d acc. = 0.625e acc. = 0.624 − − MDTU − . . . ∆ S a09m135 a acc. = 0.694 b acc. = 0.691 − . . .
50 1000 2000 3000 4000 5000 6000
MDTU − . . . ∆ S a06m310L b acc. = 0.775 c acc. = 0.754 − . . . FIG. 10. The ∆ S values computed in the accept/reject step of the HMC versus MDTU. The different colors correspond tothe different streams which are separated and shifted in MDTU for clarity. history of the ∆ S for the three ensembles. For the a15m135XL ensemble, we reduced the trajectory length significantlycompared to the a15m130 from MILC to overcome spikes in the HMC force calculations. To compensate, we loweredthe acceptance rate to encourage the HMC to move around parameter space with larger jumps, to hopefully reducethe autocorrelation time. We ran 25 HMC accept/reject steps before saving a configuration for a total trajectorylength of 5.For each accept/reject step we also measure the quark-anti-quark condensate ¯ ψψ using a stochastic estimate with5 random sources that are averaged together. We compute it for each of the quark masses am l , am s , and am c . On2 TABLE IX. Average values of the quark-anti-quark condensate ¯ ψψ with statistical errors and integrated autocorrelation times τ measured with the Γ-method analysis. Each value is averaged over all the available streams (which are all statisticallycompatible). The integrated autocorrelation time is reported in units of MDTU. The a15m135XL results are obtained fromthe second half of each stream because we have more measurements.ensemble N stream ¯ ψψ l τ l ¯ ψψ s τ s ¯ ψψ c τ c a15m135XL 4 0.02390(2) 11(2) 0.08928(2) 71(26) 0.4800580(5) < ψψ l ¯ ψψ s ¯ ψψ c HMC trajectory a15m135XL b c d e 0 100 0 2000 4000 6000 8000 10000
HMC trajectory a15m135XL b c d e 0 100 0 2000 4000 6000 8000 10000
HMC trajectory a15m135XL b c d e 0 1000 2000 4000 6000
HMC trajectory a09m135 a b 0 50 0 2000 4000 6000
HMC trajectory a09m135 a b 0 50 0 2000 4000 6000
HMC trajectory a09m135 a b 0 500 1000 2000 3000 4000 5000 6000
HMC trajectory a06m310L b c 0 50 0 1000 2000 3000 4000 5000 6000
HMC trajectory a06m310L b c 0 50 0 1000 2000 3000 4000 5000 6000
HMC trajectory a06m310L b c 0 50 ¯ ψψ l ¯ ψψ s ¯ ψψ c FIG. 11. The quark-anti-quark condensate on each configuration of the three ensembles. The different streams are plottedseparately for clarity. The plots in each column correspond to the light, strange and charm quark masses, respectively. the a15m135XL we have measured ¯ ψψ only on every saved configuration for the first half of each stream, while wemeasured it at each accept/reject step for the second half. The integrated autocorrelation time, as well as the averageand statistical errors of ¯ ψψ are computed using the Γ-method analysis [101] with the Python package unew [102]. Wereport the results in Tab. IX. In Fig. 11 we report the value of the ¯ ψψ on each saved configuration for the three quarkmasses on each ensemble.Because we observe a long autocorrelation time of the (cid:104) ¯ ψ s ψ s (cid:105) on the a15m135XL ensemble, we also studied theuncertainty on the extracted pion and kaon effective masses as a function of block size to check for possible longerautocorrelations than usual, with blocking lengths of 10, 25, and 100 MDTU. We observe that these hadronic quantitieshave a much shorter autocorrelation time as the uncertainty is independent of τ b and consistent with the unblockeddata. On this a15m135XL ensemble, while we have generated 2000 configurations, we have only utilized 1000 in thispaper (the first half from each of the four stream).Finally, in Tab. X, we list the parameters of the overrelaxed stout smearing used to measure the topological charge Q on each configuration [103] and we show the resulting Q distributions in Fig. 13. While the Q -distribution on thea06m310L ensemble is less than ideal and the intergrated autocorrelation time is long, the volume is sufficiently large3
20 22 24 26 28 t . . m π e ff ( t ) τ b = None τ b = τ b = τ b = . . . . m K e ff ( t ) a15m135XL τ b = None τ b = τ b = τ b = FIG. 12. The effective mass in the mid- to long-time region of the kaon (top) and pion (bottom) on the a15m135XL ensembleare plotted as a function of the blocking time, τ b in MDTU. For example, τ b = 10 blocks nearest neighbor configurations while τ b = 100 is blocking in groups of 20 configurations. That the uncertainty is independent of τ b indicates the autocorrelationtime for these hadronic quantities is very short.TABLE X. Values of the overrelaxed stout smearing parameters used to measure the topological charge Q and the resultingmean ( ¯ Q ) and width ( σ ) of the the distribution for each stream. The last column reports the integrated autocorrelation timein units of MDTU using the Γ-method analysis. These measurements were performed with QUDA which is now available via the su3 test test executable in the develop branch [24, 25].ensemble ρ N step ¯ Q stream ( σ ) ¯ Q all ( σ ) τ all ( σ )a b c d ea15m135XL 0 .
068 2000 – -10(34) -3(35) -2(33) 5(32) -3(33) 15(3)a09m135 0 .
065 2000 0.5(12.0) 2(12) – – – 1(12) 18(4)a06m310L 0 .
066 1800 – 4(12) -1.2(7.4) – – 1(10) 420(198) ( aL = 72 a (cid:39) . Q -values, whichnonetheless average to nearly 0. HMC trajectory − − Q a15m135XL b c d e − − HMC trajectory − − Q a09m135 a b − − HMC trajectory − − − Q a06m310L b c − − − FIG. 13. The distribution of topological charge Q measurements on each configuration. The Q -values were determined byusing the overrelaxed stout smearing technique outlined in [103] with weight parameters ρ given in Tab. X and ε = − . ρ parameter and the number of steps to perform on an ensemble-to-ensemble basis, i.e., for a handful of configurations per ensemble we choose a spread of ρ s and step numbers and observe whichcombination gives the best plateau. These values of ρ and step number ( ε is always -0.25) are then applied to the entireensemble.[1] William J. Marciano, “Precise determination of—V(us)— from lattice calculations of pseudoscalar de- cay constants,” Phys. Rev. Lett. , 231803 (2004), arXiv:hep-ph/0402299 [hep-ph].[2] C. Aubin, C. Bernard, Carleton E. DeTar, J. Osborn,Steven Gottlieb, E. B. Gregory, D. Toussaint, U. M.Heller, J. E. Hetrick, and R. Sugar (MILC), “Lightpseudoscalar decay constants, quark masses, and lowenergy constants from three-flavor lattice QCD,” Phys.Rev. D70 , 114501 (2004), arXiv:hep-lat/0407028 [hep-lat].[3] Roger Decker and Markus Finkemeier, “Short and longdistance effects in the decay tau → pi tau-neutrino(gamma),” Nucl. Phys. B438 , 17–53 (1995), arXiv:hep-ph/9403385 [hep-ph].[4] Markus Finkemeier, “Radiative corrections to pi(l2) andK(l2) decays,” , Phys. Lett.
B387 , 391–394 (1996), arXiv:hep-ph/9505434 [hep-ph].[5] Vincenzo Cirigliano and Helmut Neufeld, “A note onisospin violation in Pl2(gamma) decays,” Phys. Lett.
B700 , 7–10 (2011), arXiv:1102.0563 [hep-ph].[6] M. Tanabashi et al. (Particle Data Group), “Review ofParticle Physics,” Phys. Rev. D , 030001 (2018).[7] C. T. H. Davies et al. (HPQCD, UKQCD, MILC,Fermilab Lattice), “High precision lattice QCD con-fronts experiment,” Phys. Rev. Lett. , 022001 (2004),arXiv:hep-lat/0304004 [hep-lat].[8] S. Aoki et al. (Flavour Lattice Averaging Group),“FLAG Review 2019,” (2019), arXiv:1902.08191 [hep-lat].[9] N. Carrasco et al. , “Leptonic decay constants f K , f D , and f D s with N f = 2+1+1 twisted-mass lattice QCD,”Phys. Rev. D91 , 054507 (2015), arXiv:1411.7908 [hep-lat].[10] R. Frezzotti and G. C. Rossi, “Twisted mass latticeQCD with mass nondegenerate quarks,”
Lattice hadronphysics. Proceedings, 2nd Topical Workshop, LHP 2003,Cairns, Australia, July 22-30, 2003 , Nucl. Phys. Proc.Suppl. , 193–202 (2004), [,193(2003)], arXiv:hep-lat/0311008 [hep-lat].[11] R. Frezzotti and G. C. Rossi, “Chirally improving Wil-son fermions. 1. O(a) improvement,” JHEP , 007(2004), arXiv:hep-lat/0306014 [hep-lat].[12] R. J. Dowdall, C. T. H. Davies, G. P. Lepage, andC. McNeile, “Vus from pi and K decay constants in fulllattice QCD with physical u, d, s and c quarks,” Phys.Rev. D88 , 074504 (2013), arXiv:1303.1670 [hep-lat].[13] A. Bazavov et al. , “ B - and D -meson leptonic decay con-stants from four-flavor lattice QCD,” Phys. Rev. D98 ,074512 (2018), arXiv:1712.09262 [hep-lat].[14] E. Follana, Q. Mason, C. Davies, K. Hornbostel,G. P. Lepage, J. Shigemitsu, H. Trottier, andK. Wong (HPQCD, UKQCD), “Highly improved stag-gered quarks on the lattice, with applications to charmphysics,” Phys. Rev.
D75 , 054502 (2007), arXiv:hep-lat/0610092 [hep-lat].[15] E. Follana, C.T.H. Davies, G.P. Lepage, andJ. Shigemitsu (HPQCD, UKQCD), “High Precisiondetermination of the pi, K, D and D(s) decay con-stants from lattice QCD,” Phys. Rev. Lett. , 062002(2008), arXiv:0706.1726 [hep-lat].[16] A. Bazavov et al. (MILC), “Results for light pseu-doscalar mesons,” PoS
LATTICE2010 , 074 (2010),arXiv:1012.0868 [hep-lat].[17] S. Durr, Z. Fodor, C. Hoelbling, S.D. Katz, S. Krieg, T. Kurth, L. Lellouch, T. Lippert, A. Ramos, and K.K.Szabo, “The ratio FK/Fpi in QCD,” Phys. Rev. D ,054507 (2010), arXiv:1001.4692 [hep-lat].[18] T. Blum et al. (RBC, UKQCD), “Domain wall QCDwith physical quark masses,” Phys. Rev. D , 074505(2016), arXiv:1411.7017 [hep-lat].[19] Stephan D¨urr et al. , “Leptonic decay-constant ratio f K /f π from lattice QCD using 2+1 clover-improvedfermion flavors with 2-HEX smearing,” Phys. Rev. D , 054513 (2017), arXiv:1601.05998 [hep-lat].[20] V.G. Bornyakov, R. Horsley, Y. Nakamura, H. Perlt,D. Pleiter, P.E.L. Rakow, G. Schierholz, A. Schiller,H. St¨uben, and J.M. Zanotti (QCDSF–UKQCD),“Flavour breaking effects in the pseudoscalar mesondecay constants,” Phys. Lett. B , 366–373 (2017),arXiv:1612.04798 [hep-lat].[21] B. Blossier et al. (ETM), “Pseudoscalar decay constantsof kaon and D-mesons from N f = 2 twisted mass LatticeQCD,” JHEP , 043 (2009), arXiv:0904.0954 [hep-lat].[22] E. Berkowitz, C. Bouchard, C.C Chang, M.A. Clark,B Joo, T. Kurth, C. Monahan, A. Nicholson,K. Orginos, E. Rinaldi, P. Vranas, and A. Walker-Loud,“M¨obius Domain-Wall fermions on gradient-flowed dy-namical HISQ ensembles,” Phys. Rev. D96 , 054513(2017), arXiv:1701.07559 [hep-lat].[23] Tanmoy Bhattacharya, Vincenzo Cirigliano, Saul Co-hen, Rajan Gupta, Anosh Joseph, Huey-Wen Lin, andBoram Yoon (PNDME), “Iso-vector and Iso-scalar Ten-sor Charges of the Nucleon from Lattice QCD,” Phys.Rev.
D92 , 094511 (2015), arXiv:1506.06411 [hep-lat].[24] M.A. Clark, R. Babich, K. Barros, R.C. Brower, andC. Rebbi, “Solving Lattice QCD systems of equa-tions using mixed precision solvers on GPUs,” Com-put.Phys.Commun. , 1517–1528 (2010), https://github.com/lattice/quda , arXiv:0911.3191 [hep-lat].[25] R. Babich, M.A. Clark, B. Joo, G. Shi, R.C. Brower, et al. , “Scaling Lattice QCD beyond 100 GPUs,”(2011), arXiv:1109.2935 [hep-lat].[26] K. Symanzik, “Continuum Limit and Improved Actionin Lattice Theories. 1. Principles and phi**4 Theory,”Nucl. Phys.
B226 , 187–204 (1983).[27] K. Symanzik, “Continuum Limit and Improved Actionin Lattice Theories. 2. O(N) Nonlinear Sigma Modelin Perturbation Theory,” Nucl. Phys.
B226 , 205–227(1983).[28] Dru Bryant Renner, W. Schroers, R. Edwards,George Tamminga Fleming, Ph. Hagler, John W.Negele, K. Orginos, A. V. Pochinski, and D. Richards(LHP), “Hadronic physics with domain-wall valenceand improved staggered sea quarks,”
Lattice field the-ory. Proceedings, 22nd International Symposium, Lat-tice 2004, Batavia, USA, June 21-26, 2004 , Nucl.Phys. Proc. Suppl. , 255–260 (2005), [,255(2004)],arXiv:hep-lat/0409130 [hep-lat].[29] Richard C. Brower, Hartmut Neff, and KostasOrginos, “Mobius fermions: Improved domain wall chi-ral fermions,”
Lattice field theory. Proceedings, 22nd In-ternational Symposium, Lattice 2004, Batavia, USA,June 21-26, 2004 , Nucl. Phys. Proc. Suppl. , 686–688 (2005), [,686(2004)], arXiv:hep-lat/0409118 [hep-lat].[30] R. C. Brower, H. Neff, and K. Orginos, “Mobiusfermions,”
Hadron physics, proceedings of the Work-shop on Computational Hadron Physics, University of Cyprus, Nicosia, Cyprus, 14-17 September 2005 , Nucl.Phys. Proc. Suppl. , 191–198 (2006), arXiv:hep-lat/0511031 [hep-lat].[31] Richard C. Brower, Harmut Neff, and Kostas Orginos,“The m¨obius domain wall fermion algorithm,” (2012),arXiv:1206.5214 [hep-lat].[32] R. Narayanan and H. Neuberger, “Infinite N phase tran-sitions in continuum Wilson loop operators,” JHEP , 064 (2006), arXiv:hep-th/0601210 [hep-th].[33] M. L¨uscher and P. Weisz, “Perturbative analysis ofthe gradient flow in non-abelian gauge theories,” JHEP , 051 (2011), arXiv:1101.0963 [hep-th].[34] Martin L¨uscher, “Chiral symmetry and the Yang–Mills gradient flow,” JHEP , 123 (2013),arXiv:1302.5246 [hep-lat].[35] Martin L¨uscher, “Properties and uses of the Wil-son flow in lattice QCD,” JHEP , 071 (2010),arXiv:1006.4518 [hep-lat].[36] Robert Lohmayer and Herbert Neuberger, “Continuoussmearing of Wilson Loops,” PoS
LATTICE2011 , 249(2011), arXiv:1110.3522 [hep-lat].[37] A. Nicholson, E. Berkowitz, H. Monge-Camacho,D. Brantley, N. Garron, C. C. Chang, E. Rinaldi, M. A.Clark, B. Jo´o, T. Kurth, B. Tiburzi, P. Vranas, andA. Walker-Loud, “Heavy physics contributions to neu-trinoless double beta decay from QCD,” Phys. Rev.Lett. , 172501 (2018), arXiv:1805.02634 [nucl-th].[38] E. Berkowitz, C. Bouchard, D. B. Brantley, C.C Chang,M.A. Clark, N. Garron, B Jo´o, T. Kurth, C. Mon-ahan, H. Monge-Camacho, A. Nicholson, K. Orginos,E. Rinaldi, P. Vranas, and A. Walker-Loud, “An accu-rate calculation of the nucleon axial charge with latticeQCD,” (2017), arXiv:1704.01114 [hep-lat].[39] C. C. Chang, A. Nicholson, E. Rinaldi, E. Berkowitz,N. Garron, D. A. Brantley, H. Monge-Camacho, C. J.Monahan, C. Bouchard, M. A. Clark, B. Jo´o, T. Kurth,K. Orginos, P. Vranas, and A. Walker-Loud, “A per-cent-level determination of the nucleon axial couplingfrom quantum chromodynamics,” Nature , 91–94(2018), arXiv:1805.12130 [hep-lat].[40] Evan Berkowitz, M. A. Clark, A. Gambhir, K. McEl-vain, A. Nicholson, E. Rinalid, P. Vranas, A. Walker-Loud, C. C. Chang, B. Jo´o, T. Kurth, and K. Orginos,“Simulating the weak death of the neutron in afemtoscale universe with near-Exascale computing,”(2018), arXiv:1810.01609 [hep-lat].[41] Oliver Bar, Gautam Rupak, and Noam Shoresh, “Sim-ulations with different lattice Dirac operators for va-lence and sea quarks,” Phys. Rev.
D67 , 114505 (2003),arXiv:hep-lat/0210050 [hep-lat].[42] A. Bazavov et al. (MILC), “Scaling studies of QCD withthe dynamical HISQ action,” Phys. Rev.
D82 , 074501(2010), arXiv:1004.0342 [hep-lat].[43] A. Bazavov et al. (MILC), “Lattice QCD ensembles withfour flavors of highly improved staggered quarks,” Phys.Rev.
D87 , 054505 (2013), arXiv:1212.4768 [hep-lat].[44] David B. Kaplan, “A Method for simulating chiralfermions on the lattice,” Phys. Lett.
B288 , 342–347(1992), arXiv:hep-lat/9206013 [hep-lat].[45] Yigal Shamir, “Chiral fermions from lattice bound-aries,” Nucl. Phys.
B406 , 90–106 (1993), arXiv:hep-lat/9303005 [hep-lat].[46] Vadim Furman and Yigal Shamir, “Axial symmetries inlattice QCD with Kaplan fermions,” Nucl. Phys.
B439 , 54–78 (1995), arXiv:hep-lat/9405004 [hep-lat].[47] Paul Langacker and Heinz Pagels, “Chiral perturbationtheory,” Phys. Rev. D8 , 4595–4619 (1973).[48] J. Gasser and H. Leutwyler, “Chiral Perturbation The-ory to One Loop,” Annals Phys. , 142 (1984).[49] H. Leutwyler, “On the foundations of chiral pertur-bation theory,” Annals Phys. , 165–203 (1994),arXiv:hep-ph/9311274 [hep-ph].[50] Stephen R. Sharpe and Robert L. Singleton, Jr,“Spontaneous flavor and parity breaking with Wilsonfermions,” Phys. Rev. D58 , 074501 (1998), arXiv:hep-lat/9804028 [hep-lat].[51] Oliver Bar, Gautam Rupak, and Noam Shoresh, “Chi-ral perturbation theory at O ( a ) for lattice QCD,”Phys. Rev. D70 , 034508 (2004), arXiv:hep-lat/0306021[hep-lat].[52] Oliver Bar, Claude Bernard, Gautam Rupak, andNoam Shoresh, “Chiral perturbation theory for stag-gered sea quarks and Ginsparg-Wilson valence quarks,”Phys. Rev.
D72 , 054502 (2005), arXiv:hep-lat/0503009[hep-lat].[53] Brian C. Tiburzi, “Baryons with Ginsparg-Wilsonquarks in a staggered sea,” Phys. Rev.
D72 , 094501(2005), [Erratum: Phys. Rev.D79,039904(2009)],arXiv:hep-lat/0508019 [hep-lat].[54] Jiunn-Wei Chen, Donal O’Connell, Ruth S. Van de Wa-ter, and Andre Walker-Loud, “Ginsparg-Wilson pionsscattering on a staggered sea,” Phys. Rev.
D73 , 074510(2006), arXiv:hep-lat/0510024 [hep-lat].[55] Jiunn-Wei Chen, Donal O’Connell, and Andre Walker-Loud, “Two Meson Systems with Ginsparg-WilsonValence Quarks,” Phys. Rev.
D75 , 054501 (2007),arXiv:hep-lat/0611003 [hep-lat].[56] Kostas Orginos and Andre Walker-Loud, “Mixedmeson masses with domain-wall valence and stag-gered sea fermions,” Phys. Rev.
D77 , 094505 (2008),arXiv:0705.0572 [hep-lat].[57] Fu-Jiun Jiang, “Mixed Action Lattice Spacing Effectson the Nucleon Axial Charge,” (2007), arXiv:hep-lat/0703012 [hep-lat].[58] Jiunn-Wei Chen, Donal O’Connell, and Andre Walker-Loud, “Universality of mixed action extrapolation for-mulae,” JHEP , 090 (2009), arXiv:0706.0035 [hep-lat].[59] Jiunn-Wei Chen, Maarten Golterman, Donal O’Connell,and Andre Walker-Loud, “Mixed Action Effective FieldTheory: An Addendum,” Phys. Rev. D79 , 117502(2009), arXiv:0905.2566 [hep-lat].[60] T. Blum et al. , “Quenched lattice QCD with domainwall fermions and the chiral limit,” Phys. Rev.
D69 ,074502 (2004), arXiv:hep-lat/0007038 [hep-lat].[61] Y. Aoki et al. , “Domain wall fermions with im-proved gauge actions,” Phys. Rev.
D69 , 074504 (2004),arXiv:hep-lat/0211023 [hep-lat].[62] Sergey Syritsyn and John W. Negele, “Oscillatory termsin the domain wall transfer matrix,”
Proceedings, 25thInternational Symposium on Lattice field theory (Lattice2007): Regensburg, Germany, July 30-August 4, 2007 ,PoS
LAT2007 , 078 (2007), arXiv:0710.0425 [hep-lat].[63] Szabolcs Borsanyi et al. , “High-precision scale settingin lattice QCD,” JHEP , 010 (2012), arXiv:1203.4469[hep-lat].[64] Gabriel Amoros, Johan Bijnens, and P. Talavera, “Twopoint functions at two loops in three flavor chiral per- turbation theory,” Nucl. Phys. B , 319–363 (2000),arXiv:hep-ph/9907264.[65] B. Ananthanarayan, Johan Bijnens, and Shayan Ghosh,“An Analytic Analysis of the Pion Decay Constantin Three-Flavoured Chiral Perturbation Theory,” Eur.Phys. J. C , 497 (2017), arXiv:1703.00141 [hep-ph].[66] B. Ananthanarayan, Johan Bijnens, Samuel Friot, andShayan Ghosh, “Analytic representations of m K , F K , m η and F η in two loop SU (3) chiral perturbation the-ory,” Phys. Rev. D , 114004 (2018), arXiv:1804.06072[hep-ph].[67] B. Ananthanarayan, Johan Bijnens, Samuel Friot, andShayan Ghosh, “Analytic representation of F K /F π intwo loop chiral perturbation theory,” Phys. Rev. D97 ,091502 (2018), arXiv:1711.11328 [hep-ph].[68] J. Gasser and H. Leutwyler, “Chiral Perturbation The-ory: Expansions in the Mass of the Strange Quark,”Nucl.Phys.
B250 , 465 (1985).[69] Johan Bijnens, “CHIRON: a package for ChPT numer-ical results at two loops,” Eur. Phys. J.
C75 , 27 (2015),arXiv:1412.0887 [hep-ph].[70] S. R. Beane, P. F. Bedaque, K. Orginos, and M. J.Savage (NPLQCD), “ f K /f π in Full QCD with DomainWall Valence Quarks,” Phys. Rev. D75 , 094501 (2007),arXiv:hep-lat/0606023.[71] Jiunn-Wei Chen and Martin J. Savage, “Baryons in par-tially quenched chiral perturbation theory,” Phys. Rev.
D65 , 094001 (2002), arXiv:hep-lat/0111050 [hep-lat].[72] J. Gasser and H. Leutwyler, “Light Quarks at Low Tem-peratures,” Phys. Lett.
B184 , 83–88 (1987).[73] Gilberto Colangelo and Christoph Haefeli, “An Asymp-totic formula for the pion decay constant in a large vol-ume,” Phys. Lett. B , 258–264 (2004), arXiv:hep-lat/0403025.[74] Gilberto Colangelo, Stephan Durr, and Christoph Hae-feli, “Finite volume effects for meson masses and de-cay constants,” Nucl. Phys.
B721 , 136–174 (2005),arXiv:hep-lat/0503014 [hep-lat].[75] Johan Bijnens and Thomas R¨ossler, “Finite Volume atTwo-loops in Chiral Perturbation Theory,” JHEP ,034 (2015), arXiv:1411.6384 [hep-lat].[76] Johan Bijnens and Thomas R¨ossler, “Finite Volume forThree-Flavour Partially Quenched Chiral PerturbationTheory through NNLO in the Meson Sector,” JHEP ,097 (2015), arXiv:1508.07238 [hep-lat].[77] Johan Bijnens and Gerhard Ecker, “Mesonic low-energyconstants,” Ann. Rev. Nucl. Part. Sci. , 149–174(2014), arXiv:1405.6488 [hep-ph].[78] S. Aoki et al. , “Review of lattice results concerning low-energy particle physics,” (2016), arXiv:1607.00299 [hep-lat].[79] Robert E. Kass and Adrian E. Raftery, “Bayes factors,”Journal of the American Statistical Association , 773–795 (1995).[80] Peter Lepage, “gplepage/lsqfit: lsqfit version 11.5.1,”(2020), https://github.com/gplepage/lsqfit .[81] Johan Bijnens, Niclas Danielsson, and Timo A. Lahde,“The Pseudoscalar meson mass to two loops in three-flavor partially quenched chiPT,” Phys. Rev. D ,111503 (2004), arXiv:hep-lat/0406017.[82] Johan Bijnens and Timo A. Lahde, “Decay constants ofpseudoscalar mesons to two loops in three-flavor par-tially quenched (chi)PT,” Phys. Rev. D , 094502(2005), arXiv:hep-lat/0501014. [83] Johan Bijnens, Niclas Danielsson, and Timo A. Lahde,“Three-flavor partially quenched chiral perturbationtheory at NNLO for meson masses and decay con-stants,” Phys. Rev. D , 074509 (2006), arXiv:hep-lat/0602003.[84] D. Giusti, V. Lubicz, G. Martinelli, C.T. Sachrajda,F. Sanfilippo, S. Simula, N. Tantalo, and C. Tarantino,“First lattice calculation of the QED corrections to lep-tonic decay rates,” Phys. Rev. Lett. , 072001 (2018),arXiv:1711.06537 [hep-lat].[85] G.M. de Divitiis, R. Frezzotti, V. Lubicz, G. Martinelli,R. Petronzio, G.C. Rossi, F. Sanfilippo, S. Simula,and N. Tantalo (RM123), “Leading isospin breaking ef-fects on the lattice,” Phys. Rev. D , 114505 (2013),arXiv:1303.4896 [hep-lat].[86] D. Giusti, V. Lubicz, C. Tarantino, G. Martinelli,F. Sanfilippo, S. Simula, and N. Tantalo, “Leadingisospin-breaking corrections to pion, kaon and charmed-meson masses with Twisted-Mass fermions,” Phys. Rev.D , 114504 (2017), arXiv:1704.06561 [hep-lat].[87] A. Bazavov et al. (Fermilab Lattice, MILC), “Charmedand light pseudoscalar meson decay constants from four-flavor lattice QCD with physical light quarks,” Phys.Rev. D90 , 074509 (2014), arXiv:1407.3772 [hep-lat].[88] Tanmoy Bhattacharya, Vincenzo Cirigliano, Saul Co-hen, Rajan Gupta, Huey-Wen Lin, and Boram Yoon,“Axial, Scalar and Tensor Charges of the Nucleon from2+1+1-flavor Lattice QCD,” Phys. Rev. D , 054508(2016), arXiv:1606.07049 [hep-lat].[89] Rajan Gupta, Yong-Chull Jang, Boram Yoon, Huey-Wen Lin, Vincenzo Cirigliano, and Tanmoy Bhat-tacharya, “Isovector Charges of the Nucleon from2+1+1-flavor Lattice QCD,” Phys. Rev. D , 034503(2018), arXiv:1806.09006 [hep-lat].[90] Yong-Chull Jang, Rajan Gupta, Boram Yoon, and Tan-moy Bhattacharya, “Axial Vector Form Factors fromLattice QCD that Satisfy the PCAC Relation,” Phys.Rev. Lett. , 072002 (2020), arXiv:1905.06470 [hep-lat].[91] Andrzej Czarnecki, William J. Marciano, and Al-berto Sirlin, “Neutron Lifetime and Axial CouplingConnection,” Phys. Rev. Lett. , 202002 (2018),arXiv:1802.01804 [hep-ph].[92] Peter Lepage, “gplepage/gvar: gvar version 11.2,”(2020), https://github.com/gplepage/gvar .[93] “LALIBE: CalLat Collaboration code for lattice qcd cal-culations,” https://github.com/callat-qcd/lalibe .[94] Robert G. Edwards and Balint Joo (SciDAC Col-laboration, LHPC Collaboration, UKQCD Collabora-tion), “The Chroma software system for lattice QCD,”Nucl.Phys.Proc.Suppl. , 832 (2005), arXiv:hep-lat/0409003 [hep-lat].[95] The HDF Group, “Hierarchical Data Format, version5,” (1997-NNNN), .[96] Thorsten Kurth, Andrew Pochinsky, Abhinav Sarje,Sergey Syritsyn, and Andre Walker-Loud, “High-Performance I/O: HDF5 for Lattice QCD,” Pro-ceedings, 32nd International Symposium on LatticeField Theory (Lattice 2014): Brookhaven, NY, USA,June 23-28, 2014 , PoS
LATTICE2014 , 045 (2015),arXiv:1501.06992 [hep-lat].[97] Evan Berkowitz, “METAQ: Bundle SupercomputingTasks,” (2017), https://github.com/evanberkowitz/metaq , arXiv:1702.06122 [physics.comp-ph]. [98] Evan Berkowitz, Gustav R. Jansen, Kenneth McEl-vain, and Andr´e Walker-Loud, “Job Management andTask Bundling,” EPJ Web Conf. , 09007 (2018),arXiv:1710.01986 [hep-lat].[99] Chia Cheng Chang, Christopher K¨orber, and Andr´eWalker-Loud, “EspressoDB: A Scientific Database forManaging High-Performance Computing Workflow,” J.Open Source Softw. , 2007 (2020), arXiv:1912.03580[hep-lat].[100] “MILC Collaboration code for lattice qcd calculations,” https://github.com/milc-qcd/milc_qcd .[101] Ulli Wolff (ALPHA), “Monte Carlo errors with less er- rors,” Comput. Phys. Commun. , 143–153 (2004),[Erratum: Comput.Phys.Commun. 176, 383 (2007)],arXiv:hep-lat/0306017.[102] Barbara De Palma, Marco Erba, Luca Mantovani,and Nicola Mosco, “A Python program for the imple-mentation of the Γ -method for Monte Carlo simula-tions,” Comput. Phys. Commun. , 294–301 (2019),arXiv:1703.02766 [hep-lat].[103] Peter J. Moran and Derek B. Leinweber, “Over-improved stout-link smearing,” Phys. Rev. D77