$|V_{us}|$ from $K_{\ell 3}$ decay and four-flavor lattice QCD

Using HISQ $N_f=2+1+1$ MILC ensembles with five different values of the lattice spacing, including four ensembles with physical quark masses, we have performed the most precise computation to date of the $K\to\pi\ell\nu$ vector form factor at zero momentum transfer, $f_+^{K^0\pi^-}(0)=0.9696(15)_\text{stat}(11)_\text{syst}$. This is the first calculation that includes the dominant finite-volume effects, as calculated in chiral perturbation theory at next-to-leading order. Our result for the form factor provides a direct determination of the Cabibbo-Kobayashi-Maskawa matrix element $|V_{us}|=0.22333(43)_{f_+(0)}(42)_\text{exp}$, with a theory error that is, for the first time, at the same level as the experimental error. The uncertainty of the semileptonic determination is now similar to that from leptonic decays and the ratio $f_{K^+}/f_{\pi^+}$, which uses $|V_{ud}|$ as input. Our value of $|V_{us}|$ is in tension at the 2--$2.6\sigma$ level both with the determinations from leptonic decays and with the unitarity of the CKM matrix. In the test of CKM unitarity in the first row, the current limiting factor is the error in $|V_{ud}|$, although a recent determination of the nucleus-independent radiative corrections to superallowed nuclear $\beta$ decays could reduce the $|V_{ud}|^2$ uncertainty nearly to that of $|V_{us}|^2$. Alternative unitarity tests using only kaon decays, for which improvements in the theory and experimental inputs are likely in the next few years, reveal similar tensions. As part of our analysis, we calculated the correction to $f_+^{K\pi}(0)$ due to nonequilibrated topological charge at leading order in chiral perturbation theory, for both the full-QCD and the partially-quenched cases. We also obtain the combination of low-energy constants in the chiral effective Lagrangian $[C_{12}^r+C_{34}^r-(L_5^r)^2](M_\rho)=(2.92\pm0.30)\cdot10^{-6}$.


Abstract
Using HISQ N f = 2 + 1 + 1 MILC ensembles with five different values of the lattice spacing, including four ensembles with physical quark masses, we have performed the most precise computation to date of the K → π ν vector form factor at zero momentum transfer, f K 0 π − + (0) = 0.9696 (15) stat (11) syst . This is the first calculation that includes the dominant finitevolume effects, as calculated in chiral perturbation theory at next-to-leading order. Our result for the form factor provides a direct determination of the Cabibbo-Kobayashi-Maskawa matrix element |V us | = 0.22333 (43) f + (0) (42) exp , with a theory error that is, for the first time, at the same level as the experimental error. The uncertainty of the semileptonic determination is now similar to that from leptonic decays and the ratio f K + /f π + , which uses |V ud | as input. Our value of |V us | is in tension at the 2-2.6σ level both with the determinations from leptonic decays and with the unitarity of the CKM matrix. In the test of CKM unitarity in the first row, the current limiting factor is the error in |V ud |, although a recent determination of the nucleus-independent radiative corrections to superallowed nuclear β decays could reduce the |V ud | 2 uncertainty nearly to that of |V us | 2 . Alternative unitarity tests using only kaon decays, for which improvements in the theory and experimental inputs are likely in the next few years, reveal similar tensions and could be further improved by taking correlations between the theory inputs. As part of our analysis, we calculated the correction to f Kπ + (0) due to nonequilibrated topological charge at leading order in chiral perturbation theory, for both the full-QCD and the partially-quenched cases. We also obtain the combination of low-energy constants in the chiral effective Lagrangian [C r 12 + C r 34 − (L r 5 ) 2 ](M ρ ) = (2.92 ± 0.30) · 10 −6 .

I. INTRODUCTION
High-precision tests of the unitarity of the Cabibbo-Kobayashi-Maskawa (CKM) matrix, as predicted by the Standard Model (SM), are at the forefront of the current flavor physics program. Any violation of the unitarity of the CKM matrix, which describes flavor-changing interactions, would be evidence of the existence of physics beyond the Standard Model (BSM).
In particular, first-row unitarity, which requires that vanishes, is currently the most precisely tested condition. Even in the absence of deviations, high-precision determinations of the CKM matrix elements involved in the test in Eq. (1.1) put important constraints on the scale of the allowed new physics [1]. At the current level of precision one can neglect |V ub | 2 in Eq. (1.1). Of the other two CKM matrix elements involved, |V ud | is precisely determined from super-allowed nuclear β decays [2]. It can also be extracted from measurements of the neutron lifetime [3] and pion β decay [4], albeit with much larger errors [5]. Improved experimental measurements of these processes would be interesting because they are theoretically cleaner.
The best determinations of |V us | are from kaon decays [6]. The extraction of |V us | from semileptonic kaon (K 3 ) decay requires knowledge of the form factor at zero momentum transfer, f Kπ + (0), which is still the largest source of uncertainty on |V us |. On the experimental side, it is expected that the ongoing and forthcoming experiments (NA62, OKA, KLOE-2, LHCb and TREK E36) could reduce the experimental error to ∼ 0.12% within 5 years [7]. Reducing the theoretical error in the vector form factor calculation is therefore a crucial task: it is this task that we take up in this paper.
Determinations of |V us | from leptonic kaon and pion decays (K 2 and π 2 ), combined with f K /f π from lattice QCD, currently have somewhat smaller errors than those from K 3 . The total error in |V us | from leptonic decays is 0.25% [8][9][10][11][12][13][14][15][16], while from semileptonic decays it is 0.34% [6]. These leptonic determinations are indirect, however, because they require an external input for |V ud |, namely Ref. [2]. The direct extraction of |V us | from only kaon leptonic decays using f K as nonperturbative input gives a larger error of 0.46%. 1 Currently, the value of |V us | obtained from leptonic kaon decay is ∼ 2σ larger than the value from semileptonic kaon decay [16]. The leptonic decay is mediated by the axial-vector current while the semileptonic decay by the vector current. According to the SM, both approaches should give the same |V us |, because the W boson current has a V − A structure. Thus, any significant difference should be carefully analyzed.
In addition, if |V ud | is taken from Ref. [2], the leptonic value of |V us | is consistent with unitarity, Eq. (1.1), but the semileptonic value of |V us | leads to a ∼ 2σ disagreement with unitarity. As we were finishing this work, a paper appeared with a new calculation of the nucleus-independent electroweak radiative corrections involved in the extraction of |V ud | from superallowed β decays with a new approach based on dispersion relations [17]. If this calculation is confirmed, the resulting value of |V ud | would increase the present tension with unitarity. Investigating the origin of these tensions and performing even more stringent tests is crucial for the internal consistency of the Standard Model. It is thus necessary to reduce the error on both the experimental and the lattice-QCD inputs entering determinations of |V us |.
In this work we reduce the main sources of uncertainty in our previous calculation of f K 0 π − + (q 2 = 0) to reach a total error of 0.19%, obtaining the most precise calculation to date, and matching the current experimental uncertainty, for the first time. The main improvements over our previous calculation [24,25] are increased statistics in some key ensembles, the addition of a new (smaller) lattice spacing, and the correction of finite-volume effects at next-to-leading order (NLO) in chiral perturbation theory (ChPT). Preliminary results were presented in Refs. [32,33].
There are other ways to determine |V us |. Semileptonic hyperon decays unfortunately lack sufficiently precise knowledge of the SU(3)-breaking corrections, which precludes a competitive determination. A conservative estimate of such effects yields an uncertainty of ∼ 2% [34]. Inclusive hadronic τ decays have, in the past, yielded values of |V us | smaller than the semileptonic kaon determination and, thus, were in even more disagreement with unitarity [35]. A more recent analysis [36] uses lattice QCD to compute dimension-largerthan-4 condensates and, more importantly, employs a dispersive technique to obtain the Kπ branching fractions. It points to an inclusive-τ value of |V us | compatible with unitarity [37], although it still remains on the low side. An even more promising approach also based on inclusive strange hadronic τ decay data is presented in Ref. [38]. Its basic ingredients are replacing the operator-product expansion in the relevant sum rules by lattice hadronic vacuum polarization functions and optimizing the weight functions to suppress contributions from the high-energy region, where the experimental data have poor precision. Preliminary results in Ref. [38] are compatible with both semileptonic and leptonic determinations (and thus with unitarity), but have larger errors than either. Because the total errors on |V us | in Refs. [36][37][38] are dominated by experimental uncertainties, it is expected that these determinations will be significantly improved with new data from the Belle II experiment [39]. Determinations of |V us | from exclusive τ decays, which use the same nonperturbative inputs as the leptonic kaon decay determinations, namely f K ± and f K ± /f π ± , are still in tension with unitarity [40], but this could also change with future experimental measurements. This paper is organized as follows. In Sec. II we describe the methodology of the numerical lattice-QCD simulations and the details of the ensembles, actions, and correlation functions used. Section III shows how, following the ChPT approaches of Refs. [41] and [42], one can correct for leading-order finite-volume effects and for the effects of nonequilibrated topological charge on the ensemble with finest lattice spacing (a ≈ 0.042 fm), respectively. We discuss the joint chiral interpolation and continuum extrapolation of our data to the physical point in Sec. IV. Section V analyzes the statistical and systematic uncertainties. Final results for the form factor f K 0 π − + , as well as for the relevant O(p 6 ) low energy constants, are presented in Sec. VI. In Sec. VII, we use our form factor result to extract a value of |V us | from kaon semileptonic experimental data and discuss the implications of this value for phenomenology. Finally, we present our conclusions and the prospects for further improvement in Sec. VIII.

