$K_{l3}$ form factors at the physical point on (10.9 fm)$^3$ volume

We present the calculation of the $K_{l3}$ form factors with $N_f = 2 + 1$ nonperturbatively $O(a)$-improved Wilson quark action and Iwasaki gauge action at the physical point on a large volume of (10.9 fm)$^3$ at one lattice spacing of $a = 0.085$ fm. We extract the form factors from 3-point functions with three different time separations between the source and sink operators to confirm suppression of excited state contributions. The form factors are calculated in very close to the zero momentum transfer, $q^2 = 0$, thanks to the large volume, so that stable interpolations to $q^2 = 0$ are carried out. Using our form factors, we obtain the form factor at $q^2 = 0$, $f_+(0) = 0.9603(16)(^{+14}_{\ -4})(44)$, where the first, second, and third errors are statistical, systematic errors from fit functions and finite lattice spacing, respectively. The result of $f_+(0)$ yields the Cabibbo-Kobayashi-Maskawa (CKM) matrix element, $|V_{us}| = 0.2255(12)(4)$, where the first error comes from our calculation and the second from the experiment. This value is consistent with the ones determined from the unitarity of the CKM matrix and the $K_{l2}$ decay within one standard deviation, while it is slightly larger than recent lattice calculations by at most 1.6 $\sigma$. Furthermore, we evaluate the shape of the form factors and the phase space integral from our results. We confirm that those results are consistent with the experiment, and also $|V_{us}|$ determined with our phase space integral agrees with the one in the above.


I. INTRODUCTION
Search for signals beyond the standard model (BSM) is an important task in the field of the particle physics. In indirect search for the BSM physics, it is necessary to precisely compare physical quantities obtained from experiments and their predictions in the standard model (SM). Currently one of the indirect searches is carried out through the CKM matrix element |V us |. Its SM prediction is evaluated by the unitarity of the CKM matrix in the first row, i.e., |V ud | 2 + |V us | 2 + |V ub | 2 = 1. Using the precisely determined value of |V ud | = 0.97420 (21) [1, 2], the SM prediction is |V us | = 0.2257 (9), where |V ub | is neglected in the estimate due to |V ub | ≪ |V ud |.
Experimentally, |V us | is related to kaon decay processes, such as the K l2 decay, K → lν, and the K l3 decay, K → πlν, processes, where K, π, l, and ν are the kaon, pion, lepton and neutrino, respectively. In both cases, |V us | is not determined from the experiments only, and lattice QCD calculation also plays an important role, which is the first principle calculation of the strong interaction. For the K l2 decay, the ratio of the decay constants for the kaon and pion, F K /F π , is required to determine the value of |V us |. Using current lattice results, for example F K /F π = 1.1933 (29) in Ref.
[2], |V us | from the K l2 decay well agrees with the SM prediction in the above.
For the K l3 decay, |V us | is related to the decay rate of the K decay Γ K l3 as, where f + (0) is the value of the K l3 form factor at q 2 = 0 with q being the momentum transfer, C K l3 is a known factor including the electromagnetic correction and the SU (2) breaking effect, and I l K is the phase space integral calculated from the shape of the K l3 form factors. The experiment determines Γ K l3 and also I l K . The value of f + (0), however, is not obtained from the experiment. Currently a precise calculation of f + (0) can be performed by lattice QCD.
So far various lattice QCD calculations with dynamical quarks have been carried out to evaluate the value of f + (0) [3][4][5][6][7][8][9][10][11][12]. The most recent study in the N f = 2 + 1 + 1 QCD [10] using the staggered-type quark action reported that using their value of f + (0) there is a clear deviation of |V us | in more than 2 σ from the ones in the SM prediction and the K l2 decay. Another N f = 2 + 1 + 1 calculation using the twisted quark action [9] obtained a similar result. Therefore, it is an urgent task for the search for BSM signals to confirm those results by several lattice calculations using different types of quark actions with small statistical and systematic uncertainties.
For this purpose we calculate the K l3 form factors in the N f = 2 + 1 QCD using a nonperturbatively O(a)-improved Wilson quark action and Iwasaki gauge action at one lattice spacing of a = 0.085 fm. The gauge configurations employed in this work are a subset of the PACS10 configurations [13]. Our calculation is carried out at the physical light and strange quark masses, and on a larger physical volume of (10.9 fm) 3 than typical current lattice QCD calculations. Thus, our result significantly suppresses the uncertainties coming from the chiral extrapolation and finite volume effect. Another advantage using the large volume is that it is possible to access small q 2 region without resort to the twisted boundary condition. Thanks to the large volume, one of our data is very close to q 2 = 0, so that we can perform reliable interpolations of the form factors to q 2 = 0. Using our result of f + (0), we determine |V us | and compare it with the SM prediction, the K l2 decay, and also the previous lattice QCD results. Furthermore, since we calculate the form factors in a wide range of q 2 , the shape of the form factors and also the phase space integral are successfully evaluated from our results. Those values are compared with the experiment and previous lattice QCD results. We also determine |V us | using our result of the phase space integral.
Our preliminary result has been already reported in Ref. [14]. This paper is organized as follows. Section II explains our calculation method of the form factors from meson 2-and 3-point functions. In Sec. III simulation parameters and technical details of our calculation are presented. The result of the form factors is shown in Sec. IV. The interpolations of the form factors are discussed in Sec. V. Our results for f + (0), and the shape of the form factors, and phase space integral are also presented in this section. The result of |V us | and its comparison with other determinations are discussed in Sec. VI. Section VII is devoted to conclusion. Appendices explain interpolating functions based on the SU(3) chiral perturbation theory (ChPT) used in our analysis and tables for some interpolation results.

