$F_K / F_\pi$ from M\"{o}bius Domain-Wall fermions solved on gradient-flowed HISQ ensembles

We report the results of a lattice QCD calculation of $F_K/F_\pi$ using M\"{o}bius Domain-Wall fermions computed on gradient-flowed $N_f=2+1+1$ HISQ ensembles. The calculation is performed with five values of the pion mass ranging from $130 \lesssim m_\pi \lesssim 400$ MeV, four lattice spacings of $a\sim 0.15, 0.12, 0.09$ and $0.06$ fm and multiple values of the lattice volume. The interpolation/extrapolation to the physical pion and kaon mass point, the continuum, and infinite volume limits are performed with a variety of different extrapolation functions utilizing both the relevant mixed-action effective field theory expressions as well as discretization-enhanced continuum chiral perturbation theory formulas. We find that the $a\sim0.06$ fm ensemble is helpful, but not necessary to achieve a sub-percent determination of $F_K/F_\pi$. We also include an estimate of the strong isospin breaking corrections and arrive at a final result of $F_{K^\pm}/F_{\pi^\pm} = 1.1942(45)$ with all sources of statistical and systematic uncertainty included. This is consistent with the FLAG average value, providing an important benchmark for our lattice action. Combining our result with experimental measurements of the pion and kaon leptonic decays leads to a determination of $|V_{us}|/|V_{ud}| = 0.2311(10)$.

Leptonic decays of the charged pions and kaons provide a means for probing flavor-changing interactions of the Standard Model (SM). In particular, the SM predicts that the Cabibbo-Kobayashi-Maskawa (CKM) matrix is unitary, providing strict constraints on various sums of the matrix elements. Thus, a violation of these constraints is indicative of new, beyond the SM (BSM) physics. There is a substantial flavor physics program dedicated to searching indirectly for potential violations.
CKM matrix elements may be determined through a combination of experimental leptonic decay widths and theoretical determinations of the meson decay constants. For example, the ratio of the kaon and pion decay constants, F K , F π , respectively, may be related to the ratio of light and strange CKM matrix elements |V us |, |V ud | via [1,2], × 1 + δ EM + δ SU (2) . (1.1) In this expression, l = e, µ, the one-loop radiative QED correction is δ EM [3,4] and δ SU (2) is the strong isospin breaking correction that relates F 2 K /F 2 π in the isospin limit to F 2 K + /F 2 π + that includes m d − m u corrections [5] F 2 Using lattice QCD calculations of the ratio of decay constants in the above expression yields one of the most precise determinations of |V us |/|V ud | [6]. Combining the results obtained through lattice QCD with independent determinations of the CKM matrix elements, such as semileptonic meson decays, provides a means for testing the unitarity of the CKM matrix and obtaining signals of new physics.
F K /F π is a so-called gold-plated quantity [7] for calculating within lattice QCD. This dimensionless ratio skirts the issue of determining a physical scale for the lattices, and gives precise results due to the correlated statistical fluctuations between numerator and denominator, as well as the lack of signal-to-noise issues associated with calculations involving, for instance, nucleons. Lattice QCD calculations of F K /F π are now a mature endeavor, with state-of-the-art calculations determining this quantity consistently with sub-percent precision. The most recent review by the Flavour Lattice Averaging Group (FLAG), which performs global averages of quantities that have been calculated and extrapolated to the physical point by multiple groups, quotes a value of FK + F π + = 1.1932 (19) (1. 2) for N f = 2 + 1 + 1 dynamical quark flavors, including strong-isospin breaking corrections [8]. This average includes calculations derived from two different lattice actions, one [9] with twisted-mass fermions [10,11] and the other two [12,13] with the highly improved staggered quark (HISQ) action [14]. The results obtained using the HISQ action are approximately seven times more precise than those from twisted mass and so the universality of the continuum limit for F K /F π from N f = 2 + 1 + 1 results has not been tested with precision yet: in the continuum limit, all lattice actions should reduce to a single universal limit, that of SM QCD, provided all systematics are properly accounted for. Thus, in addition to lending more confidence to its global average, the calculation of a gold-plated quantity also allows for precise testing of new lattice actions, and the demonstration of control over systematic uncertainties for a given action. FLAG also reports averages for N F = 2 + 1, F K ± /F π ± = 1.1917 (37) from Refs. [15][16][17][18][19][20] and for N f = 2, F K ± /F π ± = 1.1205 (18) from Refs. [21], though we restrict our direct comparisons to the N f = 2 + 1 + 1 results just for simplicity.
In this work, we report a new determination of F K /F π calculated with Möbius Domain-Wall fermions computed on 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 and kaon mass extrapolation (χ), continuum limit (a), infinite volume limit (V ), physical point (phys) and model selection (M ) uncertainties, is In the following sections we will discuss details of our lattice calculation, including a brief synopsis of the action and ensembles used, as well as our strategy for extracting the relevant quantities from correlation functions. We will then detail our procedure for extrapolating to the physical point via combined continuum, infinite volume, and physical pion and kaon mass limits and the resulting uncertainty breakdown. We discuss the impact of the a ∼ 0.06 fm ensemble on our analysis, the convergence of the SU (3)-flavor chiral expansion, and the estimate of the strong isospin breaking corrections. We conclude with an estimate of the impact our result has on improving the extraction of |V us |/|V ud | and an outlook.