II. LATTICE SIMULATIONS AND CORRELATOR FITS
The methodology in this work largely follows that of our previous work in Refs. [24,25,28]. The approach, pioneered by HPQCD [43], is based on the Ward-Takahashi identity relating the matrix elements of a vector current to that of the corresponding scalar density: with Z V and Z S the lattice renormalization factors for the vector current and scalar density, respectively. Working with staggered fermions, and choosing V lat µ to be the partially conserved, taste singlet, vector current, and S lat to be its divergence, we have Z V = Z S = 1. Thus, S lat is a local, taste-singlet density, with the same flavor content as the vector current, S =su. With the above identity and the definition of the form factors in Eq. (1.4), one can extract the scalar form factor f 0 (q 2 ) at any value of the momentum transfer q 2 by using In addition, a kinematic constraint requires f Kπ + (0) = f Kπ 0 (0), so this relation can be employed to calculate f Kπ + (0) from 3-point correlation functions with a scalar insertion. As already discussed in our first work [28], the use of a local scalar density instead of a vector current has two main advantages: avoiding the use of a renormalization factor and avoiding the use of noisier correlation functions with either a non-local vector current or external non-Goldstone mesons [43,44].
We perform our calculation on the HISQ N f = 2 + 1 + 1 MILC configurations [45][46][47] with sea quarks simulated with the highly-improved staggered quark (HISQ) action [48]. We also employ the HISQ action for the valence quarks. We have already seen in our previous work that the use of the HISQ action greatly reduces discretization effects [24,25]. The charm-quark and strange-quark masses on the N f = 2 + 1 + 1 MILC configurations are always tuned to values close to the physical ones, while the light-quark masses vary between 0.2m s and m s /27, with the latter approximately the physical value. In this work, we include data generated at 5 different values of the lattice spacing down to a ≈ 0.042 fm, with sea pion masses ranging from 319 MeV to 134 MeV. Table I lists the key parameters of the ensembles analyzed here and the correlation functions calculated on them. The ensembles include four with physical quark masses and a ≈ 0.15, 0.12, 0.09, 0.06 fm. Ensembles that are new since our analysis in Ref. [24,25] are marked with a dagger in the last column; those where we have increased the statistics are marked with an asterisk. Table I also lists the pseudoscalar-taste (physical) pion mass m π,P and the root-mean-squared pion mass m RMS π for each ensemble. The difference is a measure of the dominant discretization effects, which arise from taste-changing interactions. As expected, they decrease rapidly as the lattice spacing is reduced. The data included in this analysis are graphically depicted in Fig. 1.
The structure of the three-point function with a scalar insertion that we generate to access the matrix element in Eq. (2.2) is the same as in our previous work [24,25,28]. We generate light quarks at a time slice t src and extended strange propagators at a fixed distance T from the source. For each configuration we have N src time sources placed at t src = t 0 , t 0 + L t /N src , t 0 + 2L t /N src . . . , where L t is the temporal length of the lattice. The time t 0 varies randomly from configuration to configuration in an interval [0, L t /N src ] to reduce autocorrelations. Roughly following Ref. [49] we use random-wall sources at the pion source time t src . On that spatial time-slice, we choose four stochastic color-vector fields from a Gaussian distribution, with support on all three colors, and compute light quark propagators from each of the four sources. Between the source and the kaon sink at time t src + T , we contract the extended strange propagator with a light propagator to form the scalar density. We then study the t dependence to isolate the desired matrix element.
The light-quark masses are always the same in the sea and valence sectors, while the Parameters of the N f = 2 + 1 + 1 gauge-field ensembles used in this work, and details of the correlation functions generated. A dagger at the end of a row indicates an ensemble that is new since our work in Ref. [24]; an asterisk indicates that the statistics have been increased. N conf is the number of configurations included in the analysis, N src the number of time sources used on each configuration, and L the spatial size of the lattice. The m π values are in MeV, with m π,P the Goldstone (pseudoscalar taste) π mass, and m RMS π the root-mean-squared (over all tastes) π mass. The ensemble with a ≈ 0.12 fm, m l /m s = 0.1 and m π,P L = 3.2 is used solely for the study of finite volume effects.  sea and valence strange-quark masses are slightly different in some of the ensembles 3 ; see Table I. Table I also lists the number of configurations and time sources on each ensemble. We compute three-point functions as described above for 4 or 5 different values of the source-sink separation T , corresponding to the same physical distances across ensembles. We include both even and odd values of T to disentangle the effects from oscillating states in the correlation functions. We simulate directly at zero momentum transfer, q 2 ≈ 0, by tuning the external momentum of the pion using partially twisted boundary conditions. In particular, we tune with θ 2 the twist angle of the daughter propagator going from the pion to the current. The rest of the propagators are generated with periodic boundary conditions, the same as in the sea sector. We always have diagonal twist angles, θ 2 = | θ 2 |(1, 1, 1)/ √ 3, which turn out to give smaller finite volume effects than twisting in only one direction [41]. The values of | θ 2 | for each ensemble, as well as the corresponding momentum of the pion, are given in Table II.  Table I). The area of each disk is proportional to the statistical sample size N conf × N src . Ensembles on which we have increased the statistics or we have added since our earlier work in Ref. [25] are indicated with black outlines. The three disks with a ≈ 0.12 fm and m π ≈ 200 MeV correspond to the three ensembles with m l /m sea s = 0.1 and different volumes (smaller to larger from top to bottom) in Table I. The yellow disk (smallest volume) is not included in the final analysis but used only to study finite-volume effects.
For each ensemble, we generate zero-momentum two-point π and K correlation functions and two-point π correlation functions with external momentum given by the twist angle θ 2 , defined in Eq. (2.3). These correlators are given by where the interpolating operator Φ p † P ( x, t) creates a meson P = π, K at time t with momentum p. The random wall sources automatically implement the sum over x. We also generate three-point correlation functions with the kaon at rest: where the pion recoil momentum p π is either equal to zero or set to the values listed in Table II. The scalar density is a local taste-singlet.

A. Fitting method and statistical error
The fitting strategy we follow to extract the physical quantities from our correlation functions has already been discussed in Refs. [28,51]. We fit the two-point correlation functions for a pseudoscalar meson P to the expression using Bayesian techniques. In Eq. (2.6), L t is the temporal size of the lattice. The oscillating terms with (−1) m(t+1) do not appear for a zero-momentum π. We fit the three-point correlation functions to where the pion and kaon energies and amplitudes, E n π , E n K , Z π n and Z K n , are the same as those appearing in the two-point fit functions.
We first fit the two-point functions one by one on each ensemble and check the stability of the ground state masses/energies and amplitudes under the choice of fitting range, t ∈ [t min , t max ], and the number of exponentials included in the fit function [N exp in Eq. (2.6)]. We always include the same number N of regular and oscillating states in those fits, i.e., N exp = 2N . We observe that for most of the choices of time range and for all ensembles, fits stabilize when including 2+2 or 3+3 states. From that parameter-scanning procedure, we select an optimal set of fit parameters for the two-point functions, with a common t min for all the functions on the same ensemble. Values of t min correspond to very similar physical distances, approximately 0.6-0.7 fm, on each ensemble.
We then use those t min values to perform a combined Bayesian fit including the two-and three-point functions needed to extract f Kπ 0 (0): π two-point correlation functions with and without momentum, K two-point correlation functions without momentum, N T three-point correlation functions with q 2 ≈ 0, and N T values of the source-sink separation T . In most of the cases, including more than three three-point correlators with different T values does not reduce the error but makes it more difficult to get a stable fit. This is always true for fits including more than four different values of T .
In general, we use three-point data in the combined fits with t ∈ [t min , T −t min ], where t min is the value optimized for the two-point functions. However, on some ensembles, especially the largest ones, we either shorten the three-point fit range or thin the three-point data in order to obtain a good fit.
From the combined fits, we extract the scalar form factor at zero momentum transfer via where A 00 (0) is the ground state three-point parameter in Eq. (2.7), the meson masses and energies are the values extracted from the combined fits and m s and m l are the valence strange and light-quark masses simulated.
We check the stability of the combined fit results under the variation of fit ranges, number of states, and number and values of source-sink separations included, and chose a preferred fit for each ensemble so the shift on the central value with those variations is well under the statistical error, and the error is also stable.
, on each ensemble-see Sec. III A for details of the calculation of ∆ V f + (0). The quark masses am l and am sea s are the same as in Table I, and the twisting angle θ 2 is the one for the light daughter propagator in the pion, defined in Eq. (2.3). The superscript P in the pion masses refers to the pseudoscalar taste. The errors in f + (0) are statistical only. They are generated with a 500-bootstrap distribution. We study the effect of autocorrelations by blocking the data with increasing numbers of successive configurations and redoing the analysis. For ensembles where we observe significant changes in central value and error for the form factor with blocking, those effects stabilize when blocking by four. We thus choose to block by four in all ensembles to avoid autocorrelations.
In Table II and Fig. 2, we show the raw results for the vector form factor at zero momentum transfer from the combined fits described above. The statistical errors shown in the table and the figure as a function of am l /am physical s come from 500 bootstrap resamples and range from 0.16% to 0.38%. In the figure, one can see that for a fixed value of the light quark mass, m l , the points with different shapes, which correspond to different values of the lattice spacing, lie on top of each other, with the exception of the data point for the 0.15 fm ensemble with physical quark masses (in the left-most cluster of points). This is the only ensemble where we observe statistically significant discretization effects.