II. CALCULATION METHOD
The K l3 form factors f + (q 2 ) and f − (q 2 ) are defined by the matrix element of the weak vector current V µ as, where q = p K − p π is the momentum transfer. The scalar form factor f 0 (q 2 ) is defined by f + (q 2 ) and f − (q 2 ) as, where ξ(q 2 ) = f − (q 2 )/f + (q 2 ), and m π and m K are the masses for π and K, respectively. At q 2 = 0, the two form factors f + (q 2 ) and f 0 (q 2 ) give the same value, f + (0) = f 0 (0).
The K l3 form factors are calculated from 3-point function C Kπ µ ( p, t) with the weak vector current given by where We use only the periodic boundary condition in the spatial directions for quark propagators in contrast to the recent calculations of the K l3 form factors using the twisted boundary condition [7,9,10], because the spatial extent L in our calculation is large enough to obtain the form factors near the q 2 = 0 region. Thus, p is labeled by an integer vector n p with p = | p| as p = (2π/L) n p . While we have also calculated the 3-point functions with moving K and π at rest, their form factors are much noisier than the ones from C Kπ µ ( p, t). Therefore, we will not discuss those results in this paper.
For the renormalization factor of the vector current Z V , we compute 3-point functions for the π and K electromagnetic form factors at q 2 = 0 in a similar way, which are given by where V em 4 (t) is the temporal component of the electromagnetic current. The 2-point functions for π and K are calculated as, We average the 2-point functions with the periodic and anti-periodic boundary conditions in the temporal direction to make the periodicity in the temporal direction of the 2-point functions effectively doubled. The asymptotic form of C X ( p, t) for X = π and K in the t ≫ 1 region is given by with E X (p) = m 2 X + p 2 and the temporal extent T . The mass m X and amplitude Z X are obtained from a fit of C X ( 0, t) in a large t region with the asymptotic form.
The matrix element in Eq. (2) is obtained from the ground state contribution of C Kπ µ ( p, t), which needs to avoid excited state contributions by investigating time dependences of C Kπ µ ( p, t). To do this, we define a ratio R BC µ ( p, t), which has the following time dependence as, where N 4 ( p) = 1 and N i ( p) = 1/p i with i = 1, 2, 3, and Z π and Z K are defined in Eq. (12).
C Kπ µ,BC ( p, t) is the 3-point function in Eq. (4) with the (anti-)periodic boundary condition in the temporal direction, which is represented by BC = (A)PBC in the following. In Eq. (14) it is assumed that t i ≤ t ≤ t f and two excited state contributions for the radial excited mesons and wrapping around effect, expressed by ∆A µ ( p, t) and ∆B µ ( p, t), respectively, are leading contributions of excited states in the ratio. Other excited state contributions are denoted by the dots (· · · ) term. The sign of the wrapping around effect in the temporal direction depends on the temporal boundary condition of C Kπ µ,BC ( p, t), i.e., b PBC = 1 and b APBC = −1, because in the wrapping around contribution one of the mesons in C Kπ µ,BC ( p, t) crosses the temporal boundary. A similar wrapping around effect was discussed in the 3-point function of the B K calculation [15].
The time dependence of the two excited state contributions is given by where E ′ π (p) = (m ′ π ) 2 + p 2 , and m ′ X is the mass of the radial excitation of X = π and K. In the second equation, we assume that the finite volume effect in the energy of the πK scattering state is negligible in our volume. In a small p, the first term of the right hand side in Eq. (16) has a non-negligible effect in R BC µ ( p, t), which will be presented later. This is because m π T = 7.5 is not enough to suppress the wrapping around effect at the physical m π . We remove the wrapping around effect ∆B µ ( p, t) by averaging the ratios R PBC µ ( p, t) and R APBC µ ( p, t). On the other hand, another excited state contrition ∆A µ ( p, t) remains in the averaged ratio. This contribution needs to be removed, and it will be discussed in a later section.

