The kaon semileptonic form factor in Nf=2+1 domain wall lattice QCD with physical light quark masses

We present the first calculation of the kaon semileptonic form factor with sea and valence quark masses tuned to their physical values in the continuum limit of 2+1 flavour domain wall lattice QCD. We analyse a comprehensive set of simulations at the phenomenologically convenient point of zero momentum transfer in large physical volumes and for two different values of the lattice spacing. Our prediction for the form factor is f+(0)=0.9685(34)(14) where the first error is statistical and the second error systematic. This result can be combined with experimental measurements of K->pi decays for a determination of the CKM-matrix element for which we predict |Vus|=0.2233(5)(9) where the first error is from experiment and the second error from the lattice computation.


Introduction
In the Standard Model (SM) the unitary Cabibbo-Kobayashi-Maskawa (CKM) matrix [1,2] parameterises flavour changing processes and provides the only source of CP-violation in the quark flavour sector. Determining its elements from a combination of theory predictions for flavour-changing SM processes and the corresponding experimental measurements allows for testing its unitarity and consistency. Any deviation from the expected behaviour would signpost new physics.
In this work we present a new lattice QCD computation of the hadronic contribution to the semileptonic K → π decay, the vector form factor f Kπ + (0), using a lattice fermion formulation which respects the chiral symmetry of continuum QCD. When combined with the SM analysis of the experimental results for this kaon decay channel, summarised as |V us f Kπ + (0)| = 0.2163(5) [3] (for K 0 → π − ), it determines the matrix element |V us |. In combination with the other first-row elements of the CKM matrix, |V ub | and |V ud |, this allows an immediate test for CKM unitarity in the SM and constrains models for extensions of the SM.
The history of recent efforts in predicting f Kπ + (0) [4][5][6][7][8][9] is nicely summarised in [10,11]. This quantity is currently known with an overall uncertainty of 0.3% and improving this precision further is mandatory in view of experimental progress [12]. The error budget is typically dominated by the statistical uncertainty resulting from the Monte Carlo sampling of the QCD path integral in lattice QCD. Until recently the largest systematic uncertainty arose from the fact that simulations were only feasible for QCD with unphysically heavy pions. In order to make predictions for QCD as found in nature the simulation data had to be extrapolated using effective field theory or phenomenological models. Advances in algorithmic methods and computing technology now allow us to carry out simulations directly at the physical pion mass, thereby removing the dominant systematic uncertainty from the extrapolation.
In this work we present the first prediction of the form factor f Kπ + (0) with physical valence and sea quark masses in the continuum limit of domain wall lattice QCD with N f = 2 + 1 dynamical flavours. The physics described by our simulations corresponds to nature up to isospin breaking in the light quark masses and electromagnetic effects and the contribution arising from charm (and heavier quark) vacuum polarisation effects.
We anticipate the final results presented in this paper: f Kπ + (0) = 0.9685 (34) (14) , |V us | = 0.2233(5)(9) , (1.1) for the K → π form factor at vanishing momentum transfer and the CKM-matrix element for u → s flavour changing processes, respectively. The errors are statistical and systematic, respectively. The rest of this paper is organised as follows: In Section 2 we explain the computational strategy for determining the form factor from Euclidean two-and three-point functions. In sections 3 and 4 we discuss the simulation parameters and some aspects of the computational setup and the data analysis. Section 5 details how we predict the physical results from a very small interpolation of the simulation data to the precise physical quark mass combined with an extrapolation to the continuum limit. A discussion of residual systematic errors follows in section 6 and our final results and conclusions are presented in section 7.