III. FORM-FACTOR CORRECTIONS
Before performing the chiral-continuum fit, we correct the form factor results listed in Table II and shown in Fig. 2 for the leading-order finite-volume effects and the nonequilibrated topological charge in our finest ensemble. These corrections are described in the following two subsections.

A. Finite volume
In this work we use NLO ChPT to correct our form factor results for finite volume effects, whereas in our previous calculation [24,25] we simply estimated the associated systematic error from a comparison of the lattice data at two different spatial volumes, with other parameters held fixed. The partially twisted boundary conditions used in our calculation introduce several complications in the analysis. In particular, an extra form factor h µ is required to parametrize the weak-current matrix element in finite volume: The three form factors depend on the choice of twisting angles, as well as the value of q 2 . We apply the one-loop formulae in Ref. [41] in the staggered partially-twisted partiallyquenched case to all ensembles included in our calculation for the choices of twist angles in Table II. Because we are calculating the vector form factor at zero momentum transfer via the relation in Eq. (2.2), and the quantity we obtain at finite volume on the lattice is , the FV correction to our results is given by where the meson masses in the denominators are the ones from the simulations, and quantities that are second order in the finite volume corrections have been neglected. In Eq. (3.2), the FV correction for a given quantity X, Notice that, since we extract the meson masses from correlation functions where all the propagators have zero momentum, the FV corrections to the meson masses should be calculated from the formulae in Ref. [41] with zero twisting angles. The resulting FV corrections are listed in the last column of Table II. We find that they are ≤ 0.1% on all ensembles. Some of the values for ∆ V f + (0) are particularly small due to the cancellation between the two contributions in Eq. (3.2). We subtract the ∆ V f + (0) from the finite-volume f Kπ + (0) (listed in the next-to-last column in Table II) before performing the chiral-continuum fit discussed in Section IV.

B. Nonequilibrated topological charge
The HISQ N f = 2 + 1 + 1 MILC simulations with smallest lattice spacings have reached a regime where the distribution of the topological charge Q is not properly sampled [42,52], which affects the physical observables calculated on those ensembles. The issue is relevant here for the ensemble with the finest lattice spacing, a ≈ 0.042 fm. On the other hand, the topological charge is reasonably well equilibrated on the other ensembles, which have a 0.06 fm.
In order to correct for this systematic effect, one can use ChPT to study the Q-dependence of a given observable [42,53,54]. The recent ChPT study in Ref. [42] has already been applied to the calculation of heavy-light meson decay constants and masses in Refs. [47,55]. Here, we extend the analysis of Ref. [42] to f Kπ + (0). The three-point correlation functions relevant for this study, as well as any meson mass calculated in a finite volume V and at fixed Q, satisfy [53,54] where B on the right-hand side is the infinite-volume value of the quantity of interest averaged over Q, B is its second derivative with respect to the vacuum angle θ, evaluated at θ = 0, and χ T = lim V →∞ Q 2 /V is the infinite-volume topological susceptibility. Knowing the dependence on Q or, equivalently, on θ, one can calculate the appropriate correction to B to account for the difference between the correct Q 2 and the simulation Q 2 sample . With Eq. (3.3), we follow Ref. [42] to calculate the correction as where f Kπ + (0) sample is the simulation value. Although we extract f Kπ + (0) from the scalar-density matrix element in Eq. (2.2) it is simpler to first calculate the θ dependence of the vector-current matrix element directly. In ChPT, the vector current with the relevant flavor is where Σ is the SU(3) chiral matrix. In the presence of θ, and for the m u = m d = m l and full QCD case (the case relevant for this work), the O(p 2 ) ChPT Lagrangian is where f is the chiral-limit value of the meson decay constant, and µ the low energy constant that relates meson and quark masses at leading order (LO)-see Eq. (4.3). Here M A ≡ e iθ/3 M, with M the usual quark mass matrix in the absence of θ. When θ = 0, Σ gets the vacuum expectation value The parameter α encodes the dependence on θ, with α(θ = 0) = 0. The relation between α and θ is obtained by minimizing the potential energy term in the Lagrangian, which gives the condition For the expansion of the relevant observables, one needs α , the first derivative of α with respect to θ evaluated at θ = 0. Equation (3.8) implies [42] One may expand Σ around its vacuum expectation value via with Φ the 3 × 3 matrix of meson fields. With this result substituted into Eq. (3.5), at tree level there are two possible diagrams, shown in Fig. 3, that may contribute to the matrix element π|V µ |K . The strong three-point vertex in the right-hand diagram is forbidden by parity when θ = 0, but here comes from the mass term in Eq. (3.6), which violates parity symmetry unless one also takes θ → −θ (which is called "extended parity"). The weak vertex in the right-hand diagram generates a factor of q µ , implying that that diagram contributes only to the form factor f − . From the left-hand diagram, one finds Finally, from Eq. (3.9), the result needed to adjust the form factor via Eq.
Because we actually use Eq. (2. 2) to calculate f Kπ + (0), it is important to check that we can reproduce Eq. (3.12) by calculating the matrix element of the scalar density that appears in the θ = 0 Ward identity at q 2 = 0, : Diagrams contributing to π|ρ|K at tree level, where ρ is either the vector current V µ or the scalar densityS. The squares are weak vertices with the insertion of the current or density and the black dot is a strong vertex.
, and λ f ∈ SU (N f ) the appropriate flavor matrix to select thesu current. Note that it is no longer convenient to take out a factor of m s − m l fromS, as we do for S in Eq. (2.2), because the quark masses now carry factors of exp(±iθ/3). In ChPT, (3.14) Evaluating the diagrams in Fig. 3, we find where the contributions in the first and second lines come from the propagator and vertex diagrams in Fig. 3, respectively. The Ward identity in Eq. (3.13) is then satisfied trivially at LO for θ = 0. For θ = 0, we may calculate f Kπ + (0) from Eq. (3.13) using Eq. (3.9), the fact that α = 0 (which follows from extended parity), and the second derivatives from [42]. After some algebra, we find that the result agrees with Eq. (3.12), as expected.
We have also checked analytically that the Ward identity holds for arbitrary q and θ.
Following the procedure in Ref. [42] for the partially quenched case, we may generalize Eq. (3.12) to where x and y are the active valence quarks (the valence up and strange for f K 0 π − + (0)), and m l and m s are the light and strange sea quark masses. In deriving Eq. (3.17), we have set the spectator quark mass (the d quark mass for f K 0 π − + ) equal to the light sea mass m l ; in other words, the spectator quark is unitary, not partially quenched. This has allowed us to avoid analyzing the case of three partially quenched quarks, which was not treated in Ref. [42]. Since the mass of the spectator quark does not affect f K 0 π − + (0) to LO, we believe Eq. (3.17) will remain valid even when the spectator quark is partially quenched. As expected, Eq. (3.17) reduces to Eq. (3.12) when m x = m l and m y = m s .
As discussed at the beginning of this section, the correction is only needed on the finest ensemble included in this analysis, with a ≈ 0.042 fm. We calculate the correction of Eq. (3.4) using Eq. (3.12), the value of the average of the topological charge measured on that ensemble, Q 2 sample = 27.59 [42], and the correct Q 2 as estimated by the LO ChPT expression for the topological susceptibility [70,71] where the singlet meson sea masses are defined in Eq. (4.3), below. The resulting correction Table II before performing the chiral-continuum fit.