III. SET UP
We use the configurations generated with the Iwasaki gauge action [16] and the stoutsmeared Clover quark action at the physical point on (L/a) 3 × T /a = 128 3 × 128 lattice corresponding to (10.9 fm) 4 . These configurations are a subset of the PACS10 configurations. Parameters for the gauge configuration generation are found in Ref. [13].
The bare coupling β = 1.82 corresponds to a −1 = 2.3162(44) GeV [17] determined from the Ξ baryon mass input. The hopping parameters for the light and strange quarks are (κ l , κ s ) = (0.126117, 0.124902), and the coefficient of the clover term is c SW = 1.11, which is nonperturbatively determined in the Schrödinger functional (SF) scheme [18]. It is employed the six-stout-smeared link [19] with ρ = 0.1 in the quark actions. We use the same quark actions for the measurement of the K l3 form factors. The measured π and K masses, m π = 0.13511(72) GeV and m K = 0.49709(35) GeV, in this calculation are consistent with the ones in our spectrum paper [13].
The measurements for the 2-point and 3-point functions are performed using 20 configurations separated by 10 molecular dynamics trajectories. To reduce the calculation cost of the measurements, the Z(2)⊗Z(2) random source [20] at the source time slice t i is employed, where random numbers are spread in the color and spin spaces as well as the spatial volume. For example, the operator O π ( p, t) at the source time slice t i in Eq. (4) is replaced by where N r is the number of the random source, and the color and spin indices are omitted.
The Z(2) ⊗ Z(2) random source η j ( x) satisfies the following condition as, We use the sequential source method at the sink time slice t f in C Kπ µ,BC ( p, t). The quark propagators are calculated with the periodic boundary condition in the spatial and also temporal directions in C Kπ µ,PBC ( p, t). On the other hand, in C Kπ µ,APBC ( p, t), though the spatial boundary condition is periodic for all the quark propagators, one of the three quark propagators needs to be calculated with the anti-periodic boundary condition in the temporal direction. In this work we choose the quark propagator which connects the source operator with the sink one. This choice is suitable for our purpose to remove the wrapping around effect ∆B µ ( p, t) in Eq. (14), because in this case the effect has a desirable boundary condition dependence as in Eq. (14).
In each momentum, the quark propagator of the random momentum source corresponding to the first square brackets in Eq. (17) is calculated. To improve statistical error of C Kπ µ,BC ( p, t) in a finite momentum, we average C Kπ µ,BC ( p, t) in each n p = (Lp/2π) 2 with several momentum assignments. The number of the momentum assignment is listed in Table I together with the values of q 2 = −(m K − E π (p)) 2 + p 2 calculated using the measured m K and m π with E π (p) = m 2 π + p 2 . The value of q 2 for each n p is labeled by q 2 np in the following. We vary the time separation between the source and sink operators, t sep = t f −t i = 36, 42, and 48, corresponding to 3.1, 3.6, and 4.1 fm in physical unit, to study the excited state contributions of R BC µ ( p, t) in Eq. (13). Since the statistical error increases for larger t sep , the number of the random source N r = 2 is chosen in the t sep = 42 and 48 cases, while N r = 1 in t sep = 36.
In order to increase statistics effectively, on each configuration we perform the measurements with 8 different t i equally separated by 16 time separation, 4 temporal directions by rotating the configuration, and also average C Kπ µ,BC ( p, t) with its backward 3-point function calculated in t f ≤ t ≤ t i with the same t sep . In total the numbers of the measurements are 2560 for t sep = 36 and 5120 for t sep = 42 and 48, where the different choice of N r explained above is included. The statistical errors for all the observables are evaluated by the jackknife method with the bin size of 10 trajectories.