Strategy
In this section we define the observables from which we can determine the K → π vector form factor f Kπ + (q 2 ), where q = p K − p π is the momentum transfer between the kaon and pion. The form factor is defined in terms of the QCD matrix element π(p π )|V µ |K(p K ) = f Kπ of the flavour changing vector current V µ = Z Vū γ µ s, where u and s are up-and strangequark fields and Z V is the vector current renormalisation constant. As first noted in [13] in the context of charm semileptonic decays an alternative way to determine the form factor is to consider the matrix element of the flavour changing scalar current S =ūs. From the vector Ward-Takahashi identity we derive which determines the vector form factor at zero momentum transfer by means of the identity f Kπ + (0) = f Kπ 0 (0). Note that in the above equation the renormalisation constants of the scalar current and the quark masses cancel.
In practice we determine the two matrix elements (2.1) and (2.2) from the time dependence of combinations of Euclidean QCD two-and three-point correlation functions which are the output of the actual simulation. The two-point function is defined as where i = π or K, and O s,i are pseudoscalar interpolating operators for the corresponding mesons, O s,π =ū ω s γ 5 d and O s,K =d ω s γ 5 s. The subscript s indicates the smearing type induced via the spacial smearing kernel ω s which in the simulations presented here is local (s = L, ω L (x, y) = δ x,y ) or gauge-fixed wall (s = W , ω W (x, y) = 1). In practice we employ s 1 = W and s 2 = L, W . The constants Z s,i are defined by Z s,i = P i | O † s,i | 0 where P i is a pion or a kaon. As a result of our choice of boundary conditions C i (t, p i ) is symmetric with respect to the middle of the temporal direction of length T and we have assumed that t is such that the correlation function is dominated by the lightest state (i.e. the groundstate pion or kaon) with energy E i (p i ). As explained further down meson momenta are induced using partially twisted boundary conditions. The three-point functions are defined as (exposed here for the case t i = 0) where Γ ∈ {V µ , S} is the vector/scalar current with flavour quantum numbers that allow for the P i → P f transition and Z s,f = 0 | O s,f | P f where again P f is either a pion or a kaon. We will refer to t i (t f ) as the source (sink) time plane, respectively, and t is the time plane where the vector/scalar current is inserted. The constant c Γ satisfies c V 0 = −1 for the time-component of the vector current and c V i ,S = +1 for the spatial components of the vector current and for the scalar current. In the second line of the above expression we have assumed that all time intervals are sufficiently large for the lightest hadrons to give the dominant contribution. Expressions (2.3) and (2.4) apply in slightly modified form in our earlier simulations [5][6][7] due to differences in technical details (noise-and sequential source propagators [14]). We consider two strategies for determining the form factor f Kπ + (0) from the above correlation functions. The first consists of a simultaneous fit over all two-and threepoint functions from which the matrix elements (2.1) and (2.2) are readily extracted. The other strategy is to determine both matrix elements from fits to suitably chosen ratios of correlation functions [15], e.g., (2.5) In the denominator we introduced the functionC K/π (t, p) = C K/π (t, p)− 1 2 C K/π (T /2, p) e E K/π (T /2−t) for t < T /2 to subtract the around-the-world contribution due to the periodic boundary condition the mesons experience in the time-direction.
We obtain the vector current renormalisation constant Z V by imposing conservation of the electromagnetic charge, i.e. f ππ/KK (0) = 1 and hence, (2.6) The superscript B in the denominator indicates that we take the bare (unrenormalised) current in the three-point function. While both Z π V and Z K V obtained in this way renormalise the flavour-changing vector current in (2.1) we note that they differ by mass dependent cutoff effects. We will return to this point in later sections when discussing the continuum limit.
As a result of the quantisation of quark-and hadron-momenta in a finite lattice box the desired kinematical situation q 2 = 0 is generally not directly accessible. We follow [6,[15][16][17][18] and use partially twisted boundary conditions for the quark fields, i.e. ψ(x + Lî) = e iθ i /L ψ(x), where L is the spatial extent of the lattice, i one of the spatial directions and ψ one of the up-or strange-quark fields. Varying the twist angle we expect the meson energies to follow the dispersion relation where pL ≡ ∆θ is the difference of twist angles imposed on the valence quark-and antiquark of the respective meson. The kinematical point q 2 = 0 can be reached in all our simulations by suitably tuning the quark boundary conditions [15]. In particular, enforcing we make use of two specific choices for the twist angles [15], As is evident from eqn. (2.2) the vector form factor at zero momentum transfer can be extracted directly from a fit to the three point functions of the scalar current. The vector current matrix element in eqn. (2.1) is instead parameterised in terms of two form factors. These are readily extracted from a simultaneous fit to the correlation function data for both time-and space-components of the vector current.

Simulation parameters
This report centres around two new ensembles with physical light quark masses in large volume [19]. The data analysis will however also take advantage of the information provided by ensembles with unphysically heavy pions which we generated and discussed earlier [20][21][22].
The gauge field ensembles represent QCD with N f = 2 + 1 dynamical flavours at two different lattice spacings. The basic parameters of these ensembles are listed in Table 1. All ensembles have been generated with the Iwasaki gauge action [23,24]. For the discretisation of the quark fields we adopt the domain wall fermion (DWF) action with the Möbius kernel [25][26][27] for the ensembles with physical quark mass and otherwise Shamir DWF [28,29]. The difference between both kernels in our implementation corresponds to a rescaling such that Möbius domain wall fermions are equivalent to Shamir domain wall fermions at twice the extension in the fifth dimension [19]. Möbius domain wall fermions are hence cheaper to simulate while providing the same level of lattice chiral symmetry. Results from both formulations of domain wall fermions lie on the same scaling trajectory towards the continuum limit with cutoff-effects starting at O(a 2 ). Even these O(a 2 ) cutoff-effects themselves are expected to agree between our Mobius and Shamir formulations, with their relative difference at or below the level of 1% [19].
We have used eqn. (2.9) to determine our choice of twist angles. The values for the unphysical ensembles are given in [7]. In our previous studies of the K → π form factor we found that the kinematical situation with the kaon at rest provides for a better signal-tonoise ratio. On the physical point ensembles we therefore only twist the up quark in the pion while leaving the kaon at rest. Our choice is θ T π = 2π L (0.5893, 0.5893, 0.5893) on A phys and θ T π = 2π L (0.5824, 0.5824, 0.5824) on C phys with θ T K = (0, 0, 0) in each case.

Simulation results
In this section we present all steps involved in the analysis towards the determination of the form factor on the physical point ensembles. Apart from details which we mention in the text, the analysis for the unphysical-point ensembles follows [5][6][7].

Correlation functions and AMA
The substantial cost of solving for light quark propagators on the large volume, physical pion mass A phys and C phys ensembles required us to make several algorithmic refinements to our measurement strategy. All correlation functions associated with these ensembles were computed using Coulomb gauge-fixed wall source propagators, together with the allmode averaging (AMA) technique introduced in [30]. In the AMA formalism one replaces a direct calculation of an expensive lattice observable O with a less-expensive approximation O and a correction term ∆O. The lattice action and ensemble averages O , O , and ∆O are all assumed to be invariant under a group G of lattice symmetries. We define the AMA estimator by averaging the inexpensive approximation O over some number N of transformations g ∈ G, and applying the correction term ∆O: where the notation O g denotes O computed after g is applied. We find in practice that the statistical error per unit of computer time can be markedly reduced using AMA with a judicious choice of O and ∆O, relative to computing O directly.
In the context of this calculation the relevant lattice symmetry is the group of translations in the temporal direction. Quark propagators were computed using a deflated mixed-precision conjugate gradient (CG) solver, with 600 (1500) single-precision low-mode deflation vectors obtained from the EigCG algorithm applied to a four dimensional volume source on the A phys (C phys ) ensemble. We further distinguish between exact and sloppy light quark propagators. Exact light quark propagators were computed using a tight CG stopping residual r = 10 −8 for 7 (8) time slices. To avoid bias associated with the even-odd preconditioning used in the CG solves we randomly shifted the time slices used to compute exact propagators on each configuration. Sloppy light quark propagators were computed using a reduced precision r = 10 −4 and for all time slices. Strange quark propagators were sufficiently inexpensive that exact solves were computed for all time slices. For a given two-or three-point function we then constructed a sloppy estimate (O ) for all time slices with the sloppy light quark propagators, and a correction term (∆O) using the exact light quark propagators on time slices for which these are available. We then compute the AMA estimator according to (4.1), after averaging O over all time translations. The full measurement package, which also computes observables related to the K → (ππ) I=2 decay [31] and other low-energy QCD observables [19] from the same propagators, took 5.5 days per configuration on the A phys ensemble using 1 rack (1024 nodes) of IBM Blue Gene/Q hardware, and 5.3 hours per measurement on the C phys ensemble using 32 racks of Blue Gene/Q sustaining 1.2PFlop/s. Additional details of the calculation can be found in [19].
The above set of quark propagators allows generating AMA three-point functions where at least one of the quarks coupling to the external current is a strange quark, for all possible source-sink-separations ∆t = |t f − t i | up to T /2a (e.g. K → π and K → K); in each case based on T /a different positions of the source plane. Results at constant ∆t from different source planes were averaged into one bin on every configuration. For our choice of the 7 (8) source planes for the exact light quark solves on the A phys (C phys ) ensembles the π → π three-point function entering e.g. the determination of Z π V from eqn. (2.6) can be computed on every fourth (fifth) source-sink-separation following the AMA prescription.
In all cases we use the bootstrap resampling technique with 500 samples to determine the statistical errors.

Extracting the form factor
We consider two ways of extracting the form factors: a) simultaneously fit the pion and kaon two-point functions -C π,K (t, p), through the vector and scalar current, as well as their time-reversed counterparts -C Γ,Kπ (t, t snk , p K , p π ) and C Γ,πK (t, t snk , p K , p π ), by minimising a single, global χ 2 . Two-point functions are fit to the right-hand side of eqn. b) determine the result for the ratios R V 0 ,Kπ , R V i ,Kπ and R S,Kπ by taking their value at t = ∆t/2 for even ∆t and the weighted average of the two time-slices adjacent to ∆t/2 for odd ∆t. In this way we avoid having to choose a suitable plateau region subjectively by hand. With the values for the ratios at hand we can solve (2.1) and (2.2) for the vector form factor. In the next section we describe how remaining contaminations by excited states are dealt with.
We find that the values for the form factor obtained from the two analyses agree. All results presented in the following are based on method (b) for the following reasons: we observe as much as a factor of 5 difference in the statistical error on Z π,K V between the ratio fit approach and the global fit approach; for this particular quantity the ratio (2.6) is clearly superior since the measurement of Z π,K V is not contaminated by excited states and hence the operator can be placed closer to the source/sink leading to reduced statistical errors. We also observe that the three-point functions C Γ,Kπ and C Γ,πK are not symmetric between the source and sink walls, since the initial and final states are different, making it a priori difficult to decide on sensible fit ranges for extracting the form factors.