IV. CHIRAL-CONTINUUM INTERPOLATION/EXTRAPOLATION
We follow a methodology very similar to that in our previous analyses [24,25] in order to combine our simulation data into physical results in the continuum limit and with the correct quark/meson masses. Here we summarize the main ingredients and then discuss in more detail the new features added in order to accurately account for finite-volume and isospin-breaking corrections. Accounting for these effects turns out to be essential, given the improvements in the simulation data.
Our methodology is developed in the framework of chiral perturbation theory (ChPT), which allows us to incorporate effects due to mass dependence, discretization, finite-volume, and isospin-breaking in a systematic way. In particular, in the isospin limit, we can write f Kπ + (0) as a chiral expansion where the functions f i are chiral corrections of O(p i ). The Ademollo-Gatto (AG) theorem [56] ensures that the vector form factor goes to one in the limit m s → m u , and that corrections to this limit are second order. That means that the functions f i are proportional to (m s − m u ) 2 or, equivalently, (m 2 K − m 2 π ) 2 . The theorem thus implies that, in the continuum, the O(p 2 ) (one-loop) contribution, f 2 , is completely fixed in terms of experimental quantities: the decay constant f π and meson masses.
The specific fit function we employ for the extrapolation to the continuum and interpolation to the physical quark masses is the same as in Ref. [24]. It consists of an NLO partially quenched staggered ChPT (PQSChPT) expression [57] f PQSChPT 2 (a), plus NNLO continuum ChPT terms [58] f cont 4 , plus extra analytic terms to parametrize higher order discretization  [46,61]. The absolute scale r 1 is from Ref. [62]. The value of the decay constant f π is taken from Ref. [3]; its error, though shown, is negligible in our calculation. Taste splittings r 2 1 a 2 ∆ Ξ are taken from Ref. [46] and more recent updates; slopes aµ come from the analysis presented in Ref. [63], alhough they were not published there. We do not consider errors either on the taste splittings or on the slopes because they also have a negligible effect on the final results. Notice that taste splittings for the a ≈ 0.042 fm ensemble are not measured but obtained from the 0.06 fm results, by applying the expected scaling factor α 2 s a 2 . The LECs L 7 and L 8 , both central values and errors, are taken from fit BE14 in Ref. [64].
where the functions g 1,a and g 2,a account for higher-order discretization effects, and the function h mπ includes analytical terms that parametrize higher order chiral effects. We have taken the pure counterterm contribution at two loops out of f cont 4 and written it separately. This contribution corresponds to the term proportional toC 4 , which is given by the combination of low energy constants (LECs) C 12 + C 34 − L 2 5 . The O(p 4 ) LEC L 5 can be extracted from global fits or from lattice-QCD calculations of light-light quantities, but the O(p 6 ) LECs C 12 and C 34 [59,60] are not known. (Only model-based estimates and imprecise global fit values exist.) We therefore takeC 4 as a constrained fit parameter. All dimensionful quantities entering in the fit function in Eq. (4.2) are converted into r 1 units by using the values of r 1 /a in Table III.
Since our simulations are performed in the isospin limit, m u = m d , f 2 and f 4 are evaluated for degenerate up and down quarks. The explicit NLO PQSChPT function f PQSChPT 2 (a) can be found in Ref. [57]. It incorporates the dominant discretization effects coming from the taste-symmetry breaking of staggered fermions. The function f PQSChPT 2 (a) depends on the HISQ taste splittings ∆ Ξ through the sea meson masses  2), as well as the posterior values obtained for those parameters in our preferred fit. The dimensionless χPT parameter s is given by the quantity 1/(16π 2 (r 1 f π ) 2 ) ≈ 0.3. The priors listed for the hairpin parameters are for the a ≈ 0.12 fm ensembles, and those for the other lattice spacings are obtained by rescaling these numbers, assuming that the hairpin parameters scale like the average of the ∆ Ξ . These values are obtained from fits to light-light quantities using two-loop PQChPT [65]. The uncertainty includes statistical and systematic errors. The prior central values for the NLO LECs are from fit BE14 in Ref. [64] with Λ χ = 0.77 GeV, while the prior widths are twice the errors in Ref. [64]. We fix the LECs L 7 and L 8 and give their values in Table III, as explained in the text. The entries "−0.000" denote small negative numbers that round to zero.

Fit parameters Gaussian priors
ChPT fit (central value ± width) posteriors 1.01 ± 0.12 1.00 ± 0.12 L r 6 (Λ χ ) × 10 3 0.14 ± 0.10 0.13 ± 0.09 with m i , m j sea quark masses, the slope µ to be determined by fits of the ChPT expressions to experimentally measured meson masses, and Ξ labelling the meson taste. Values of ∆ Ξ for each ensemble are given in Table III. The function f PQSChPT 2 (a) also depends on the tasteviolating hairpin parameters, δ V and δ A , which come from ChPT disconnected diagrams. We fix the taste splittings in the fit function to their values in Table III since they are precisely enough known that the corresponding errors do not affect our results significantly. The values are from Ref. [46], as well as unpublished updates with better statistics and the inclusion of new ensembles not previously analyzed. The uncertainty in the hairpin parameters is, however, quite large. We therefore treat them as constrained fit parameters with central values and widths equal to those in Table IV, determined from fits to light-light meson quantities [65]. Their uncertainty is thus propagated to the final fit errors.
For some of the meson masses that appear in f 2 there are no experimental measurements or lattice results, as for example, for m valence ss or for the sea-valence meson masses involving strange quarks. Because we use values given by NLO ChPT for these masses, our f 2 function has some dependence on the corresponding O(p 4 ) LECs L i . This is the best approximation we have, and we find that different implementations of higher-order corrections result in changes to the central values that are significantly smaller than the statistical errors.
The continuum NNLO ChPT function f cont 4 also depends on the O(p 4 ) LECs. We take most of them as constrained fit parameters with prior central values equal to the posteriors obtained in the O(p 6 ) global fit BE14 in Ref. [64]. We take as an input parameter the combination 2L 6 − L 4 instead of L 4 because the fit is more sensitive to that combination and because L 4 is fixed in fit BE14. The prior widths are set to twice the errors in Ref. [64]. The chiral scale, at which the LECs and chiral logarithms in the ChPT expression of Eq. (4.2) are evaluated, is set equal to the mass of the ρ meson, i.e., Λ χ = M ρ . The O(p 4 ) LECs from Ref. [64], used as priors here, agree within errors with (but are more precise than) the only realistic lattice calculations available at the moment: the N f = 2 + 1 MILC [13] and the N f = 2 + 1 + 1 HPQCD [12] calculations. The prior central values and widths used in our chiral-continuum fit to Eq. (4.2) are listed in Table IV.
The O(p 4 ) LECs L 7 and L 8 appear only in the isospin corrections and in the NLO expressions for some of the meson masses in f 2 and f 4 . Their effect on f 2 and f 4 in the isospin limit is, however, negligibly small, and their main impact is via the isospin-breaking corrections, which are added after performing the chiral-continuum fit (see Sec. V F). We chose then to take L 7 and L 8 as fixed parameters in the chiral-continuum fit and include their uncertainties in the total error as described in Sec. V B.
Once the O(α 2 s a 2 ) taste-violating discretization errors for staggered fermions are removed through the explicit dependence on a of f PQSChPT . We take these errors into account through the functions g 1,a and g 2,a in Eq. (4.2): where the K i are fit parameters,∆ = 1 16 (∆ P + 4∆ A + 6∆ T + 4∆ V + ∆ I ) is the average taste splitting, and r 2 1 a 2∆ is a proxy for α 2 s a 2 . Table IV lists the priors employed for the K i s. The terms proportional to K 2 and K 2 are generic terms parametrizing discretization effects of O(α s a 2 ) and O(α 2 s a 2 ), respectively, obeying the AG theorem. We include the terms proportional to K 1 and K 3 to account for O(α s a 2 ) and a 4 violations of the AG theorem at finite lattice spacing arising from symmetry-breaking discretization effects in the form factor decomposition, Eq. (1.4), and in the continuum dispersion relation. We find that adding an O(a 4 ) term instead of the one proportional to K 2 yields fit results that are nearly identical.
As in Refs. [24,25], we also add generic analytical terms corresponding to higher orders in the chiral expansion until the error of the chiral-continuum fit saturates. That happens at N 4 LO [O(p 8 )]-see Sec. V. The function h mπ in Eq. (4.2), which collects these effects, therefore takes the form h mπ =C 6 r 2 1 m 2 π +C 8 r 4 1 m 4 π . (4.5) The terms proportional toC 6 andC 8 are O(p 6 ) and O(p 8 ), respectively. TheC i are constrained fit parameters; the priors for them can be found in Table IV Table II and are corrected for the one-loop finite-volume effects also listed in that table. Different symbols and colors denote different lattice spacings. The data point at the a ≈ 0.042 fm ensemble includes the correction given in Sec. III B. The errors bars on the data points are statistical only, obtained from 500 bootstrap resamples. Data points at the same light-quark mass but different lattice spacings are offset horizontally. The gray continuous line shows the continuum extrapolation as a function of the light-quark mass, and the yellow star is the continuum result interpolated to the physical light-quark mass. The cyan error band, as well as the error bar on the physical point is the statistical chiral-continuum fit error (obtained using 500 bootstrap resamples), which includes discretization and higher order chiral errors, as well as the uncertainty from some of the input parameters, as discussed in the text. The continuum extrapolation line is obtained by setting the valence and sea light quark masses equal, setting m s to its physical value, turning off all discretization effects, and adding isospin-breaking effects.