IV. K l3 FORM FACTORS
In this section we discuss the two kinds of excited state contributions in the ratio of the 3-point function R BC µ ( p, t) as explained in Sec. II, which are the wrapping around effect and the radial excited state contributions. We also present the results for the K l3 form factors, f + (q 2 ) and f 0 (q 2 ). In the following discussions, we choose t i = 0 so that t sep = t f .

A. Wrapping around effect
A typical example of the wrapping around effect in R PBC 4 ( p, t) and R APBC 4 ( p, t) defined in Eq. (13) is shown in Fig. 1, where the ratios with t sep = 42 at q 2 = q 2 0 are plotted. A clear discrepancy between R PBC 4 ( p, t) and R APBC 4 ( p, t) is observed in the region of t > t sep /2, where the first term in the right hand side of Eq. (16) is expected to have a large contribution.
The averaged ratio, has a milder t dependence than the two ratios, because the wrapping around effect ∆B µ ( p, t) in Eq. (14) cancels in the average.
Since the effect decreases as p 2 increases expected from Eq. (16), the discrepancy between the two ratios becomes smaller at q 2 = q 2 1 as shown in Fig. 2. Although the effect is small in large p 2 , we always adopt the averaged ratio R µ ( p, t) in the following analyses.   Figure 3 shows t sep dependence of the averaged ratio R µ ( p, t) in Eq. (19) at q 2 = q 2 1 . For R 4 ( p, t), we observe a reasonable consistency of the data with the different t sep in flat regions between the source at t = 0 and sink at t = t sep . For R i ( p, t), flat regions are shorter than those of R 4 ( p, t), and their shapes are non-symmetric. In contrast to R 4 ( p, t) the central values of R i ( p, t) in the flat region of t =10-15 with t sep = 36 are about 1% smaller than those with t sep = 42 and 48. We consider that it is caused by excited state contributions in R i ( p, t), and at first assume that it is the radial excitation of the mesons as explained in Sec. II.
To remove the contribution and extract the matrix element π( p) |V µ | K( 0) corresponding to the constant part in R µ ( p, t), we fit R µ ( p, t) with a fit form given by, where R µ (p), A π µ (p), and A K µ (p) are fit parameters, and E ′ π (p) = (m ′ π ) 2 + p 2 . Since our simulation is carried out at the physical point, the masses for the radial excited mesons, m ′ π and m ′ K , are fixed to the experimental values m ′ π = 1.3 GeV and m ′ K = 1.46 GeV in PDG18 [2]. We examine if these masses are appropriate in our calculation by effective masses for the first excited states in the 2-point functions. The effective mass is evaluated from C X ( 0, t) without the ground state contribution defined as where with Z X and m X obtained from a fit using the asymptotic form in Eq. (12). As shown in Fig. 4, we observe that the effective masses for the first excited states show reasonable consistencies with these experimental values.
The simultaneous fit results using all the t sep data with the fit form in Eq. (20) are presented for µ = 4 and i at q 2 1 in Fig. 3. We employ uncorrelated fits in this analysis, because our statistics is not enough to determine the covariance matrix precisely. The minimum time slice t min of the fit range is fixed for all t sep , while the maximum time slice t max is changed for each t sep as t max = t sep − t fit . In the q 2 = q 2 1 case as shown in the figure, (t min , t fit ) = (7, 12) and (6,18) are chosen for µ = 4 and i, respectively. The fit result of R µ (p) represented by the shaded band in the figure agrees with the data in the flat region with the larger t sep in both cases.
Although the above fit using the experimental m ′ π and m ′ K works well in our data, their contributions might not be the leading excited state contributions in the ratios. In order to test the possibility, we also fit R µ (p) with E ′ π (p) and m ′ K as fit parameters in the fit form Eq. (20), and compare the results from the two analyses. In this case we can choose wider fit ranges as (t min , t fit ) = (5, 7) and (4, 14) for µ = 4 and i, respectively, than the ones with the experimental m ′ π and m ′ K . The fit curves are presented in Fig. 5. The results of R µ (p) agree with those in Fig. 3, although the error of R i (p) becomes larger.
For later convenience, we call the data obtained from the fit with the fixed m ′ π and m ′ K as "A1", and ones from another fit as "A2". In the following, we use these two data to estimate a systematic error originating from the choice of the fitting form. In each q 2 , we carry out similar analyses to obtain R µ (p) for µ = 4 and i, except for at q 2 = q 2 0 where only R 4 (p) is available.
We also perform the same analysis without the data of t sep = 36 to study effects from the smallest t sep data in our analysis. It is found that the effect is not significant in our result, because the fit result agrees with the above ones within the error.