II. DETAILS OF THE LATTICE CALCULATION
A. MDWF on gradient-flowed HISQ There are many choices for discretizing QCD, with each choice being commonly referred to as a lattice action. These actions correspond to different UV theories that share a common low-energy theory, QCD. Sufficiently close to the continuum limit, the discrete lattice actions can be expanded as a series of local operators known as the Symanzik Expansion [26,27], the TABLE I. Input parameters for our lattice action. The abbreviated ensemble name [23] indicates the approximate lattice spacing in fm and pion mass in MeV. The S, L, XL which come after an ensemble name denote a relatively small, large and extra-large volume with respect to mπL = 4. a Additional ensembles generated by CalLat using the MILC code. The m350 and m400 ensembles were made on the Vulcan supercomputer at LLNL while the a15m135XL, a09m135, and a06m310L ensembles were made on the Sierra and Lassen supercomputers at LLNL and the Summit supercomputer at OLCF using QUDA [24,25]. These configurations are available to any intersted party upon request, and will be available for easy anonymous downloading-hopefully soon.
low-energy Effective Field Theory (EFT) for the discrete lattice action. The Symanzik EFT contains a series of operators having higher dimension than those in QCD, multiplied by appropriate powers of the lattice spacing, a. For all lattice actions, the only operators of massdimension ≤ 4 are those of QCD, such that the explicit effects from the various discretizations are encoded only in higher-dimensional operators which are all irrelevant in the renormalization sense. There is a universality of the continuum limit, a → 0, in that all lattice actions, if calculated using sufficiently small lattice spacing, will recover the target theory of QCD, provided there are no surprises from non-perturbative effects.
Performing LQCD calculations with different actions is therefore valuable to test this universality, to help ensure a given action is not accidentally in a different phase of QCD, and to protect against unknown systematic uncertainties arising from a particular calculation with a particular action. In this work, we use a mixed-action [28] in which the discretization scheme for the valence-quarks is the Möbius Domain-Wall Fermion (MDWF) action [29][30][31] while the discretization scheme for the sea-quarks is the Highly Improved Staggered Quark action [14]. Before solving the MDWF propagators, we apply a gradientflow [32][33][34] smoothing algorithm [35,36] to the gluons to dampen UV fluctuations, which also significantly improves the chiral symmetry properties of the MDWF action [22] (for example, the residual chiral symmetry breaking scale of domain-wall fermions m res is held to less than 10% of m l for reasonable values of L 5 and M 5 , see Tab. I). Our motivation to perform this calculation is to improve our understanding of F K /F π and to test the MDWF on gradient-flowed HISQ action we have used to compute the π − → π + neutrinoless double beta decay matrix elements arising from prospective higherdimension lepton-number-violating physics [37], and the axial coupling of the nucleon g A [38,39]. As there is an otherwise straightforward path to determining g A to sub-percent precision with pre-exascale computing such as Summit at OLCF and Lassen at LLNL [40], it is important to ensure this action is consistent with known results at this level of precision.
There are several motivations for chosing this mixedaction (MA) scheme [28,41]. The MILC Collaboration provides their gauge configurations to any interested party and we have made heavy use of them. They have generated the configurations covering a large parameter space allowing one to fully control the physical pion mass, infinite volume and continuum limit extrapolations [42,43]. The good chiral symmetry properties of the Domain Wall (DW) action [44][45][46] significantly suppress sources of chiral symmetry breaking from any sea-quark action, motivating the use of this mixed-action setup. While this action is not unitary at finite lattice spacing, we have tuned the valence quark masses such that the valence pion mass matches the taste-5 HISQ pion mass within a few percent, so as the continuum limit is taken, we recover a unitary theory.
EFT can be used to understand the salient features of such MALQCD calculations. Chiral Perturbation Theory (χPT) [47][48][49] can be extended to incorporate discretization effects into the analytic formula describing the quark-mass dependence of various hadronic quantities [50]. The MAEFT [51] for DW valence fermions on dynamical rooted staggered fermions is well developed [52][53][54][55][56][57][58][59]. The use of valence fermions which respect chiral symmetry leads to a universal form of the MAEFT extrapolation formulae at NLO (next-to-leading order) in the joint quark mass and lattice spacing expansions [55,58], which follows from the suppression of chiral symmetry breaking discretization effects.

B. Correlation function construction and analysis
The correlation function construction and analysis follows closely the strategy of Ref. [22] and [38,39]. Here we summarize the relevant details for this work.
The pesudoscalar decay constants F can be obtained from standard two-point correlation functions by making use of the 5D Ward-Takahashi Identity [60,61] where q 1 and q 2 denote the quark content of the meson with lattice input masses m q1 and m q2 respectively. The point-sink ground-state overlap-factor z 0p and groundstate energy E 0 are extracted from a two-point correlation function analysis with the model 2) where n encompasses in general an infinite tower of states, t is the source-sink time separation, T is the temporal box size and we have both smeared (s) and point (p) correlation functions which both come from smeared sources. From Ref. [22], we show that gradientflow smearing leads to the suppression of the domainwall fermion oscillating mode (which also decouples as M 5 → 1, at least in free-field [62]), and therefore this mode is not included in the correlator fit model. Finally, the residual chiral symmetry breaking m res is calculated by the ratio of two-point correlation functions evaluated at the midpoint of the fifth dimension L 5 /2 and bounded on the domain wall [31] m res (t) = x π(t, x, L 5 /2)π(0, 0, 0) x π(t, x, 0)π(0, 0, 0) , where π(t, x, s) ≡q(t, x, w)γ 5 q(t, x, w) is the pseudoscalar interpolating operator at time t, space x and fifth dimension s. We extract m res by fitting Eq. (2.3) to a constant.

Analysis strategy
For all two-point correlation function parameters (MDWF and mixed MDWF-HISQ), we infer posterior parameter distributions in a Bayesian framework using a 4-state model which simultaneously describes the smeared-and point-sink two-point correlation functions (the source is always smeared). The joint posterior distribution is approximated by a multivariate normal distribution (we later refer to this procedure as fitting). The two-point correlation functions are folded in time to double the statistics. The analysis of the pion, kaon,sγ 5 s, and mixed MDWF-HISQ mesons are performed independently, with correlations accounted for under bootstrap resampling.
We analyze data of source-sink time separations between 0.72 fm to 3.6 fm for all 0.09 fm and 0.12 fm lattice spacing two-point correlation functions, and 0.75 fm to 3.6 fm for all 0.15 fm lattice spacing two-point correlation functions.
We choose normally distributed priors for the groundstate energy and all overlap factors, and log-normal distributions for excited-state energy priors. The groundstate energy and overlap factors are motivated by the plateau values of the effective masses and scaled correlation function, and a prior width of 10% of the central value. The excited-state energy splittings are set to the value of two pion masses with a width allowing for fluctuations down to one pion mass within one standard deviation. The excited-state overlap factors are set to zero, with a width set to the mean value of the ground-state overlap factor.
Additionally, we fit a constant to the correlation functions in Eq. (2.3). For the 0.09 fm to 0.12 fm ensembles, we analyze source-sink separations that are greater than 0.72 fm. For the 0.12 fm ensemble, the minimum sourcesink separation is 0.75 fm. The prior distribution for the residual chiral symmetry breaking parameter is set to the observed value per ensemble, with a width that is 100% of the central value. The uncertainty is propagated with bootstrap resampling.
We emphasize that all input fit parameters (i.e. number of states, fit region, priors) are chosen to have the same values in physical units for all observables, to the extent that a discretized lattice allows. Additionally, we note that the extracted ground-state observables from these correlation functions are insensitive to variations around the chosen set of input fit parameters.