A. Fit results
We fit our finite-volume corrected form factor data to the functional form in Eq. (4.2) with the functions g i,a and h mπ given in Eqs. (4.4) and (4.5), respectively. All the results listed in Table II are included in our central-value fit, except for the ensemble with a ≈ 0.12 fm and m π,P L = 3.2, which we use only to check finite-volume effects. We then extrapolate to the continuum limit and interpolate to the pure-QCD meson masses, i.e., with electromagnetic effects removed) using the parameters determined from the fit together with the continuum isospin-breaking NNLO ChPT expressions in Ref. [66]. We take the pure-QCD masses from Ref. [47]: 4 m QCD  4 The π 0 QCD mass is just the experimental one, and the π + QCD mass includes the estimate of the small isospin-breaking correction, which comes from Ref. [69]. and m QCD π + , and each sea meson mass to the relevant QCD isospin-breaking mass above, we find where the error is from the fit only, and does not yet include all systematic effects. The statistical error is estimated by fitting to a set of 500 bootstrap samples for each ensemble. The plot in Fig. 4 shows, as a function of the light quark mass, the central interpolation curve as well as its error band in the continuum and with the strange quark mass adjusted to its physical value. The result in Eq. (4.6) includes isospin corrections up to NNLO-see Sec. V F for more details. For K 0 decays, isospin corrections enter only at NLO (f 2 ) and beyond, and are small, < 0.15%. It also includes corrections for the leading-order finite-volume effects as described in Sec. III A.
The second column in Table IV shows the posteriors for the fit parameters of the chiralcontinuum fit that leads to the result in Eq. (4.6). We cannot determine the coefficients K i accurately since there is very little a 2 dependence in our results. In fact, if we remove the a ≈ 0.15 fm point we could fit our data without including discretization effects at all. Our lattice data also provide little constraint on the individual values of the O(p 4 ) LECs. As seen in Table IV, the posterior fit values of the L r i are generally the same as the priors.

V. SYSTEMATIC ERROR ANALYSIS
The error in Eq. (4.6) includes statistical, chiral-extrapolation, and discretization errors, as well as the uncertainties associated with the inputs that are treated as constrained fit parameters: O(p 4 ) LECs (except L 7,8 ) and taste-violating hairpin parameters. The uncertainties of the data and constrained input parameters are propagated through the fit via 500 bootstrap resamples.
In this section, we further study these sources of uncertainty, perform tests of the stability of our preferred fit strategies, and estimate the other sources of systematic error entering in our calculation of f Kπ + (0): uncertainty in the inputs, scale error, partial-quenching effects, higher-order finite volume effects, isospin-breaking corrections, and the effects of non-equilibrated topological charge.

A. Fit function, discretization error and chiral interpolation
Because we have data at the physical light quark masses, the chiral fit is an interpolation, and is largely independent of the precise form of the fit function and the values of the ChPT parameters. We have performed a number of tests to check this stability under variations in the fit function and to estimate the effect of higher-order terms in the chiral and Symanzik expansions. Figure 5 shows the tests performed and discussed in this subsection, together with the results from additional fits discussed in the following subsections, which we use to estimate several systematic uncertainties.
First of all, in order to be sure that effects due to higher-order terms in the chiral expansion are properly included in the Bayesian analysis that leads to the fit error shown in Eq. (4.6), we need to check that this error stabilizes as we add higher-order chiral terms. The point labelled NNLO in Fig. 5 includes only terms up to NNLO, i.e., without the h mπ function in Eq. (4.2). Only minimal changes in the central value and errors are produced by the addition of the N 3 LO termC 6 (m 2 K − m 2 π ) 2 m 2 π from Eqs. (4.5) and (4.2). The difference between this and the base fit is the N 4 LO termC 8 (m 2 K − m 2 π ) 2 m 4 π , and we see that it has a negligibly small effect.
Different values of the decay constant used in the two-loop (NNLO) term f cont 4 are equivalent up to omitted higher-order terms, and therefore should have a negligible effect on our analysis. In contrast, the decay constant at one loop has to be set equal to the physical pion decay constant, f π , in order to be consistent with the particular expression chosen for f cont 4 .
In the base fit, we use f π as the chiral expansion parameter in f cont 4 . We check that other possible choices, such as f K = (155.6 ± 0.4) MeV [16] (labeled "f K vs. f π at two loops" in Fig. 5) and an estimate of the decay constant in the chiral limit, f 0 = (113.5±8.5) MeV [13], shift the central value in Eq. (4.6) by less than 0.06%, well under the statistical error.
As another test of the ChPT fit and errors, we replace the continuum two-loop ChPT expression in Eq. (4.2), f cont 4 , by an analytic function, and consecutively add N 3 LO and N 4 LO analytic terms-see results labeled "NNLO analyt.," "N 3 LO analyt.," and "N 4 LO analyt.," respectively, in Fig. 5. All the results agree very well with our base fit within statistics, being nearly identical once the N 3 LO analytic term is included.
The three results labeled "no analyt. a 2 " (which corresponds to a fit without including g 1,a and g 2,a in the fit function), "α 2 s a 2 (m 2 π − m 2 K )", and "α 2 s a 2 (m 2 π − m 2 K ) + α s a 2 " in Fig. 5 represent a check that the discretization errors are properly included in the fit error of Eq. (4.6). Once we include the term of order α 2 s a 2 (m 2 K − m 2 π ) 2 (proportional to K 2 ) in Eq. (4.4), which is required to get a good fit, the central value and errors barely change with the addition of α s a 2 corrections (the result labeled "α 2 s a 2 (m 2 π − m 2 K ) + α s a 2 "). Adding the two remaining discretization terms in g 1,a and g 2,a , which returns us to the base fit, makes no noticeable difference. The rapid stabilization of the fit reflects the tiny lattice-spacing dependence of our data.
Of all data in Fig. 2, only the point at a ≈ 0.15 fm shows what appear to be significant discretization effects. Dropping that data point has the effect of increasing the errors (see result labeled "no a ≈ 0.15 fm"), since the other ensembles provide very little constraint of the analytical a 2 fit parameters. In fact, after dropping that point, we can fit our remaining data with a continuum fit function, although we see from Fig. 5 (result labeled "continuum, no a ≈ 0.15 fm") that the result is larger than our central result by about two standard deviations, measured in terms of the fit errors. Adding analytical discretization corrections via the functions g 1,a and g 2,a to the continuum fit function allows us to fit all our data, giving a result that is consistent with the base fit (see result labeled "continuum + analyt. a 2 "), although with a larger error.
In contrast to the noticeable effect of the coarsest ensemble on the total error, the effect of our finest lattice spacing, a ≈ 0.042 fm, on the central value and the error is very small since statistics in this ensemble is limited and, in addition, it has m l = 0.2m s , and thus it is relatively far from the physical point.
As shown in Fig. 5, both the ensembles with physical quark masses and those with unphysical masses are important in fixing the central value and reducing the fit error. The larger error of the fit including only physical-quark-mass ensembles reflects the weaker constraints on the higher-order discretization terms and the lack of constraints on the higher-order chiral terms, which can have an effect on the results from nominally "physical" ensembles due to mistunings of the strange and light quark masses. On the other hand, the larger error of the fit including only the unphysical-quark-mass ensembles reflects primarily the error of the chiral extrapolation.
For the reasons discussed above, the statistical fit error shown in Eq. (4.6), which is obtained with our base fit using Eq. (4.2) together with the higher-order chiral and discretization terms in Eqs. (4.4) and (4.5), properly includes the errors from higher-order discretization effects and chiral corrections in addition to the statistical errors. The inclusion of the unphysical light-quark-mass data in our ChPT description gives us a handle on these higher-order effects and allows us to robustly correct for mass mistunings and estimate the error associated with the truncation of the corresponding series.