C. Form factors
For the renormalization of the local vector current, the renormalization factor Z V is calculated from the 3-point functions for π and K with the electromagnetic current as presented in Eqs. (8) and (9). To determine Z V , a ratio R Z V (t) is defined as whose value in a plateau region corresponds to Z V . Figure 6 shows that the data of R Z V (t) with t sep = 36 and 48 agree with each other. Thus, we determine Z V from a constant fit with R Z V (t) of t sep = 36 in the middle t region of 10 ≤ t ≤ 24. The result of Z V = 0.95587 (18) is 0.45% larger than the value obtained by the SF scheme [21], Z SF V = 0.95153(76), which is also shown in the figure. From the discrepancy we will estimate a systematic error of the form factors in a later section.
Combining Z V , Z π and Z K from the 2-point functions, and the results for R 4 (p) and R i (p), we calculate the matrix elements π( p) |V 4 | K( 0) and π( p) |V i | K( 0) /p i , and then evaluate f + (q 2 ) and f 0 (q 2 ) with Eqs. (2) and (3) at each q 2 , except at q 2 = q 2 0 where only f 0 (q 2 0 ) is obtained. The form factors f + (q 2 ) and f 0 (q 2 ) obtained from the two data sets, A1 and A2, are plotted in Fig. 7 as a function of q 2 together with the ratio of the form factors ξ(q 2 ) defined in Eq. (3). Their numerical values are presented in Table II. The results for f + (q 2 ) and f 0 (q 2 ) with the A1 data, which are taken to be the central values in our analysis, are obtained within less than 0.3% statistical error.