Excited state contamination
The availability of data for a large range of source-sink separations in the three-point functions allows us to study excited states in detail and to reduce their contribution to a minimum. Figure 1 shows a non-trivial dependence on the source-sink separation ∆t of R V 0 ,Kπ , R V i ,Kπ and R S,Kπ as defined in eqn. (2.5). R V 0 ,Kπ is stable starting from surprisingly small values for ∆t. R V i ,Kπ and R S,Kπ on the other hand show a rather noticeable dependence on the source-sink separation. The dependence of R S,Kπ on ∆t can be parameterised very well with a single exponential term added to the constant term which we expect to dominate for large ∆t. This procedure is stable under variation of the smallest and largest values for ∆t included in the fit. While this provides for a good way of determining f Kπ + (0) an alternative approach is to determine from the exponential fit the smallest value of ∆t where excited state contributions are well below the statistical error. Fitting a constant to R S,Kπ above this minimum value for ∆t results in a compatible central value and approximately the same statistical error. The latter fit shows more sensitivity to the range of ∆t included in the fit and we therefore decide to take the result from the exponential parameterisation.
The determination of f Kπ + (0) from R V 0,i ,Kπ proceeds differently since the corresponding ground state matrix element is parameterised in terms of two form factors. In order to take advantage of the data from the large number of source-sink separations at hand we minimise the corresponding global χ 2 -function.
The contamination at small ∆t seen in particular in R V i ,Kπ could in principle be  Table 2. Simulation results that enter the data analysis. For the form factor we consider the results obtained from the vector current matrix element renormalised with Z π V and Z K V , respectively, and we also consider the result determined from the scalar current matrix element. parameterised as above by adding an exponential. This however turned out to result in unstable χ 2 -minimisations. Whilst the ∆t dependence in R S,Kπ was sufficiently pronounced for stable fit results, the combination of a milder ∆t dependence together with a worse signal to noise ratio did not allow us to proceed in this way for R V i ,Kπ . Instead we vary the minimum and maximum values of ∆t entering the above function minimisation for both R V 0 ,Kπ and R V i ,Kπ over a wide range. We were able to find a stable region from which we determine the final results. We illustrate the fit ranges and fit results as shaded bands in the second and third row of plots in figure 1.
The determination of the renormalisation constants Z π,K V proceeds along very similar lines. There is very little dependence of the fit result on ∆t and we determine it by fitting a constant.