B. Inputs for the fixed parameters in the chiral function
The values and errors of the fixed inputs we use in our chiral-continuum interpolation/extrapolation are listed in Table III. The HISQ taste splittings are known precisely enough that their errors have no impact on the final uncertainty. Similarly, when we change the pion decay constant within its error and repeat the fit, results for the form factor are unchanged at the precision we quote. The uncertainty is small because the dependence on f π enters through the coefficients and parameters in the ChPT fit function, which, as discussed above, already have little effect on the results. Finally, by varying Λ χ in the range M ρ ± 0.5 GeV, we have checked that our results are independent of the chiral scale, as they should be. We therefore do not need to add an uncertainty due to the errors in the inputs or the choice of chiral scale to the statistical fit error.
However, the LECs L 7 and L 8 (which we treat as fixed input parameters, unlike the other LECs), do have an impact on the form factor error, mainly through their effect on the isospin corrections. We estimate this uncertainty by varying their central values by their respective standard deviations, repeating the fit, and recalculating the form factor (including isospin corrections). We take the shift that this variation produces as the uncertainty associated with these LECs, and add it in quadrature to the fit error, as shown in Table V. The above procedure does not underestimate the error due to these LECs, since if we treat L 7 and L 8 as constrained fit parameters instead, the same as the other O(p 4 ) LECs, we obtain a slightly smaller total error.

C. Lattice scale
We rewrite all the dimensionful quantities entering in the two-loop ChPT fit function in r 1 units, where the r 1 scale is obtained from the static-quark potential [67,68]. The lattice parameters are converted to r 1 units by multiplying by the values of the relative scales r 1 /a in Table III, while the physical parameters are converted by using r 1 = 0.3117 (22) fm [62].
The form factor f Kπ + (0) is a dimensionless quantity, and thus the effect of the error in the lattice scale is small. When we change the scale r 1 by its error, the central value only shifts by ±0.0008. We include this variation as a systematic error in Table V. The errors in the relative scales r 1 /a, on the other hand, have no significant impact on our results.

D. Partial-quenching effects in m s at NNLO
The valence and sea strange quark masses differ on some of the ensembles as explained in Sec. II, leading to partial-quenching effects-see Table I. These effects can be exactly treated at NLO within the PQSChPT framework, but at NNLO only the full-QCD ChPT expressions are available. We then have the choice of using either the sea or the valence m s at NNLO and beyond. In practice, this ambiguity only affects f cont 4 in Eq. (4.2) since the factor (m 2 π − m 2 K ) 2 in that equation comes from the valence sector. The result in Eq. (4.6) is obtained using the valence strange quark masses at NNLO. If we use the sea strange quark masses at NNLO instead, f K 0 π − + (0) shifts by 0.013%, which we include on the line labeled "m val s = m sea s " in the error budget. This systematic effect is small because the sea strange quark masses are generally well tuned on the HISQ N f = 2 + 1 + 1 MILC ensembles, and m val s = m sea s on the most relevant ensembles in the chiral-continuum interpolation/extrapolation, the ensembles with physical quark masses and a ≈ 0.09, 0.06 fm.

E. Higher-order finite-volume corrections
In our previous calculation [24,25], the uncertainty due to finite-volume effects was one of the two dominant sources of error. (The other was the fit error.) The finite-volume error was estimated to be of the same order as the statistical error from a comparison of the lattice data from two different volumes, with other parameters held fixed. Then, although very small, 0.2%, this error turned out to be a limiting factor for precision. In this work we have increased the statistics on the ensembles analyzed in Ref. [24,25] to check finite size effects. We have also sharpened this direct comparison by generating data on a third, smaller, volume. The three ensembles are those with a ≈ 0.12 fm and m l /m sea s = 0.1 in Table I and Fig. 1. Table II gives the values for f Kπ + (0) on these three volumes. The results on the two largest volumes are essentially the same, while that on the smallest volume differs from the others by less than the statistical error. From this comparison alone we could conclude that finite-volume effects are smaller than 0.17%, the smallest statistical error on the three ensembles.
To reduce the error further we use NLO staggered partially-twisted partially-quenched ChPT [41] to correct the form factor prior to the chiral-continuum fit, as described in Sect. III A. The resulting finite-volume corrections are ≤ 0.1% on all ensembles. If we did not correct our data for finite-volume effects at one loop, the result for f K 0 π − + (0) would shift by 0.00051. We use a power-counting estimate for the error due to neglected higherorder finite-volume effects, by taking the calculated one-loop correction and multiplying it by a typical chiral-loop suppression factor. For quantities involving a strange quark, we may estimate this factor to be m 2 K /(8π 2 f 2 π ) ≈ 0.18. The size of the ratio of NNLO and NLO contributions to f K 0 π − + (0) that we obtain in this work is a bit larger, ≈ 0.26. We conservatively multiply the calculated finite-volume correction by the larger number, which yields a 0.014% uncertainty that we include in the error budget in Table V.

F. Isospin-breaking corrections
Isospin-breaking corrections accounting for the difference between the up and down-quark masses can be calculated in the ChPT framework and thus written as a chiral expansion starting at NLO for neutral kaons S,K 0 π − . . . , (5.1) where the parameters ζ  [69] and [66], respectively. These corrections depend on the lowest-order π 0 −η mixing angle ε (2) , or, alternatively, the quantity R ≡ (m s −m)/(m d −m u ) withm ≡ (m u + m d )/2. In order to arrive at the number in Eq. (4.6), we use the expressions in Ref. [66], the QCD meson masses quoted in Sec. IV A, and the values of the LECs obtained from our fits and shown in Table IV (for L r 7 and L r 8 we take the input values in Table III). The only combination of O(p 6 ) LECs that enters at this order in the isospin-breaking terms for K 0 → π − ν decays is C 12 + C 36 . This combination, which we obtain from our fitting procedure, is the same one that appears in the isospin limit.
As with higher-order finite-volume corrections, we use a power-counting estimate for the error due to isospin corrections not included in our result, N 3 LO and higher. Taking the calculated NNLO isospin-breaking correction, −0.00057, and multiplying it by the chiral loop suppression factor discussed in Sec. III A gives an error of ±0.015%.
Another source of error is the parametric uncertainty in the isospin-breaking quantity R The analysis that yields to this result is the same as in Ref. [47], except that we have included more configurations at the ensembles with a ≈ 0.06 fm and a ≈ 0.042 fm, and included the a ≈ 0.15 fm data in the central fit. The electromagnetic errors are estimated as in Ref. [72]. We estimate the error on the form factor coming from the uncertainty on R by varying this quantity within its error and redoing the fit. As expected, the impact on the form factor for the neutral mode is nearly negligible, 0.002%. Nevertheless, we include it in our error budget.

G. Nonequilibrated topological charge
As described in Sec. III B, a correction due to improper sampling of the topological charge is needed only on the a ≈ 0.042 fm ensemble with m l = 0.2m s , where we obtain ∆ Q f Kπ + (0) = 0.00018. Not surprisingly, given that (1) this ensemble has little influence on the chiral-continuum interpolation/extrapolation (see Fig. 5 for the effect of removing the ensemble completely), and (2) the correction is much smaller than the statistical error on the ensemble (see Table II), the effect of the correction on the physical value of f Kπ + (0) is negligible. We therefore do not add an uncertainty due to this effect to our error budget.