V. q 2 DEPENDENCE OF FORM FACTORS
It is necessary for the determination of |V us | to extract f + (0) from f + (q 2 ) and f 0 (q 2 ).
Although in our calculation q 2 4 is very close to zero, we need a small interpolation using some fit function. In this section we explain the interpolation procedures and give the results for f + (0) with systematic errors. All the interpolations are carried out with uncorrelated fits due to a lack of enough statistics to determine a precise covariance matrix. Furthermore, we also discuss the shape of the form factors and the phase space integral evaluated with our form factors. A. Interpolations to q 2 = 0 For the form factors, f + (q 2 ) and f 0 (q 2 ), the next-to-leading order (NLO) formulae are available in the SU(3) ChPT [22,23]. We employ the following fit functions for the interpolations to q 2 = 0, which are based on the NLO ChPT formulae: where F 0 is the pion decay constant 1 in the chiral limit, and L 9 , L 5 , c 0 , and c +,0 2 are free parameters. The constraint of f + (0) = f 0 (0) is required, so that the same c 0 appears in 1 We adopt the normalization of F π ∼ 132 MeV at the physical point. the two fit functions. The functions K + (q 2 ) and K 0 (q 2 ) are given in Appendix A, which depend on m π , m K , q 2 , F 0 , and the scale µ. In our analyses we fix µ = 0.77 GeV. The last two terms in each fit function can be regarded as a part of the NNLO analytic terms in the SU(3) ChPT. In this point of view, the constant term c 0 represents a sum of m 4 π , m 4 K , and m 2 π m 2 K terms, because m 2 π and m 2 K are constant in our analysis due to the physical point calculation.
We also carry out another fit using the same fit forms of Eqs. (24) and (25), while setting in the above fit. This observation indicates that a systematic error due to the fixed F 0 should be small in the above fit. The fit results are summarized in Table III. The table   also contains the fit results using the A2 data. The results of f + (0) with the A2 data agree with the ones with the A1 data, while they have larger errors than those with the A1 data.
Our results for L 9 and L 5 show similar values to the previous lattice QCD results for the low energy constants in the SU(3) ChPT summarized in the FLAG review 2019 [26], i.e., L 9 ∼ (2.4-3.8)×10 −3 and L 5 ∼ (0.9-1.5)×10 −3 . Note that these values are given in the chiral limit for all the quark masses, so that they cannot be directly compared with our results obtained at the physical quark masses.
We also employ several different fit forms for the interpolation, such as a mono-pole function, a simple quadratic function of q 2 , and variations of the z-parameter expansion [28].
The fit forms and the fit results are summarized in Appendix B. The results of f + (0) obtained from these fits are also consistent with the ones obtained from the above ChPT analyses.
Furthermore, we confirm that the result is not changed, when the fit range of q 2 is squeezed as q 2 2 ≤ q 2 ≤ q 2 5 in the NLO ChPT fit with the fixed F 0 using the A1 data, which gives f + (0) = 0.9604 (16). Based on these fit analyses we conclude that the systematic error originating from the fit form dependence for the interpolation is as small as the statistical error in our result of f + (0).

B. Result of f + (0)
From the fit results discussed in the previous subsection, whose values are tabulated in Table III and Tables in Appendix B, we obtain the result of f + (0) as where the central value and statistical error (the first error) are determined from the fit result based on the ChPT formulae in Eqs. (24) and (25) with the fixed F 0 using the A1 data. The second error is the systematic one for the fit form dependence, which is estimated from the deviation of the various fit results from the central value. The third error is the systematic one for the discrepancy of Z V and Z SF V , 0.45%, discussed in Sec. IV C. We consider that it is regarded as order of a systematic error due to the finite lattice spacing effect, because the discrepancy should vanish in the continuum limit. In our calculation this estimate is more conservative than the one employed in the previous studies [5,7], where the finite lattice We do not include the systematic error of the finite volume effect, because our physical volume is large enough to suppress the effect. The estimate based on ChPT, (1 − f + (0)) × e −mπ L , gives 0.002%, which is much smaller than other errors. In the following we will not discuss this systematic error.  (3) ChPT formulae in Eqs. (24) and (25) together with the value of the uncorrelated χ 2 /dof. A1 and A2 data sets are explained in Sec. IV B. fit-1 and fit-2 denote fits with the fixed and free F 0 , respectively, as explained in Sec. V A. We also list the results for f + (0), slope, curvature, and phase space integral.   Figure 9 shows comparison of our result with the previous dynamical lattice QCD calculations [3][4][5][6][7][8][9][10][11][12]. Our result is reasonably consistent with the previous N f = 2 [11,12] and N f = 2 + 1 [3][4][5][6][7] calculations, while it is slightly smaller than the recent N f = 2 + 1 + 1 results [8][9][10]. The largest discrepancy in comparison with the previous results is 1.8 σ from the one in Ref. [10] in the total error. At present the reason of the discrepancy is not clear.
However, an analysis using only the physical point data in Ref. [10] gives a smaller value than their result in the figure, so that the discrepancy would become smaller with larger systematic errors [10]. In order to understand the source of the discrepancy, it is important to reduce our uncertainties, especially, the finite lattice spacing effect, which is the largest error in our calculation. For this purpose, in the next step we will calculate the form factors using other sets of PACS10 configurations with the finer lattice spacings at the physical point to evaluate f + (0) in the continuum limit.  [3,8,10], twisted [9,12], overlap [7], and domain wall [4][5][6]11] quark calculations, respectively. The filled symbols represent the results at a finite lattice spacing. The inner and outer errors express the statistical and total errors. The total error is evaluated by adding the statistical and systematic errors in quadrature.