Simulation results for the form factor
The valence quark boundary conditions in our simulations were determined using equations (2.9) with estimates for the pion and kaon masses which were computed on a small subset of configurations. The central values of the masses on the full set differ mildly from the ones estimated initially. As a result the value for q 2 may differ significantly from zero. We corrected for this small mistuning as follows: from prior studies [5] we know that the momentum dependence of the form factor is well described by a pole ansatz, with M the pole mass. The latter can be determined from the data for the form factor computed with our choice of twist angle supplemented with the result at q 2 max = (m K − m π ) 2 , i.e. with the pion and kaon at rest. We then inter-or extrapolate to q 2 = 0 and use the corrected data in all subsequent steps of the analysis. Using instead a linear ansatz for the correction leads to a numerically small difference far below the statistical error. This dependence on the parameterisation can therefore safely be neglected. Table 2 shows all accumulated simulation results and in figure 2 we provide a first  visual impression of the data. Both the table and the plot also show the data previously analysed in [5][6][7] where, at variance with what is done here, the geometric mean of the renormalisation constants Z π V Z K V was considered and data for the scalar current matrix element was not considered at all. As can be seen from table 2 the renormalisation constants Z π V and Z K V differ significantly on all ensembles leading to differing results for the vector current matrix element. The form factors renormalised with either choice of Z V differ by mass dependent cutoff effects and hence follow two independent trajectories in the approach to the continuum limit where they are expected to agree by universality.