Our final result for the vector form factor is
where the first error in the middle expression is the combined statistical, discretization and chiral interpolation error discussed in Sec. IV A, and the second the sum in quadrature of all the systematic errors discussed in Sec. V. Table V summarizes all sources of error in our calculation. The total uncertainty is the smallest achieved to date. We compare our result for f K 0 π − + (0) with the results from the most recent lattice calculations and phenomenological approaches in Table VI, and with the results entering the FLAG  TABLE VI: Form factor f K 0 π − + (0) as extracted from the most recent lattice calculations (first half of the table), from phenomenological approaches using of two-loop ChPT, and from the 1984 calculation by Leutwyler and Roos, which uses one-loop ChPT and a quark-model for higher order terms. For those calculations based on two-loop ChPT, we also indicate the method used in the estimate of the O(p 6 ) LECs.
Group (15)(11) staggered fermions (N f = 2 + 1 + 1) ETM [27] 0.9709(45)(9) twisted-mass fermions (N f = 2 + 1 + 1) Fermilab Lattice/MILC [24] 0.9704(24) (22) staggered fermions (N f = 2 + 1 + 1) JLQCD [31] 0.9636(36)( +49 −19 )(29) overlap fermions (N f = 2 + 1) RBC/UKQCD [26] 0.9685 (34)(14) domain-wall fermions (N f = 2 + 1) Bijnens & Ecker [64,73] 0.970 (8) ChPT + NNLO global fit Kastner & Neufeld [74] 0.986 (8) ChPT + large N c + dispersive Cirigliano et al. [75] 0.984 (12) ChPT + large N c Jamin, Oller, & Pich [76] 0.974 (11) ChPT + dispersive (scalar form factor) Bijnens & Talavera [58] 0.976 (10) ChPT + Leutwyler & Roos Leutwyler & Roos [77] 0.961 (8) One-loop ChPT + quark model average and those from phenomenological approaches in Fig. 6. Our value for f K 0 π − + (0) agrees within errors with previous N f = 2 + 1 and N f = 2 + 1 + 1 lattice calculations. In particular, the value is close to the other N f = 2 + 1 + 1 results, but with significantly smaller errors. It also agrees with the most recent phenomenological determinations [64,73], which are based on two loop ChPT with LECs determined by NNLO global fits. The lattice results in Table VI We can thus use the values ofC 4 and L r 5 from our fit output in Table IV to extract the combination of O(p 6 ) LECs involved. Taking correlations into account, we find 3) The first error in Eq. (6.3) includes statistics, chiral extrapolation and discretization errors, as well as the uncertainty from the LECs (except L 7 and L 8 ) and the taste-violating hairpin parameters, as discussed in Sec. V. The second error is the sum in quadrature of the rest of the systematic uncertainties. The detailed error budget is in Table VII. We obtain all the errors in the same way as for f Kπ + (0). Isospin corrections do not apply to this quantity  6: Comparison of f K 0 π − + (0) from this analysis with previous lattice results entering in the FLAG averages [6] together with those averages for N f = 2 + 1 + 1 and N f = 2 + 1, as well as nonlattice determinations based on ChPT. The beige band corresponds to our result. The references and numerical results for all determinations are given in Table VI. since it is defined in the isospin limit. In practice, the values of LECs coming from a fit may be significantly affected by the presence or absence of higher-order chiral terms in the fit function. Therefore, applications of our result in Eq. (6.3) should allow the same type of corrections as in (the continuum limit of) Eq. (4.2). The complete error budget for this quantity can be found in Table VII. Our result in Eq. (6.3) agrees with non-lattice determinations in Refs. [75][76][77]. In those papers, the contribution to f + (0) from C 12 + C 34 was calculated using the large N c approximation, a coupled-channel dispersion relation analysis, and a quark model, respectively. However, the value for C 12 + C 34 − L 2 5 found in Ref. [74], which is based on ChPT, large N c estimates of the LECs, and dispersive methods, is ∼ 3σ smaller than our value, [C r 12 + C r 34 − (L r 5 ) 2 ] (M ρ ) = (2.92 ± 0.30) × 10 −6 . The result in Eq. (6.3) also agrees very well with our previous calculation of this combination of LECs in Ref. [28], on the MILC N f = 2 + 1 asqtad configurations, although with greatly reduced errors. In fact, all sources of error are reduced due to several factors: the use of the MILC N f = 2 + 1 + 1 HISQ configurations with smaller discretization errors than the asqtad action, data at smaller lattice spacings, data with physical light-quark masses, better tuning of the strange sea quark masses, and including NLO finite volume corrections explicitly. The agreement with the JLQCD result in Ref. [31] is borderline, but the JLQCD calculation relies on simulations at a single lattice spacing, although a systematic error is quoted for it, and it does not include data at the physical light-quark masses. Those systematics could affect more strongly the value of the combination of LECs than the form factor itself. where the first error is from the uncertainty on the form factor, and the second is the experimental uncertainty. Both errors are now of the same size. The experimental error in Eq. (7.1) includes the uncertainty on the long-distance electromagnetic and strong isospin-breaking corrections, δ Kl EM and δ Kπ SU (2) , which are taken into account when doing the experimental average of the neutral and charged modes [7]. This uncertainty is however dominated by the errors in the lifetime and branching-ratio measurements of the neutral-kaon modes [7]. Other uncertainties such of those from the phase-space integrals are insignificant [7].

VII. PHENOMENOLOGICAL IMPLICATIONS
In Fig. 7 we compare our extraction of |V us | from K semileptonic decays with other determinations using K semileptonic and leptonic decays, and hadronic τ decays. Our semileptonic determination of |V us | is the most precise to date not relying on an external input for |V ud |. The central value agrees very well with the most recent lattice and nonlattice semileptonic calculations, as well as with those based on hadronic tau decay, the latter have much larger errors. Our result, however, is in tension with the leptonic determination using f K /f π and with the unitarity prediction given by |V us | = 1 − |V ud | 2 with |V ud | from Ref. [2]. The agreement with the leptonic determination using f K is borderline. The sizes of the disagreements-2.6σ with unitarity and 2.2σ with the leptonic determination using f K /f π -are similar to those using other recent lattice calculations for the semileptonic vector form factor.
As a consistency check of the semileptonic extraction of |V us |, we can consider the neutral-and charged-kaon modes separately. Using our result in Eq. (6.1) together  [24,26,27], respectively. The leptonic determinations, labeled K l2 , use as inputs the 2 + 1-flavor lattice-QCD average f K from FLAG [6], which only includes calculations where the lattice scale is set from physical inputs other than f π , and the most recent and precise determination of f K ± /f π ± from Ref. [47]. The inclusive hadronic τ -decay determinations are the most recent ones, from Boyle et al. 2018 [38] and Hudspith et al. 2017 [37]. The second value from Ref. [38] comes from relating the τ → K ν branching fraction to the K µ2 branching fraction to get the experimental contribution from the K pole. The two values in Ref. [37] correspond to using the normalization for τ decays into Kπ modes as obtained in Ref. [36] or as given by HFLAV [40]. For the exclusive τ determination we follow the calculation by the HFLAV group [40], but we update the value of the ratio f K ± /f π ± to that in Ref. [47]. The unitarity value is taken to be |V us | = 1 − |V ud | 2 with |V ud | from Ref. [2]. RC stands for radiative corrections. The dotted magenta vertical lines correspond to this unitarity value. The gray vertical band corresponds to our result in Eq. (7.1).
with the experimental average for neutral modes only [7], |V us |f K 0 π − + (0) = 0.2163(5), 5 we can compare |V us | as extracted exclusively from neutral kaon decays: |V us | K 0 π − = 0.22309 (43) (66). In this case, we can disentangle the purely experimental error from the uncertainty in the long-distance electromagnetic corrections, δ Kl EM , which is the same for all neutral modes. This result is in very good agreement with the value in Eq. (7.1) within errors, which constitutes a good test of the ChPT calculation of isospin (larger for the charged modes) and EM (larger for the neutral modes) corrections included in the total experimental average.