III. EXTRAPOLATION FUNCTIONS
We now turn to the extrapolation/interpolation to the physical point. We have three ensembles at the physical 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 mass extrapolation is an interpolation. Nevertheless, we explore how the ensembles with heavier pion masses impact the physical point prediction and we use our dataset to explore uncertainty arising in the SU (3)-flavor chiral expansion.
We begin by assuming a canonical power counting scheme for our MA LQCD action [52] in which O(m 2 π ) ∼ O(m 2 K ) ∼ O(a 2 Λ 4 QCD ) are all treated as small scales. For the quark mass expansion, the dimensionless small parameters (m P /4πF ) 2 naturally emerge from χPT where P ∈ {π, K, η}. For the discretization corrections, while  (17) aΛ 2 QCD is often used to estimate the relative size of corrections compared to typical hadronic mass scales, it is a bit unatural to use this in a low-energy EFT as Λ QCD is a QCD scale that does not emerge in χPT.
We chose to use another hadronic scale to form a dimensionless parameter with a, that being the gradient flow scale w 0 ∼ 0.17 fm [63]. This quantity is easy to compute, has mild quark mass dependence, and the value is roughly w −1 0 4πF π . We then define the dimensionless small parameters for controlling the expansion to be We leave F ambiguous, as we will explore taking F = F π , F = F K and F 2 = F π F K in our definition of Λ χ . This particular choice of a is chosen such that the range of values of this small parameter roughly correspond to 2 π 2 a 2 K as the lattice spacing is varied, similar to the variation of 2 π itself over the range of pion masses used, see Tab. II. As we will discuss in Sec. IV, this choice of a seems natural as determined by the size of the 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 orders in the expansion are defined to be  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( 2 a ) correction as it must accompany terms which vanish in the SU (3) V limit, such as 2 K − 2 π . Therefore, at NLO, there can not be any counterterms proportional to 2 a and the only discretization effects that can appear at NLO come through modification of the various meson masses that appear in the MA EFT.
We find that the precision of our results requires including terms higher than NLO, and we have to work at a hybrid N 3 LO order to obtain a good description of our data. Therefore, we will begin with a discussion of the full N 2 LO χPT theory expression for F K /F π in the continuum limit [64][65][66][67].
The analytic expression for F K /F π up to N 2 LO is [67] The first line is the LO (1) plus NLO terms, while the next three lines are the N 2 LO terms. Several non-unique choices were made to arrive at this formula. Prior to discussing these choices, we first define the parameters appearing in Eq. (3.3). First, the small parameters were all defined as where F π (m P ) is the "on-shell" pion decay constant at the masses m P . The quantities P are defined as where µ is a renormalization scale. The coefficientL 5 = (4π) 2 L r 5 (µ) is one of the regulated Gasser-Leutwyler LECs [68] which has a renormalization scale dependence that exactly cancels against the dependence arising from the logarithms appearing at the same order. In the following, we define all of the Gasser-Leutwyler LECs with the extra (4π) 2 for convenience: The η mass has been defined through the Gell-Mann-Okubo (GMO) relation with the corrections to this relation being propagated into Eq. (3.3) for consistency at N 2 LO. The logs are The ln 2 terms are encapsulated in the F F (x) function, defined in Eqs. (8)(9)(10)(11)(12)(13)(14)(15)(16)(17) of Ref. [67], 1 and theK r i λ P λ P terms whose coefficients are given by 2 The single log coefficientsĈ r 1−3 are combinations of the NLO Gasser-Leutwyler coefficientŝ 1 They also provide an approximate formula which is easy to implement, but our numerical results are sufficiently precise to require the 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 have provided a Python interface as well so that the function can be called from our main analysis code, which is provided with this article. 2 We correct a typographical error in the K r 6 term presented in Ref. [67]: a simple power-counting reveals the ξ 2 K = 4 K accompanying this term should not be there. where Finally,Ĉ r 4 is a combination of these L r i coefficients as well as counterterms appearing at N 2 LO. At N 2 LO, only two counterterm structures can appear due to the SU (3) V constraints: which are linear combinations of the N 2 LO counterterms and contributions from the Gasser-Leutwyler LECs (Eq. (7) of Ref. [67]) There were several non-unique choices that went into the determination of Eq. (3.3). When working with the full N 2 LO χPT expression, the different choices one can make result in different N 3 LO or higher corrections and exploring these different choices in the analysis will expose sensitivity to higher-order contributions that are not explicitly included. The first choice we discuss is the Taylor expansion of the ratio of F K /F π Where the · · · represent higher order terms in the expansion and δF 3) has been derived from this standard Taylor-expanded form with the choices mentioned above: the use of the onshell renormalized value of F → F π and the definition of the η mass through the GMO relation. The NLO expressions are the standard ones [68] δF NLO The δF N 2 LO P terms have been determined in Ref. [64] and cast into analytic forms in Refs. [65,66]. The NLO terms are of O(20%) and so Taylor expanding this ratio leads to sizeable corrections from the δF NLO π 2 − δF NLO π δF NLO K contributions. Utilizing the full ratio expression could in principle lead to a noticeable difference in the analysis (a different determination of the values of the LECs for example). Rather than implementing the full δF N 2 LO P expressions for kaon and pion, we explore this convergence by instead just resumming the NLO terms which will dominate the potential differences in higher order corrections. A consistent expression at N 2 LO is (3.18) and the ratio correction is given by Another choice we explore is the use of F → F π in the definition of the small parameters. Such a choice is very convenient as it allows one to express the small parameters entirely in terms of observables one can determine in the lattice calculation (unlike the bare parameters, such as χPT's F 0 and Bm q , which must be determined through extrapolation analysis). Equally valid, one could have chosen F → F K or F 2 → F π F K . Each choice induces explicit corrections one must account for at N 2 LO to have a consistent expression at this order. The NLO corrections in Eq.
plus higher order corrections.
Related to this choice, Eq. (3.3) is implicitly defined at the standard renormalization scale [67] While F K /F π of course does not depend upon this choice, the numerical values of the LECs do. Further, a scale setting would be required to utilize this or any fixed value of µ. Instead, as was first advocated in Ref. [70] to the best of our knowledge, it is more convenient to set the renormalization scale on each ensemble with a lattice quantity. For example, Ref. [70] used µ = f latt is the lattice-determined value of the pion decay constant on a given ensemble. The advantage of this choice is that the entire extrapolation can be expressed in terms of ratios of lattice quantities such that a scale setting is not required to perform the extrapolation to the physical point.
At NLO in the expansion, one is free to make this choice as the corrections appear at N 2 LO. In the present work, we must account for these corrections for a consistent expression at this order, which is still defined at a fixed renormalization scale. To understand these corrections, we take as our fixed scale where F 0 is the decay constant in the SU (3) chiral limit. Define µ π = 4πF π and consider the NLO expression where we have introduced the notation If we chose the renormalization scale µ π and add the second term of the last equality, then this expression is equivalent to working with the scale µ 0 through N 2 LO. The convenience of this choice becomes clear as µ π /µ 0 has a familiar expansion Using the GMO relation Eq. (3.7) and expanding ln(1+x) for small x, this expression becomes (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 which are made more convenient if one also makes the replacements F 2 π → {F π F K , F 2 K } in the definition of the small parameters plus the corresponding N 2 LO corrections that accompany these choices.
If we temporarily expose the implicit dependence of the expression for F K /F π on the choices of F and µ, such that Eq. (3.3) is defined as then the following expressions are all equivalent at N 2 LO and the LECs in these expressions are related to those at the standard scale by evolving them from µ ρ 0 → µ 0 with their known scale dependence [68]. Implicit in these expressions is the normalization of the small parameters (3.30) We have described several choices one can make in parameterizing the χPT formula for F K /F π . The key point is that if the underlying chiral expansion is well behaved, the formula resulting from each choice are all equivalent through N 2 LO in the SU (3) chiral expansion, with differences only appearing at N 3 LO and beyond. Therefore, by studying the variance in the extrapolated answer upon these choices, one is assessing some of the uncertainty arising from the truncation of the chiral extrapolation formula.

B. Discretization Corrections
We now turn to the discretization corrections. We explore two parametrizations for incorporating the corrections arising at finite lattice spacing. The simplest approach is to use the continuum extrapolation formula and enhance it by adding contributions from all allowed powers of 2 P and 2 a to a given order in the expansion. This is very similar to including only the contributions from local counter-terms that appear at the given order. At N 2 LO, the set of discretization corrections is given by where A 4 s and A 4 α S are the LECs and α S is the running QCD coupling that emerges in the Symanzik expansion of the lattice expansion through loop corrections. Each contribution at this order must vanish in the SU (3) V limit because the discretization corrections are flavor-blind and so we have the limiting constraint lim m l →ms at any lattice spacing. From a purist EFT perspective, we should instead utilize the MA EFT expression. Unfortunately, the MA EFT expression is only known at NLO [52] and our results require higher orders to provide good fits. Nevertheless, we can explore the utility of the MA EFT by replacing the NLO χPT expression with the NLO MA EFT expression and using the continuum expression enhanced with the local discretization corrections at higher orders, Eq. (3.31).
Using the parameterization of the hairpin contributions from Ref. [55], the NLO MA EFT expressions are In these expressions, we use the partially quenched flavor notation [71]

35)
TABLE III. Extracted masses of the mixed MDWF-HISQ mesons. We use the notation from Ref. [71] in which mπ and mK 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 andss mesons to have the same mass as the HISQ sea pion andss mesons within a few percent, the quantities m 2 ju − m 2 π and other splittings provide an estimate of the additive mixed-meson mass splitting due to discretization effects, a 2 ∆Mix [52] and additional additive corrections [59]. At LO in MAEFT, these splittings are predicted to be quark mass independent, which we find to be approximately true, with a notable decrease in the splitting as the valence-quark mass is increased as first observed in Ref. [56] as well as a milder decrease as the seq-quark mass is increased. ensemble amju amjs amru amrs amss where m ju is the mass of a mixed valence-sea pion. The partial quenching parameters δ ju 2 and δ 2 rs provide a measure of the unitarity violation in the theory. For our MDWF on HISQ action, at LO in MA EFT, they are given by the splitting in the quark masses plus a discretization correction arising from the taste-Identity splitting For the tuning we have done, setting the valence-valence pion mass equal to the taste-5 sea-sea pion mass, these parameters are given just by the discretization terms as m u = m j and m s = m r within 1-2%. The sea-sea eta mass in this tuning is given at LO in MA EFT as These parameters, and the corresponding meson masses are provided in Tab. III. The expressions for d π , K φ1φ2 , K (2,1) φ1φ2 and K φ1φ2φ3 are provided in Appendix B.
At NLO in the MA EFT, the LECs which contribute to δF K and δF π are the same as in the continuum, L 4 and L 5 , plus a discretization LEC which we have denoted L a . Just like the L 4 contribution, the contribution from L a exactly cancels in δF K − δF π . At N 2 LO, beyond the continuum counterterm contributions, (3.12), there are the two additional LECs contributions, (3.31).

C. Finite Volume Corrections
We now discuss the corrections arising from the finite spatial volume. The leading FV corrections arise from the tadpole integrals which arise at NLO in both the χPT and MA expressions. The well known modification to the integral can be expressed as [72][73][74]  the finite volume correction to these integrals is The full finite volume corrections to the continuum formula are also known at N 2 LO [75] as well as in the partially quenched χPT [76]. In this work, we restrict the corrections to those arising from the NLO corrections as our results are not sensitive to higher-order FV corrections. This is because, with the ensembles used in this work, all ensembles except a12m220S satisfy m π L 4 (see Tab. II). MILC generated three volumes for this a12m220 ensemble series to study FV corrections. Fig. 1 shows a comparison of the results from the a12m220L, a12m220, and a12m220S along with the predicted volume corrections arising from NLO in χPT. The uncertainty band arises from an N 3 LO fit using the full N 2 LO continuum χPT formula enhanced with discretization LECs and N 3 LO corrections arising from continuum and finite lattice spacing corrections. Even with one of the most precise fits, we see that the numerical results are consistent with the predicted NLO FV corrections.

D. N 3 LO Corrections
The numerical data set in this work requires us to add N 3 LO corrections to obtain a good fit quality. At this order, we only consider local counterterm contributions, of which there are three new continuum like corrections and three discretization corrections. A non-unique, but complete parameterization is In principle, we could also add counterterms proportional to higher powers of α S but with four lattice spacings, we would not be able to resolve the difference between the complete set of operators including all possible additional α S corrections. The set of operators we do include is sufficient to parametrize the approach to the continuum limit.

IV. EXTRAPOLATION DETAILS AND UNCERTAINTY ANALYSIS
We now carry out the extrapolation/interpolation to the physical point, which we perform in a Bayesian framework. To obtain a good fit, we must work to N 3 LO in the mixed chiral and continuum expansion. The results from the a06m310L ensemble drive this need, in particular, for higher-order discretization corrections to parametrize the results from all the ensembles. We will explore the impact of the a06m310L ensemble in more detail in this section. First, we discuss the values of the priors we set and the definition of the physical point.

A. Prior widths for LECs
The number of additional LECs we need to determine at each order in the expansion is order N Li N χ N a NLO 1 0 0 N 2 LO 7 2 2 N 3 LO 0 3 3 Total 8 5 5 . N Li is the number of Gasser-Leutwyler coefficients, N χ the number of counterterms associated with the continuum χPT expansion and N a is the number of counterterms associated with the discretization corrections. In total, there are 18 unknown LECs. While we utilize 18 ensembles in this analysis, the span of parameter space is not sufficient to constrain all the LECs without prior knowledge. In particular, the introduction of all 8 L i coefficients requires prior widths informed from phenomenology.
In the literature, the L i are typically quoted at the renormalization scale µ ρ = 770 MeV while in our work, we use the scale µ F0 = 4πF 0 . We can use the BE14 values of the L i LECs from Ref. [77] and the known scale dependence [68] to convert them from µ ρ to µ F0 :   (Table 3), using their known scale dependence [68], Eq. (4.1). We assign the following slightly more conservative uncertainty as a prior width in the minimization: If a value of Li is less than 0.5 × 10 −3 , we assign it a 100% uncertainty at the scale µ = 770 MeV; If the value is larger than 0.5 × 10 −3 , we assign it the larger of 0.5 or 1/3 of the mean value.  (14) -0.34 (34) 0.29 (47) with the values of Γ i listed in Table V for convenience. We use F 0 = 80 MeV, which is the value adopted by FLAG [8]. We set the central value of all the L i with this procedure and the widths are set as described in Tab. V. Next, we must determine priors for the N 2 LO and N 3 LO local counterterm coefficients, A n K,π,s . We set the central value of all these priors to 0 and then perform a simple grid search varying the widths to find preferred values of the width, as measured by the Bayes Factor. Our goal is not to optimize the width of each prior individually for each model used in the fit, but rather find a set of prior widths that is close to optimal for all models. To this end, we vary the width of the χPT LECs together at each order (N 2 LO, N 3 LO) and the discretization LECs together at each order (N 2 LO, N 3 LO) for a four-parameter search. We apply a very crude grid where the values of the widths are taken to be 2, 5, or 10.
We find taking the width of all these A n K,π,s LECs equal to 2 results in good fits with near-optimal values. This provides evidence the normalization of small parameters we have chosen for 2 P and 2 a , Eq. (3.1), is "natural" and supports the power counting we have assumed, Eq. (3.2). The N 2 LO LECs mostly favor a width of 2 while the N 3 LO discretization LECs prefer 5 and the N 3 LO χPT LECs vary from model to model with 5 a reasonable value for all. As a result of this search, we pick as our priors A 4 K,π = 0 ± 2,Ã 4 s = 0 ± 2, A 6 K,π = 0 ± 5,Ã 6 s = 0 ± 5. (4.2)

B. Physical Point
As our calculation is performed with isospin symmetric configurations and valence quarks, we must define a physical point to quote our final result. We adopt the definition of the physical point from FLAG. FLAG[2017] [78] defines the isospin symmetric pion and kaon masses to be (Eq.  The isospin symmetric physical point is then defined by extrapolating our results to the values (for the choice (4.5)

C. Model Averaging Procedure
Our model average is performed under a Bayesian framework following the procedure described in [39,79]. Suppose we are interested in estimating the posterior distribution of Y = F K /F π , ie. P (Y |D) given our data D.
To that end, we must marginalize over the different models M k .
Here P (Y |M k , D) is the distribution of Y for a given model M k and dataset D, while P (M k |D) is the posterior distribution of M k given D. The latter can be written, per Bayes' theorem, as .

(4.7)
We can be more explicit with what the latter is in the context of our fits. First, mind that we are a priori agnostic in our choice of M k . We thus take the distribution P (M k ) to be uniform over the different models. We calculate P (D|M k ) by marginalizing over the parameters (LECs) in our fits: (4.8) After 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 Bayes Factor as reported by lsqfit [80]. Thus with K the number of models included in our average. We emphasize that this model selection criterium not only rates the quality of the decription of data, it also penalizes parameters which do not improve this description. This helps rule out models which overparameterize data. Now we can estimate the expectation value and variance of Y .
The variance V [Y ] results from the total law of variance; the first term in brackets is known as the expected value of the process variance (which we refer to as the model averaged variance), while the latter is the variance of the hypothetical 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 extrapolation/interpolation to the physical point. The  Based upon the quality of fit (gauged by the Bayesian analogue to the p-value, Q, or the reduced chi square, χ 2 ν ) and/or the weight determined as discussed in the previous section, we can dramatically reduce the number of models used in the final averaging procedure. First, any model which does not include the FV correction from NLO is heavily penalized. This is not surprising given the observed volume dependence on the a12m220 ensembles, Fig. 1. However, even if we remove the a12m220S ensemble from the analysis, the Taylor-expanded fits have a relative weight of e −6 or less compared to those that have χPT form at NLO.
If we add FV corrections to the Taylor expansion fits (pure counterterm) and use all ensembles, c n m π L|n| K 1 (m π L|n|) + · · · (4.12) they still have weights which are ∼ e −8 over the normalized model distribution and also contribute negligibly to the model average. We observe that the fits which use the MA EFT at NLO are also penalized with a relative weight of ∼ e −8 , and fits which only work to N 2 LO have unfavorable weights by ∼ e −5 (and are also accompanied by poor χ 2 ν values). Cutting all of these variations reduces our final set of models to be N 3 LO χPT with the following variations ×3 : use The finite volume uncertainty is assessed by removing the a12m220S ensemble from the analysis, repeating the model averaging procedure and taking the difference. The final probability distribution broken down into the three choices of F 2 is shown in Fig. 2.

Impact of a06m310L ensemble
Next, we turn to understanding the impact of the a06m310L ensemble on our analysis. The biggest dif- Bayes Model Avg PDF Final probability distribution giving rise to Eq. (1.3), separated into the three choices of F 2 = {F 2 π , FπFK , F 2 K } in the definition of the small parameters, Eq. (3.1). The parent "gray" distribution is the final PDF normalized to 1 when integrated.
ference upon removing the a06m310L ensemble is that the data are not able to constrain the various terms contributing to the continuum extrapolation as well, particularly since there are up to three different types of scaling violations: and thus, the statistical uncertainty of the results grows as well as the model variance, with a total uncertainty growth from ∼ 0.0044 to ∼ 0.0057, and the mean of the extrapolated answer moves by approximately half a standard deviation. Furthermore, N 2 LO fits become acceptable, though they are still grossly outweighed by the N 3 LO fits. Including both effects, the final model average result shifts from (57) . • (Right): no a06m310L, N 2 LO χPT with only counterterms at N 2 LO and F = F π .
As can be seen from the middle plot, the a15, a12 and a09 ensembles prefer contributions from both 2 a and 4 a contributions and are perfectly consistent with the result on the a06m310L ensemble. They are also consistent with an N 2 LO fit (no 4 a contributions) as can be seen in the right figure. However, the weight of the N 3 LO fits is still significantly greater than the N 2 LO fits even without the a06m310L data.
We conclude that the a06m310L ensemble is useful, but not necessary to obtain a sub-percent determination of F K /F π with our lattice action. A more exhaustive comparison can be performed with the analysis notebook provided with this publication.
In Fig. 4, we show the stability of our final result for various choices discussed in this section.

Convergence of the chiral expansion
While the numerical analysis favors a fit function in which only counterterms are used at N 2 LO and higher, it is interesting to study the convergence of the chiral expansion by studying the fits which use the full χPT expression at N 2 LO.
In Fig. 5, we show the resulting light quark mass dependence using the N 3 LO extrapolation with the full N 2 LO χPT formula. After the analysis is performed, the results from each ensemble are shifted to the physical kaon mass point, leaving only dependence upon 2 π and 2 a as well as dependence upon the η mass defined by the GMO relation. The magenta band represents the full 68% confidence interval in the continuum, infinite volume limit. The different colored curves are the mean values as a function of 2 π at the four different lattice spacings. We also show the convergence of this fit in the lower panel plot. From this convergence plot, one sees that roughly that at the physical pion mass (vertical gray line) the NLO contributions add a correction of ∼ 0.16 compared to 1 at LO, the N 2 LO contributions add another ∼ 0.04 and the N 3 LO corrections are not detectible by eye. The band at each order represents the sum of all terms up to that order determined from the full fit. The reduction in the uncertainty as the order is increased is due on large part to the induced correlation between the LECs at different orders through the fitting procedure.
In Fig. 2, we observe that the different choices of F are all consistent, indicating higher-order corrections (starting at N 3 LO in the non-counterterm contributions) are smaller than the uncertainty in our results. It is also interesting to note that choosing F πK or F K is penalized by the analysis, indicating the numerical results prefer larger expansion parameters. In Tab. VI, we show the resulting χPT LECs determined in this analysis for the two choices F = {F π , F Kπ }, as well as whether the ratio form of the fit is used, Eq. (3.17). For the Gasser-Leutwyler LECs, we evolve the values back from µ 0 → µ ρ for a simpler comparison with the values quoted in literature. For most of the L i , we observe the numerical results have very little influence on the parameters as they mostly return the prior value (also listed in the table for convenience). The only LECs influenced by the fit are L 5 , L 7 , and L 8 with L 5 getting pulled about one sigma away from the prior value and L 7 and L 8 only shifting by a third or half of the prior width. What is most interesting from these results is that our fit prefers a value of L 5 that is noticeably smaller than the value obtained by MILC [16]   and HPQCD [12] and also smaller than the BE14 result from Ref. [77]. In Fig. 6, we show the impact of using the fully expanded expression, Eq. (3.15) versus the expression in which 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 2 LO χPT expression. We see that fits without the ratio form are preferred, but the central value of the final result depends minimally upon this choice.
In Fig. 7, we show that the results strongly favor the use of only counterterms at N 2 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 fit favors only counterterms at N 2 LO and higher. While the linear 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 penalized for the counterterm redundancy, Eq. (3.12). Using the priors for L i from Tab. V, we find the contribution from the Gasser-Leutwyler LECs to these N 2 LO counterterms, Eq. (3.14) are given by (94) . (4.14) As the A 4 P terms are priored at 0(2), it is sufficient to re-run the analysis by simply setting L 4 P = 0. We find this result marginally improves the Bayes Factors but not statistically significantly, leaving us with the puzzle that the optimal fit is a hybrid NLO χPT plus counterterms at higher orders.
If the Taylor expansion fits (pure counterterm) were good and favored over the χPT fits, this could be a sign that the SU (3) χPT formula was failing to describe the lattice results. However, we have to include the NLO χPT expression, including its predicted (counterterm free) volume dependence to describe the numerical results. It would be nice to have the full N 2 LO MA EFT expression to understand why the hybrid MA EFT fits are so relatively disfavored in the analysis. There may be compensating discretization effects that cancel against those at NLO to some degree that might allow the full N 2 LO MA EFT to better describe the results. It is unlikely, however, that this expression will be derived as at two loops in χPT, the universality of MA EFT expressions [58] breaks down and they can no longer be "derived" from the corresponding partially quenched formula, which is known for F K /F π at two loops [76,[81][82][83].

E. QCD isospin breaking corrections
Finally, we discuss the correction to our result to obtain a direct determination of F K + /F π + including strong isospin breaking corrections, but excluding QED corrections. This is the standard value quoted in the FLAG reviews [8,78]. Our calculation, like most, are performed in the isospin symmetric limit, and therefore, the strong isospin breaking correction must be estimated, rather than having a direct determination. The optimal approach is to incorporate both QED and QCD isospin breaking corrections into the calculations such that the separation is not necessary, as was done in Ref [84] by incorporating both types of corrections through the perturbative modification of the path integral and correlation functions [85,86]. In this work, we have not performed these extensive computations and so we rely upon the SU (3) χPT prediction to estimate the correction due to strong isospin breaking. As we have observed in Sec. IV D, the SU (3) chiral expansion behaves and converges nicely, so we expect this approximation to be reasonable.
The NLO corrections to F K and F π including the strong isospin breaking corrections are given by K ± , (4.15) where we have kept explicit the contribution from each flavor of meson propagating in the loop. There are three points to note in these expressions: 1. At NLO in the SU (3) chiral expansion, there are no additional LECs that describe the isospin breaking corrections beyond those that contribute to the isospin symmetric limit. Therefore, one can make 0.00 0.02 0.04 0.06 0.08 2 π = (m π /4πF π ) 2   The PDFs are taken from the parent PDF, Fig. 2 without renormalizing such that height in this figure reflects the relative weight comopared to the total PDF. PDFs from N2LO ct or full XPT FIG. 7. Comparison of N 3 LO χPT analysis with F = Fπ using the full N 2 LO χPT expression (smaller histogram) versus only counterterms at N 2 LO, Eq. (3.12) and (3.31). As in Fig. 6, the PDFs are drawn from the parent PDF.
a parameter-free prediction of the isospin breaking corrections using lattice results from isospin symmetric 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 free of isospin breaking corrections at this order; 3. We have used the kaon mass splitting in place of the quark mass splitting, which is exact at LO in χPT The estimated shift of our isospin-symmetric result to incorporate strong isospin breaking is then Ref. [5] suggested replacingL 5 with the NLO expression equating it to the isospin symmetric F K /F π which yields In this expression, we have utilized the two relations At this order, both of these expressions, Eqs. (4.16) and (4.17) are equivalent. However, they can result in shifts that differ by more than one standard deviation. Further, the direct estimate of the strong isospin breaking corrections [9] is larger in magnitude than either of them. Therefore, to estimate the strong isospin breaking corrections, we take the larger of the two corrections as the mean and the larger uncertainty of the two, and then add an additional 25% uncertainty for SU (3) truncation errors. In Sec. IV D 2 we observe the N 2 LO correction is ∼ 25% of the NLO correction (while NLO is ∼ 16% of LO). In order to evaluate these expressions, we have to define the physical point with strong isospin breaking and without QED isospin breaking. We employ the values from FLAG[2017] [78] (exceptM π 0 = 134.6(3) MeV): With this definition of the physical point, we find (under the same model average as Tab. VII) 20) resulting in our estimated strong isospin breaking correction

V. SUMMARY AND DISCUSSION
The ratio F K /F π may be used, in combination with experimental input for leptonic decay widths, to make a prediction for the ratio of CKM matrix elements, |V us |/|V ud |. Using the most recent data, Eq. (1.1) becomes [6], Result for the ratio of CKM matrix elements, |Vus|/|V ud |, extracted from the ratio FK /Fπ reported in this work (red band). The global lattice value for |Vus| extracted from a semileptonic decay form factor, f+(0) [8], is shown as a horizontal blue band, while the global experimental average for |V ud | from nuclear beta decay [6] is given as a vertical green band. Note that the intersection between the red and green bands agrees well with the unitarity constraint for the CKM matrix, while the intersection between the red and blue bands shows ∼ 2σ tension.
Alternatively, rather than using the experimental determination of |V ud | as input for our test of unitarity, we may instead use the global lattice average for |V us | = 0.2231(7) [8], extracted via the quantity f + (0), the zero momentum transfer limit of a form factor relevant for the semileptonic decay K 0 → π − lν. This leads to, leading to a roughly 2σ tension with unitarity. Our result, along with with the reported experimental results for |V ud | and lattice results for |V us |, are shown in Fig. 8. One could also combine our results with the more precise average in the FLAG review which would lead to a slight reduction their reported uncertainties, but we will leave that 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 nucleon structure calculations [38][39][40] can be used to reproduce an accepted value from other lattice calculations that are known at the sub-percent level. Our result provides the first sub-percent cross-check of the universality   [87] result (LEFT) with the continuum extrapolation in the present work with the MDWF on gradient-flowed HISQ action (right). In the right plot, we also include the FLAG[2020] average value from the N f = 2 + 1 + 1 calculations [8]. While the a06m310L ensemble is not necessary for us to extrapolate to a consistent value as this FLAG average (see Fig. 3), the overall size of our discretization effects are larger. This is not necessarily surprising as the HISQ action used by MILC has perturbatively removed all O(a 2 ) corrections such that the leading scaling violations begin at O(αSa 2 ), as implied by the x-axis of the left plot.
of the continuum limit of this quantity with N f = 2+1+1 dynamical flavors, albeit with the same sea-quark action as used by MILC/FNAL and HPQCD [12,13].
Critical in obtaining a sub-percent determination of any quantity is control over the continuum extrapolation. This is relevant to our pursuit of a sub-percent determination of g A as another calculation, utilizing many of the same HISQ ensembles but with a different valence action (clover fermions), obtains a result that is in tension with our own [88,89]. While there has been speculation that this discrepancy is due to the continuum extrapolations [89], new work suggests the original work underestimated the systematic uncertainty in the correlation function analysis, and when accounted for, the tension between 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 neutron lifetime discrepancy [91], it is important to understand the scaling violations of our lattice action. While a smooth continuum extrapolation in one observable does not guarantee such a smooth extrapolation in another, it at least provides some reassurance of a well-behaved continuum extrapolation. Furthermore, the determination of F K /F π involves the same axial current that is relevant for the computation of the nucleon matrix element used to compute g A . Fig. 3 shows the continuum extrapolation of F K /F π from our analysis. The size of the discretization effects are noticeably larger than we observed in our calculation of g A [39]. In Sec. IV D 1, we demonstrated that, while helpful, the a06m310L ensemble is not necessary to achieve a sub-percent determination of F K /F π . This is in contrast to the determination by MILC which requires the a ∼ 0.06 fm (or smaller) lattice spacings to control the continuum extrapolation (though we note, the HPQCD calculation [12], also performed on the HISQ ensembles, does not utilize the a ∼ 0.06 fm ensembles but agrees with the MILC result). It should be noted, the MILC result does not rely on the heavier mass ensembles except to adjust for the slight mistuning of the input quark masses on their near-physical point ensembles. In Fig. 9, we compare our continuum extrapolation to that of MILC [87].
In Ref. [87], they also utilize the same four lattice spacings as in this work (they have subsequently improved their determination with an additional two finer lattice spacings [13].) A strong competition between the O(a 2 ) and O(a 4 ) corrections was observed in that work, such that the a ∼ .06 fm ensemble is much more instrumental for a reliable continuum extrapolation than is the case in our setup. At the same time, the overall scale of their discretization effects is much smaller than we observe in the MDWF on gradient-flowed HISQ action for this quantity. This is not entirely surprising as the HISQ action has been tuned to perturbatively remove all O(a 2 ) corrections such that the leading corrections formally begin as O(α S a 2 ).

ACKNOWLEDGMENTS
We would like to thank V. Cirigliano, S. Simula, J. Simone, and T. Kaneko for helpful correspondence and discussion regarding the strong isospin breaking corrections to F K /F π . We would like to thank J. Bijnens for helpful correspondence on χPT and a C++ interface to CHIRON [69] that we used for the analysis presented in this work. We thank the MILC Collaboration for providing some of the HISQ configurations used in this work, and A. Bazavov, C. Detar and D. Toussaint for guidance on using their code to generate the new HISQ ensembles also used in this work. We would like to thank P. Lepage for enhancements to gvar [92] and lsqfit [80] that enable the pickling of lsqfit.nonlinear fit objects.
Computing time for this work was provided through the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program and the LLNL Multiprogrammatic and Institutional Computing program for Grand Challenge allocations on the LLNL su-percomputers. This research utilized the NVIDIA GPUaccelerated Titan and Summit supercomputers at Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725 as well as the Surface, RZHasGPU, Pascal, Lassen, and Sierra supercomputers at 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] and status of tasks logged with EspressoDB [99]. The HMC was performed with the MILC Code [100], and for the ensembles 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 analysis and data for this work can be found at this git repo: https://github.com/callat-qcd/project_fkfpi.  which has been regulated and renormalized with the standard χPT modified dimensional-regularization scheme [48]. The finite volume corrections to δ π are given by The expression for K φ1φ2 arises from the integral Similarly, K (2,1) φ1φ2 is given by Finally, K φ1φ2φ3 is given by In each of these expressions, the corresponding expression including FV corrections are given by replacing φ → FV φ , Eq. (3.38).  Table IV of Ref. [43]), we list the abbreviated ensemble name, number of streams Nstream and total number of configurations N cfg . For a given ensemble, each stream has an equal number of configurations. The gauge coupling, light, strange, and charm quark masses on each ensembles are given as well as the tadpole factor, u0 and the Naik-term added to the charm quark action, N . The total length in molecular dynamics time units (MDTU) between each saved configuration is s while the length between accept/reject steps is Len. The microstep size used in the HMC is provided as Len./Nsteps which was input with single precision. The average acceptance rate over all streams is listed as well as the number of streams. history of the ∆S for the three ensembles. For the a15m135XL ensemble, we reduced the trajectory length significantly compared to the a15m130 from MILC to overcome spikes in the HMC force calculations. To compensate, we lowered the acceptance rate to encourage the HMC to move around parameter space with larger jumps, to hopefully reduce the autocorrelation time. We ran 25 HMC accept/reject steps before saving a configuration for a total trajectory length of 5. For each accept/reject step we also measure the quark-anti-quark condensateψψ using a stochastic estimate with 5 random sources that are averaged together. We compute it for each of the quark masses am l , am s , and am c . On  the a15m135XL we have measuredψψ only on every saved configuration for the first half of each stream, while we measured it at each accept/reject step for the second half. The integrated autocorrelation time, as well as the average and statistical errors ofψψ are computed using the Γ-method analysis [101] with the Python package unew [102]. We report the results in Tab. IX. In Fig. 11 we report the value of theψψ on each saved configuration for the three quark masses on each ensemble. Because we observe a long autocorrelation time of the ψ s ψ s on the a15m135XL ensemble, we also studied the uncertainty on the extracted pion and kaon effective masses as a function of block size to check for possible longer autocorrelations than usual, with blocking lengths of 10, 25, and 100 MDTU. We observe that these hadronic quantities have a much shorter autocorrelation time as the uncertainty is independent of τ b and consistent with the unblocked data. On this a15m135XL ensemble, while we have generated 2000 configurations, we have only utilized 1000 in this paper (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 the a06m310L ensemble is less than ideal and the intergrated autocorrelation time is long, the volume is sufficiently large    [103] with weight parameters ρ given in Tab. X and ε = −0.25. We cross-checked a sample of our stout smeared measurements with the more expensive Symanzik flow technique and saw good agreement between the two. We determined the ρ parameter and the number of steps to perform on an ensemble-toensemble basis, i.e., for a handful of configurations per ensemble we choose a spread of ρs and step numbers and observe which combination gives the best plateau. These values of ρ and step number (ε is always -0.25) are then applied to the entire ensemble.