Corrections towards the physical point
In order to compare our results to experimental measurements of the semi-leptonic kaon decay the data needs to be interpolated to physical quark masses, extrapolated to the continuum limit and finite volume corrections estimated. There are further systematic effects which we will discuss later.
The precise physical process which we are considering is the decay of a neutral kaon into a charged pion, K 0 → π − . While our sea quark masses represent isospin symmetric QCD we approximate the situation found in nature by carrying out a small interpolation of the data in the pion and kaon mass to their values m π − = 139.6MeV and m K = 496.6MeV [32]. The results of this interpolation will then be subject to a continuum extrapolation.
An ansatz for the continuum-extrapolation, mass-interpolation and infinite-volume limit is of the generic form 1) The subscripts NLO and NNLO indicate an expansion equal or similar in spirit to chiral perturbation theory augmented by additional terms that parameterise the cutoff and volumedependence of the lattice data. The case A = 1 corresponds to the SU (3)-symmetric point of the form factor and there are no contributions ∝ (m 2 K − m 2 π ) (Ademollo-Gatto [33]). In a first attempt at finding a suitable model for eqn. (5.1) we try to parameterise the simulation results at constant lattice spacing and defer the extrapolation to the continuum limit to later in this section. In particular, we consider four ansätze for the mass dependence, We have studied these ansätze previously in [7]. They model the basic properties of chiral perturbation theory at NLO [33,34] and NNLO [35], respectively. Fits A and B, which we will refer to as fits based on NLO chiral perturbation theory, adopt the NLO-term [34], where we will employ the leading order relation m η = (4m 2 K − m 2 π )/3 for the η-mass. In fit E and F we replace f 2 by its Taylor expansion in m 2 K − m 2 π . The latter class of fits will be referred to as fits based on polynomial models. In fits B and F the term proportional to A 1 models the the mass dependence expected in the NNLO-term in chiral perturbation theory [35].
In the SU (3)-symmetric limit the form factor is constrained to unity at q 2 = 0, f Kπ + (0) = 1. This relation is exactly reproduced at finite lattice spacing and in finite volume when computed from the matrix element of the vector current renormalised with Z π V or Z K V as defined in eqn. (2.6). The statistical error on the form factor therefore tends towards zero when approaching this limit as can be seen from table 2. Fits are therefore mostly constrained by data points further away from the experimental value of ∆M 2 but close to the SU (3)-symmetric limit. While this can be a virtue if one knows the correct functional form of (5.1) it constitutes a danger when one only has a phenomenological  model at hand. We note that the matrix element of the scalar current does not share this property at finite lattice spacing since cutoff effects can contribute when relating it to the vector current in the discretised theory. In fact, the statistical signal of the form factor as obtained from the scalar current deteriorates compared to the signal we observe for the vector current in the vicinity of the SU (3)-symmetric point. We decided to include data from the scalar current only close to the physical point where the statistical properties are similar to the ones from the vector current.