C. Shape of form factors
The slopes for the form factors are defined by the Taylor expansion in a vicinity of q 2 = 0 as, where s = + and 0, and m π − = 0.13957061 GeV.
The fit of the form factors discussed in Sec. V A gives the slope λ ′ s and curvature λ ′′ s , whose results are presented in Table III and Tables in Appendix B. For λ ′ s , we obtain where the central value, the first and second errors are determined in a similar way to the f + (0) case shown in Sec. V B. Since the systematic error coming from Z V affects only the overall constant f + (0), there is no corresponding systematic error for λ ′ s and λ ′′ s . Those results are well consistent with the experimental ones [29] 3 , λ ′ + = 2.575(36) × 10 −2 and λ ′ 0 = 1.355(71) × 10 −2 , and the previous lattice ones [7,9,12] as shown in Fig. 10.
For λ ′′ s , we obtain where the results and errors are determined in a similar way to λ ′ s . Those results agree with the experimental ones calculated with the dispersive representation [30] [7,9,12]. The experimental results [29] are denoted by the shaded bands.

D. Phase space integral
Since our results for the slopes and the curvatures of the form factors agree with the experiment, we evaluate the phase space integral, I l K in Eq. (1), which is usually calculated using the q 2 dependence of the experimental form factors. The phase space integral [32] is given by f s (−t)/f + (0) with s = + and 0, t = −q 2 , and m l is the mass of the lepton l. Substituting the fit results for the form factors into the equation, we calculate I l K for the K 0 → π − e + ν e and K 0 → π − µ + ν µ processes and obtain each integral as, I e K 0 = 0.15476 (18) and I µ K 0 = 0.10253 (16), in Ref. [31]. We also show the results for the unnormalized phase space integrals as, which will be used for evaluation of |V us | in the next section. The first error represents the statistical one. The second and third errors are the systematic ones due to the choice of the fit forms and the finite lattice spacing effect, respectively. Their numerical values for each fit form are summarized in Table III and Tables in Appendix B.