B. Tests of CKM unitarity
Using our main result for |V us | in Eq. (7.1), the value |V ud | = 0.97420(21) from superallowed nuclear β decays [2], and noting that |V ub | 2 is negligible, we find that the measure of deviation from first-row CKM unitarity in Eq. (1.1) is which is ∼ 2.1σ away from the unitarity prediction, with an error dominated by the uncertainty on |V ud |. This makes revisiting the determination of |V ud | a priority for CKM tests.
In this vein, one should examine not only superallowed β decays but also other approaches. At present, the precision in the extraction of |V ud | from the measurement of the neutron lifetime [3] or pion β decays [4] is still far from that obtained from superallowed β decays. In the case of superallowed β decays, additional measurements will have a small effect on |V ud |. At the moment, the greatest improvement would come from a calculation of the short-distance radiative correction, which is the main source of uncertainty [2]. A very recent calculation of the nucleus-independent contribution to those corrections, following a new methodology based on dispersion relations [17], obtains a value around 2σ larger than the current best determination by Marciano and Sirlin [78] and with a significant reduction of the error. The increased electroweak radiative correction, when combined with the superallowed β decay results [2], results in a lower value of |V ud |. The authors of Ref. [17] quote |V ud | = 0.97366 (15). Together with our result for |V us |, this value of |V ud | considerably increases the tension with unitarity: 3) a more than 5σ discrepancy. We discuss further phenomenological implications of this new calculation in Sec. VII D. For the remainder of this section, we use the result by Marciano and Sirlin [78], which leads to |V ud | = 0.97420 (21) and Eq. (7.2). To avoid using |V ud | as an input, we can instead perform a unitarity test relying only on experimental kaon-decay measurements [7], on the lattice input from the most recent determination of f K + /f π + [47], and on our result in Eq. (6.1) for f K 0 π − + (0). The result of the unitarity test using those inputs, noting again that |V ub | is negligible, is 6 where the 2.2σ deviation from unitarity is a reflection of the tension between the leptonic and semileptonic determinations of CKM matrix elements. These results are shown in Fig. 8, together with the test that takes |V ud | from superallowed β decays as input. 6 The disentanglement of the EM and experimental errors in Eq. (7.4) is approximate, and intended only to indicate the relative size of these errors. The separation of the sources of error is precise for leptonic decays, but for semileptonic decays we assume an overall 0.11% EM error in the uncertainty of the experimental average. This should be a fairly good approximation, however, since the average is dominated by the neutral modes for which the error is indeed 0.11%. One can perform another test of the unitarity of the CKM matrix by comparing |V us | with |V cd |, which in the SM should be equal up to corrections of O(λ 5 ), with |V us | = λ + O(λ 7 ). Including the O(λ 5 ) corrections, which only affect the last significant digit, the most precise determination of |V cd | = 0.2151(6) f D (49) expt (6) EM from leptonic decays [47] implies the value |V us | |V cd | = 0.2158 (52). This value of |V us | agrees at the 1.4σ level with our result in Eq. (7.1), although with an uncertainty that is an order of magnitude larger. The uncertainty is dominated by the experimental error on the leptonic decay rate D + → + ν, which is expected to be reduced by BESIII and Belle II. This result is also depicted in Fig. 8. As is the case for our main result, |V us | |V cd | is in tension with first-row CKM unitarity by about 2σ when it is used together with |V ud | from superallowed nuclear β decays in Eq. (1.1).
Note that, in order to perform this test, we change the normalization of the decay constant f D + obtained in Ref. [47] to account for a change in the scale-setting quantity in that work, f π + , from the PDG value f π + = 130.50 ± 0.13 MeV [3] to the FLAG average f π + = 130.2 ± 0.8 MeV [6]. That gives us f D + = 212.
The reason for that change is that the PDG value relies on an external input for |V ud |, which is taken from superallowed nuclear β decays, which obscures the comparison. The FLAG number, however, is an average of direct lattice determinations of f π + . With this choice of f π + , the errors are fairly large, and the value of |V ud | extracted from experimental data on pion leptonic decays [16] agrees within ∼ 1.5σ with both |V ud | from superallowed nuclear β decays and the value from kaon decays only that we discussed above.  9: Comparison of the unitarity point using |V ud | = 0.97366 (15) with the results in this work, and with the unitarity point corresponding to |V ud | = 0.97420 (21). RC stands for radiative corrections.
The unitarity test comparing |V us |/|V ud | with |V us | |V cd | /|V ud |, again including corrections up to O(λ 5 ), and taking the decay constants f K + /f π + and f D + /f π + from Ref. [47] and the experimental data on leptonic experimental data from Ref. [16], fails at the 2σ level. This test is limited by the experimental error on the D + leptonic decay rate.
D. Implications of the new extraction of |V ud | If the decrease of the central value and uncertainty of the nucleus-independent electroweak radiative corrections involved in the extraction of |V ud | from superallowed β decays in Ref. [17] is confirmed, the new value |V ud | = 0.97366 (15) would exacerbate some of the tensions we have just discussed.
First, as shown above, this value of |V ud | and our semileptonic result for |V us | would imply a greater than 5σ violation of first-row CKM unitarity. The tension between our semileptonic value of |V us | and the one extracted from kaon leptonic decays and f K ± /f π ± , however, would be slightly reduced to 2σ, since a smaller value of |V ud | would give a smaller value of the leptonic |V us |, closer to our semileptonic extraction. For the same reason, the tension between the ratios involving |V ud | in Eqs. (7.6) and (7.7) would be slightly lessened.
In Fig. 9, as an example, we show the comparison of the unitarity prediction 1 − |V ud | 2 for |V us | using both |V ud | = 0.97366 (15) and |V ud | = 0.97420 (21), together with the results in this work. Given the important implications of a value of |V ud | with a smaller error and a smaller central value, it is very important to confirm the new calculation of radiative corrections in Ref. [17], and to understand the discrepancy with the previous best determination in Ref. [78].

VIII. CONCLUSIONS AND OUTLOOK
Using the HISQ N f = 2 + 1 + 1 MILC ensembles, we have performed the most precise computation to date of the vector form factor at zero momentum transfer, f K 0 π − + (0), and the first one to include the dominant FV effects, as calculated in ChPT at NLO. Our result for the form factor enables a direct determination of the CKM matrix element |V us | from semileptonic kaon decays with a theory error that is, for the first time, at the same level as the experimental error. Further, the uncertainty in this direct determination is now similar to those from indirect determinations based on leptonic decays with |V ud | as input.
A key to achieving this level of precision is simulating at near-physical values of the quark masses, which drastically reduces the systematic errors associated with the chiral extrapolation (replacing it with an interpolation), as well as the error coming from the chiral LECs that are inputs to the analysis. The finite-volume effects, one of the main sources of uncertainty in our previous analyses, have also been significantly reduced by explicitly including them at NLO (the leading nontrivial order) in ChPT. The dominant remaining source of error is now statistics, which could be reduced by extending the key ensembles with physical quark masses, and including the existing MILC physical-mass ensemble with a finer lattice spacing of a ≈ 0.042 fm.
Another important error arises from the uncertainty in the ChPT LECs of order p 6 . That uncertainty could be reduced by performing a combined analysis of form-factor data together with light meson masses and decay constants, which would put more constraints on the ChPT LECs. In particular, the error from L 8 is comparable to, but greater than, that from L 7 , and the combined analysis could significantly reduce the L 8 error. Errors from L 4 , L 5 , and L 6 would also be reduced, but they have a much smaller effect on the total error here.
We find that the extraction of |V us | from semileptonic kaon decays is in tension both with the extraction from leptonic kaon decays and with unitarity at the ∼ 2-2.6σ level. In particular, the unitarity test based only on kaon decay data, without any external input for |V ud |, and having as nonperturbative inputs f K 0 π − + (0) from this work and f K ± /f π ± from Ref. [47], shows a ∼ 2.2σ tension. While unitarity tests based on |V ud | are currently limited by the uncertainty in that matrix element, the tension with unitarity would raise to the 5σ level if the new calculation of radiative corrections involved in the extraction of |V ud | from superallowed β decays [17] is confirmed.
The test based on kaon-decay data has similarly-sized uncertainties arising from both theory and experiment. In order to shed light over these tensions, improvements from both the theoretical and experimental sides are urgently needed, as are improvements in other approaches. A new round of experiments are expected to reduce the experimental error to ∼ 0.12% in the next few years [7]. More importantly, the new high-statistics data will help to check the consistency of current fits, and to perform a more thorough study of systematic errors on the experimental averages.
For the experimental determination of |V us |f K 0 π − + , electromagnetic and isospin effects are currently being estimated using phenomenology and ChPT techniques. Although they are not yet a dominant source of error (EM effects make a 0.11% correction to the individual neutral channels), with the reduction of other sources of error and the forthcoming improvement in the experimentally measured branching ratios and lifetimes, they will eventually need to be included directly in the lattice-QCD simulations. Recent efforts in that direction can be found in Refs. [79][80][81][82].
Isospin corrections are numerically important for the charged kaon channels, where those effects enter already at LO through π 0 -η mixing. The NNLO ChPT estimate of the corrections for the charged modes has large errors [66] due to the unknown value of the O(p 6 ) LECs. Fortunately, the experimental average is dominated by the neutral kaon channels, so the charged-mode uncertainty does not have a large effect on the final experimental average. The strong isospin-breaking correction δ SU(2) used in the experimental average is a NLO ChPT estimate that partially includes NNLO effects; it does not include the uncertainty associated with higher-order terms in the chiral expansion. However, the fact that the value used in the average and the one extracted from experiment are so close (2.45(19)% vs. 2.82(38)% [7]), that the result for |V us | using only the neutral modes agrees with the one using all decay modes-see Sec. VII A, and that neutral modes are the dominant ones in the average, indicate that the experimental average using this estimate is robust.
The uncertainties from the phase-space integrals are insignificant at present in the final error for the experimental average. It is therefore not crucial at present to have a better representation of those, i.e., to have the q 2 dependence of the form factors. In the future, however, lattice calculations of f Kπ + (q 2 ) could provide better determinations of the form factor slope than those relying on experimental data [7,27].
An important future step in the investigation of the tensions observed in the first-row unitarity relation, and in the value of |V us | extracted from different sources, will be to perform a correlated analysis of semileptonic and leptonic kaon decays. That analysis would provide a more precise value of the ratio [f K ± /f π ± ] / |V ud |f Kπ + (0) and potentially give an insight into the tensions. Another key point in the study of those tensions is clarifying the role of the electroweak radiative corrections in the extraction of |V ud | from superallowed β decays, as well as reducing the error of that CKM matrix element as extracted, not only from superallowed β decays, but from other sources.