Fits based on NLO chiral perturbation theory
Here we consider fits A and B where the decay constant f (and in the latter case also A 1 ) are the parameters to be determined. In figure 3 we show representative fits to the form factor renormalised with Z π V to ensembles A and C, respectively. By simple visual inspection and also following from the observation of large values of χ 2 /dof of 4.6 and 1.7, respectively, these ansäzte are at variance with the simulation data. The situation is qualitatively the same when fitting ensembles C, when fitting the other two definitions of the form factor (vector current renormalised with Z K V and scalar current) and also when varying the mass cut-off in any of these cases. Adding a model for the NNLO mass dependence in fit B improves the situation slightly. Large χ 2 /dof values however indicate a strong tension between the ansatz and the data.
Based on these findings we discard fits A and B from the further analysis. While the above data strongly suggests that NLO chiral perturbation theory is unable to describe our results in the range between the SU (3)-symmetric point and the physical quark mass point, it is still possible that the data is better described by the full NNLO-expression [35]. Our data set is however too limited to allow for a study of NNLO fits independent of experimental or model input.  Table 3. Results for global fit E on ensembles A and C for various mass cut-offs. Both the coefficients A and A 0 are allowed to depend on the cutoff. The first column indicates the form factor (vector matrix element (ME) renormalised with Z π V or Z K V , or scalar ME). The superscript in the top line indicates the lattice spacing for which this result was obtained (0.11fm for A and 0.08fm for C).  Table 4. Results for fit F for ensembles A and C for various mass cut-offs. The superscripts on the coefficients A, A 0 and A 1 indicate the set of ensembles for which this coefficient was determined. The first block of data is for the form factor renormalised with Z π V and the next block for the form factor renormalised with Z K V .

Fits based on polynomial models
The results shown in figure 2 show little non-trivial structure in the mass dependence of the form factor when plotted against ∆M 2 . This suggests that a polynomial model should work well in describing the data. We therefore repeat the study of the previous section with fits E and F which are simple polynomials in the SU (3)-breaking parameter ∆M 2 . Some representative plots are shown in figures 4 and 5 and numerical results for various mass cut-offs are summarised in tables 3 and 4. These fits turn out to be of very good quality as indicated by the low value of χ 2 /dof. In the fits with ansatz E the only sensitivity to corrections beyond the linear term (in ∆M 2 ) is indicated by the steep increase in χ 2 /dof when including the data closest to the SU (3)-symmetric point into the fit with the form factor renormalised using Z K V . Indeed, in these cases ansatz F performs better as can be seen in particular by the reduced value of χ 2 /dof when fitting to the more comprehensive  Table 5. Results for fit E for ensembles A and C for various mass cut-offs. The superscripts on the coefficient A indicates the set of ensembles for which this coefficient was determined. The first four data lines are for the form factor as renormalised with Z π V and the next four lines for the form factor renormalised with Z K V . The bottom two lines are for the results from the scalar matrix element.
set of ensembles A. We do not have enough data to allow for a meaningful fit with ansatz F for the scalar form factor. As illustrated in figure 4 (see also table 3), the correction towards the physical point by means of fit E constitutes only a small shift with respect to the results on our physical pion mass ensembles A phys and C phys . The statistical error is almost constant as long as the cut in the pion mass m cut π is chosen well below 600MeV. We conclude that ensembles A phys and C phys play a dominant role in determining the final result and the main role of the ensembles with heavier pion mass is to determine the slope of the interpolation.