VI. RESULT OF |V us |
Using our result of f + (0) in Sec. V B and the experimental value |V us |f + (0) = 0.21654(41) [29], the result of |V us | in our study is given by where the first, second, and third errors are the statistical error, systematic error of the fit form in our analyses, and systematic error of the finite lattice spacing effect estimated from the difference of Z V and Z SF V , respectively. The fourth error comes from the experimental one. The result can be expressed as |V us | = 0.2255(12)(4), where the first error is given by the combined error of the three errors in our calculation. Figure 11 shows that our result is consistent with the value estimated by assuming the unitarity condition of the first row of the CKM matrix: where |V ub | is neglected due to |V ub | ≪ |V ud | and we use |V ud | = 0.97420 (21) [1, 2]. Furthermore, our result agrees with the results determined from the K l2 decay process through |V us |/|V ud | × F K /F π = 0.27599(38) [2]. In the figure we plot two data: one is obtained with the use of the value of F K /F π in PDG18 [2] and the other is from the result of F K /F π calculated with the same configuration as in this work [13]. These observations suggest that our result is consistent with the SM prediction within the error. Using a new evaluation of |V ud | = 0.97370(14) [33], however, the value from the unitarity condition significantly changes as |V us | = 0.2278 (6), while the ones from the K l2 decay do not move within the error. In this case, our result is smaller than the unitarity condition by 1.7 σ. More recent evaluation of |V ud | = 0.97389(18) [34] leads to the unitarity value of |V us | = 0.2270 (8), which is consistent with our result within 1.1 σ. Figure 11 also presents the comparison of our result with the recent N f = 2 + 1 [5][6][7] and N f = 2 + 1 + 1 [9,10] calculations. Our result is reasonably consistent with all the results, although it is 1.6 σ larger than the recent result in Ref. [10] as for the result of f + (0) presented in Sec. V B. To understand the difference, it is an important future work to reduce the uncertainties in our calculation. We will remove our largest systematic error by measuring the form factors at a finer lattice spacing in the next calculation. We also determine |V us | using the phase space integral calculated with our form factors, where we use |V us |f + (0) I e K 0 = 0.08510 (22) and |V us |f + (0) I µ K 0 = 0.06935 (21), which are evaluated from the experimental results and correction factors in Ref. [31]. The meaning of the errors is the same as in the above |V us |. A weighted average of the two decay processes using the experimental error gives This value is well consistent with that in Eq. (37) including the sizes for the uncertainties.
It is encouraging that the K l3 form factors calculated in the lattice QCD can be used for not only the determination of f + (0), but also the evaluation of the phase space integral. To obtain the value of the form factors at q 2 = 0, which is essential to evaluate |V us |, we have interpolated the form factors to q 2 = 0 employing several fit forms. These interpolations also contribute to determination of the shape of the form factors as a function of q 2 . The chiral extrapolation is not necessary in our analysis thanks to the calculation at the physical point. The central value of the form factor at q 2 = 0 is determined with less than 0.2% statistical error. Since one of our data is very close to q 2 = 0, the interpolations are fairly stable, and the systematic error from the interpolations is as small as the statistical error.
The largest systematic error in our result comes from a finite lattice spacing effect, which is estimated from the different determination of Z V . The final result of f + (0) in this work is where the first error is statistical and the second and third errors are the systematic ones for the choice of the fit forms and the finite lattice spacing effect, respectively. Our result is reasonably consistent with the recent N f = 2 + 1 and N f = 2 + 1 + 1 QCD calculations. The largest deviation from the recent results is 1.8 σ. It is important to reduce the uncertainties in our calculation to understand the source of the deviation. Thus, in the next calculation, we will measure the form factors at a finer lattice spacing.
The slope and curvature for the form factors at q 2 = 0 are determined from the inter- where the first three errors inherit those of f + (0). The fourth error comes from the experimental one. This result agrees with |V us | determined from the unitarity condition of the CKM matrix and also from the K l2 decay. On the other hand, using a new evaluation of |V ud |, our result differs from the unitarity value by 1.7 σ. To make the comparison with the SM predictions more stringent, we would need to reduce the uncertainties in our calculation.
Thus, our next calculation with the finer lattice spacing is important also from this point of view.
It is encouraging that another determination of |V us | using the phase space integral evaluated with our form factors completely agrees with the above conventional determination.
This suggests that the lattice calculation could contribute to not only a precise determination of f + (0), but also the phase space integral. with t = −q 2 and and variations of the z-parameter expansion [28], c + i z i (q 2 ) and f 0 (q 2 ) = f + (0) + Nz i=1 c 0 i z i (q 2 ) (B3) where N z = 1 or 2 and z(q 2 ) = (m K + m π ) 2 + q 2 − (m K + m π ) (m K + m π ) 2 + q 2 + (m K + m π ) . (B4) Our choice of z(q 2 ) corresponds to the one with t 0 = 0 in the general representation of z(q 2 ) [28]. The parameters f + (0), c +,0 i , and M +,0 are fit parameters. These fit results are summarized in Tables IV-VI.   with N z = 1 and 2.