Continuum extrapolation
The results for the slope parameter A 0 in table 3 are compatible between ensembles A and C. This indicates that we are not sensitive to a cut-off dependence in the slope with ∆M 2 . In fact, in the O(a)-improved theory we expect the lattice-spacing dependence of the slope-parameter to be A 0 (a) = A 0 (0) 1 + α (aΛ QCD ) 2 + ... , (5.4) with the parameter α parameterising the size of a 2 cut-off effects. Assuming conservatively Λ QCD = 500MeV, (aΛ QCD ) 2 is at most around 8%. Based on these considerations the expected difference in A 0 between the A and C ensembles is of O(3%) and hence smaller than the statistical error on A 0 itself. Our observation that A 0 does not show any significant cutoff effects is therefore not surprising. We confirmed that numerically the parameter α is compatible with zero, in agreement with the observation that the results for A 0 on ensembles A and C are compatible.  f Kπ physical quark masses C-ensembles Figure 5. Illustration for fit F to all data for the form factor renormalised with Z π V . The l.h.s. plot shows the fit to ensembles A, the r.h.s. plot the fit to ensembles C.
Based on these findings we also attempt a modified version of fit E which assumes the slope coefficient A 0 to be independent of the lattice spacing. This fit is illustrated in figure 6 and the numerical results are given in table 5. The plot shows two error bands for ensembles A and C, respectively. The fit is better constrained from ensembles A owing to the larger range of quark masses for which simulation data is available. The fit is of very good quality as indicated by the acceptable values of χ 2 /dof except for the largest mass cut-off in the case of the form factor renormalised with Z K V . We note that when imposing cut-off independence of A 0 the statistical error on the result for the form factor reduces as m cut π is increased. We attribute this to the smaller statistical error on the heavy pion mass ensembles which help constraining the fit. For the most stringent pion mass cut m cut π = 355MeV, though, the statistical error is unaffected by the assumption of cut-off independence of A 0 . In this case the results on A phys and C phys dominate the fit-result.
While our data would allow for taking three independent continuum limits for the form f Kπ physical quark masses A-ensembles C-ensembles Figure 6. Illustration for fit E to all data for the form factor renormalised with Z π V . The coefficient A 0 is assumed to agree for ensembles A and C. Note the two sets of error bands, one for ensemble A and one for ensemble C.
V , C-ensembles S-ME, A-ensembles S-ME, C-ensembles continuum limit Figure 7. Continuum extrapolation for results from fit E with mass cut-off 600MeV. Left: Coefficients A and A 0 differ between ensembles A and C. Right: A 0 assumed to be the same for ensembles A and C. factors as determined from the vector current renormalised with Z π V and Z K V and from the scalar current, respectively, we instead analyse their joint continuum limit assuming universality: We impose that all three extrapolations have to agree in the continuum limit. The combined extrapolation is shown in figure 7 once without and once with the assumption of cutoff independence on A 0 . In table 6 we only show fits for which the χ 2 /dof in the mass interpolation was below one. The result is very stable under variation of the fit ansatz. To underline the stability of our fit ansatz we also show the final result from fits F where either A 1 or A 0 and A 1 are assumed to be cut-off independent. The gain in statistical error from assuming A 0 to be cut-off independent carries over to the continuum limit.  Table 6. Continuum limit results for the form factor based on fit E and F for mass interpolation.

Systematic errors and final result
In this section we explain our choice for the final central value and present the full error budget.
In the previous section we showed a number of approaches for the mass interpolation and eventually combined it with the continuum extrapolation. Table 6 summarises the results for the form factor in the continuum limit based on fits with acceptable values of χ 2 /dof. First of all we see that all variants of the fit lead to compatible results. The major difference between the results is a reduction in statistical error when m cut π is increased while assuming A 0 cut-off independent. As discussed in the previous section this originates from the comparatively small statistical error on the heavy pion ensembles. Choosing m cut π = 355MeV the results from the three fits in table 6 are essentially the same. Without a clear preference we just pick fit E with A 0 fixed as our final result.
We now discuss the remaining sources of systematic errors: finite volume errors: The largest remaining source of systematic uncertainty is expected to be due to finite volume effects (FVE). Since there are only single initial and final states in the K → π matrix element we expect the dominant finite volume effects to be exponentially suppressed with m π L. Moreover, SU (3)-symmetry enforces these effects to disappear exactly in the limit of equal quark masses, also in a finite volume. With the smallest value of m π L = 3.8 on our ensembles we would therefore naively expect finite volume effects on the form factor to be of order 1 − f Kπ + (0) e −mπL = 0.0007. Based on [36] one obtains FVE about twice this number on our ensembles A phys and C phys 1 . However, a conceptually clean calculation of the finite volume effects in chiral perturbation theory including the effect of twisted boundary conditions (see [18]) would be timely and useful in order to improve the quality of this estimate. In the meantime we take twice our naive estimate as the systematic error due to the finite volume. partial quenching: The strange quarks on ensemble A 3 5 and on the C-ensembles are partially quenched, i.e. the mass of the valence strange quark differs mildly from the mass of the sea strange quark. We checked that dropping A 3 5 from the analysis does not change the outcome of our study in any significant way. On the C-ensembles the partial quenching is small and we do not expect significant systematic effects.
isospin symmetry: The unitary light quarks in our simulations are isospin symmetric. We approximate the isospin broken theory by interpolating in the valence sector to the value of ∆M 2 corresponding to the physical point [32]. This still leaves a systematic uncertainty due to the sea-quark isospin breaking which is difficult to quantify in our setup. We expect however that these effects are small compared to the other components of our error budget. Techniques to include such effects in future calculations are being developed [37][38][39][40].

Discussion and conclusions
Simulations of lattice QCD are now feasible with physical light quark masses. This step change in simulation quality leads to the reduction if not removal of the often dominant systematic uncertainty due to the chiral extrapolation. In this paper we have demonstrated how this can be achieved in practice in the case of the K → π form factor at vanishing momentum transfer. This is a phenomenologically important quantity allowing for unitarity tests of the CKM matrix and therefore for stringent constraints of beyond SM physics. Lattice QCD is the only first principles computational tool that can predict this form factor. An important strategic decision that has been made is in which way to make use of our previous results for unphysically heavy light quark masses. We have chosen an intermediate path, i.e. we have used the information from the heavier ensembles to correct for a small mistuning in the average up-and down-quark mass and the strange quark mass to the physical point. Our choice of fit ansatz and fit range gives the result at the physical point the heaviest weight and uses earlier simulation results with heavier pion masses merely for guiding small corrections towards the physical point. In this way any model dependence in the fit ansatz is reduced to a minimum. We note that by restricting the set of ensembles entering the fit less (i.e. including ensembles with heavier pion mass) the statistical error on our final result could have been reduced by around 30%. This would however have come at the cost of an increased model dependence which we find difficult to quantify. The remaining dominant systematic is due to finite volume effects for which we provide an estimate based on effective theory arguments.
Whilst having results with physical quark masses and in the continuum limit represents a huge leap forward there are a number of improvements which we wish to address in the future. As mentioned above the largest systematic uncertainty is now the one due to finite size effects. We would like to understand them better, for example within the framework of partially twisted chiral perturbation theory or by generating results for the form factor on different volumes with otherwise constant simulation parameters (quark masses and lattice spacing). Moreover, in the analysis of experimental results in [3] several approximations were made when averaging the individual decay channels. In particular, electromagnetic and isospin breaking effects were treated within chiral perturbation theory. With the advances made in this work it seems appropriate to rethink this approach and try to treat these effects in a fully nonperturbative fashion in lattice field theory. For an overview of progress in this direction we direct the reader to [37][38][39][40].