$B_s\to K\ell\nu$ decay from lattice QCD

We use lattice QCD to calculate the form factors $f_+(q^2)$ and $f_0(q^2)$ for the semileptonic decay $B_s\to K\ell\nu$. Our calculation uses six MILC asqtad 2+1 flavor gauge-field ensembles with three lattice spacings. At the smallest and largest lattice spacing the light-quark sea mass is set to 1/10 the strange-quark mass. At the intermediate lattice spacing, we use four values for the light-quark sea mass ranging from 1/5 to 1/20 of the strange-quark mass. We use the asqtad improved staggered action for the light valence quarks, and the clover action with the Fermilab interpolation for the heavy valence bottom quark. We use SU(2) hard-kaon heavy-meson rooted staggered chiral perturbation theory to take the chiral-continuum limit. A functional $z$ expansion is used to extend the form factors to the full kinematic range. We present predictions for the differential decay rate for both $B_s\to K\mu\nu$ and $B_s\to K\tau\nu$. We also present results for the forward-backward asymmetry, the lepton polarization asymmetry, ratios of the scalar and vector form factors for the decays $B_s\to K\ell\nu$ and $B_s\to D_s \ell\nu$. Our results, together with future experimental measurements, can be used to determine the magnitude of the Cabibbo-Kobayashi-Maskawa matrix element $|V_{ub}|$.


I. INTRODUCTION
Semileptonic decays of hadrons can be used to determine elements of the Cabibbo-Kobayashi-Maskawa (CKM) matrix. However, since the quarks that participate in the underlying electroweak transition are constituents of bound states, it is necessary to understand the effects of the strong interactions on the decay. These effects are encapsulated in form factors for hadronic matrix elements of the weak currents that govern the decay.
Lattice QCD has allowed us to calculate the form factors with increasing precision, making possible stringent tests of the Standard Model and the CKM paradigm. Should there be a violation of unitarity of the CKM matrix, or should two decay processes that depend on the same CKM matrix element imply different values for that CKM matrix element, we would have evidence for physics beyond the Standard Model. The decay studied here, B s → K ν depends on the same matrix element V ub as the decay B → π ν. Indeed, the only difference between the two decay processes is that the light spectator up (u) or down (d) quark in the latter process is replaced by a strange (s) quark in the case at hand. Since in lattice QCD, strange quarks generally yield smaller statistical errors and are easier to deal with computationally, a lattice calculation of the form factors for B s → K ν decay can enable a precise |V ub | determination. This, in turn, can provide a useful test of |V ub | determinations quarks. Other previous calculations of the B s → K ν decay form factors are based on the relativistic quark model [33], light-cone sum rules [34,35], and next-to-leading-order (NLO) perturbative QCD [36]. In Sec. VI C, we compare our results with the prior results. Preliminary reports on this study can be found in Refs. [37] and [38], where the vector current renormalization factors were still multiplied by a blinding factor. This factor was disclosed only after the analysis was finalized.
The rest of this paper is organized as follows. In Sec. II, we define the continuum decay form factors and the hadronic matrix elements needed to calculate them. In Sec. III, we introduce the lattice QCD operators and the form factors most convenient to calculate on the lattice. We detail how to calculate the needed lattice matrix elements and enumerate the MILC asqtad 2+1 flavor ensembles we have used. Section IV discusses our analysis of the two-and three-point functions needed to construct the lattice form factors. We also explain how we take the chiral-continuum limit. Section V contains our analysis of systematic errors in the range of momentum transfer accessible in our calculation. To construct the continuum form factors over the entire range of momentum transfer, we present the functional z expansion in Sec. VI. We then apply it to obtain our final results for the form factors. Section VII presents some of the phenomenological implications of the results. Appendix A contains details of our application of SU(2) chiral perturbation theory to perform the chiral extrapolation in Sec. IV. Appendix B details how we construct the continuum form factors in Sec. VI B. Appendix C contains the binned differential decay rates, as well as the full correlation matrices.

II. MATRIX ELEMENTS AND FORM FACTORS
To lowest order in the weak coupling constant, the semileptonic B s → K ν decay can be described via the Feynman diagram shown in Fig. 1. The relevant hadronic matrix element form factors corresponding to the exchange of 1 − and 0 + particles. These two form factors are subject to a kinematic constraint: f + (0) = f 0 (0), (2.2) which eliminates the spurious pole at q 2 = 0 in Eq. (2.1). The tensor form factor f T parametrizes the hadronic matrix element of the tensor current T µν = iūσ µν b. Since it does not contribute to the Standard Model decay rate, we do not include it in this calculation.
In the Standard Model, the angular-dependent differential decay rate for the B s → K ν can be written as in the B s meson rest frame. Here G F is the Fermi constant, V ub is an element of the CKM matrix, m is the lepton mass, and θ is the angle between the final charged-lepton and the B s meson momenta in the rest frame of the final state leptons. Thus, to determine |V ub | from a measurement of the differential decay rate, it is necessary to compute the form factor f + (q 2 ). If the charged lepton is the τ , however, the lepton mass cannot be neglected and f 0 (q 2 ) is also necessary.

III. LATTICE-QCD CALCULATION
In this section, we present the ingredients of our lattice-QCD calculation. The square of the lepton momentum transfer q 2 can then be expressed as where E K = p K · v is the kaon energy. Defining as the projection of the kaon momentum in the direction perpendicular to v µ and using Eq. and f ⊥ (E K ) as The relations to the original form factors f + and f 0 are given by The kinematic constraint, Eq. (2.2), is automatically satisfied in Eq. (3.5).
In the B s rest frame, which we use throughout the lattice-QCD calculation, the form factors f and f ⊥ are related to the temporal and spatial components of the matrix element of the vector current V µ via Note that there is no summation over the superscript i in Eq. (3.6b). The continuum-QCD current is related to the lattice current operator V µ by a multiplicative renormalization factor, i.e., The lattice current V µ is defined in Sec. III C, below. We use a mostly nonperturbative method to compute Z Vµ . The details are explained in Sec. III C.
The desired matrix elements (and hence form factors) can be calculated from suitably defined two-and three-point correlation functions: where O Bs and O K are lattice interpolating operators, which are defined in Sec. III C, below.
Further, p K is the kaon spatial momentum, whose components in a finite volume are integer multiples of 2π/N s , where N s is the lattice spatial dimension in lattice units.
The basic procedure for calculating the continuum form factors f + and f 0 in Eq. (2.1) in lattice QCD is the following: 1. For each ensemble: (i) Determine the lattice B s meson masses, kaon masses and energies from the lattice two-point correlation functions.
(ii) Determine the lattice form factors f lat and f lat ⊥ at several discrete kaon momenta p K from the two-and three-point correlation functions.
(iii) Obtain the renormalized form factors by matching the lattice current to the continuum as in Eq. (3.7).
2. Use chiral perturbation theory together with Symanzik effective theory to perform a combined chiral-continuum fit to the renormalized form factors and extrapolate them to the physical quark masses and continuum (zero lattice spacing) limits. This yields the continuum form factors f and f ⊥ as functions of the kaon recoil energy E K in the interval covered by the simulation, roughly 0.5 GeV E K 1 GeV.
3. Construct the continuum form factors f + and f 0 from f and f ⊥ via Eq. (3.5) and employ a z expansion to parametrize their shapes and to extrapolate them from the low-recoil range to the entire kinematically allowed region, which extends at high recoil to q 2 = 0.
The asqtad fermion action is also used for the valence u, d, and s quarks. The heavy valence bottom (b) quarks use the Sheikholeslami-Wohlert (SW) Wilson-clover action [47] with the Fermilab interpretation [48].
Some of the parameters used to generate the configurations are listed in Table I [16] and [18]. The tadpole factor u 0 appearing in the one-loop improved Lüscher-Weisz gauge action and in the asqtad fermion action are determined from the fourth root of the average plaquette.
The parameters used in the valence quarks and in generating correlation functions are listed in Table II. The valence light quarks are degenerate with the sea quarks, i.e., am l = am l ; the valence s quark masses are set to our best determination of the s quark mass on each ensemble, based on all of our analysis of the asqtad ensembles. In general am h < am h .
The heavy b quark Wilson fermions with SW lattice action are controlled by the hopping parameter κ and the clover coefficient of the SW action c sw . We use κ b to denote the values 1 In this paper, we use primed quantities to denote the sea quarks and the unprimed for the valence quark. the sea-quark mass ratios am l /am h , the gauge coupling β, the tadpole improvement factor u 0 , the number of gauge-field configurations N conf , and the pion mass times the box linear spatial size M π L (L = N s a). The gauge-field configurations can be downloaded using the DOI links provided in Refs. [49][50][51][52][53][54][55][56][57][58]. used in the computation. We use the tadpole-improved tree-level value for c sw = u −3 0 , with u 0 listed in Table I. The parameter d 1 is used for the correlation function generation and will be explained later in Sec. III C.  [14] for the SW action are listed in the fourth column, where the first error is the statistics plus fitting error and the second one is due to the uncertainty in the lattice spacing. The critical κ crit for the SW action are listed in the fifth column. The errors of κ crit are in the last digit. We also list the Goldstone pion mass (M π ) and root-mean-square (RMS)  Table III lists the parameters derived from the lattice simulation. The relative lattice scale is set by calculating r 1 /a on each ensemble, where r 1 is related to the force between static quarks, r 2 1 F (r 1 ) = 1.0 [59,60]. A mass-independent procedure is used to set r 1 /a. We use the r 1 /a to convert all lattice quantities to r 1 units. The physical value of r 1 is determined from f π : r 1 = 0.3117 (22) fm [22,61]. The physical value κ b [14], corresponding to the physical b-quark mass, and the critical value κ crit , corresponding to the zero quark masses in the SW action on each ensemble, are also listed in Table III. They will be used only for correcting the b-quark masses as will be discussed in Sec. IV C. The Goldstone pion mass M π and the root-mean-square (RMS) pion mass M RMS π are listed in the last two columns of Table III.

C. Interpolating operators, currents, and correlation functions
Here we specify the interpolating operators for the kaon and B s meson and the lattice vector current needed for the correlation functions in Eq. (3.8). For the kaon, the local pseudoscalar interpolating operator is used where χ(t, x) is the one-component staggered fermion field.
The B s meson interpolating operator contains a b-quark field, simulated with the improved Wilson action, and a light staggered field for the s-quark [11,62,63] O ) where ψ(t, y) is the four-component b-quark field, and S(x, y) is a spatial smearing function.
We use two smearing functions for the B s meson. One is the local S(x, y) = δ(x − y). The other one is the ground-state 1S wave function of the Richardson potential [61].
The lattice vector current operator in Eqs. (3.7) and (3.8c) is defined as in Refs. [11] and [62] where the rotated b-quark field Ψ, defined by removes O(a) discretization effects from the current [48]. Here D lat is a symmetric nearestneighbor covariant difference operator. The coefficient d 1 , shown in Table II, is set to its tadpole-improved tree-level value so that the lattice vector current is tree-level O(a) improved.
The renormalization constant Z Vµ , needed to match the lattice vector current to its continuum counterpart (see Eq. (3.7)), is determined using a mostly nonperturbative renormalization procedure [64,65]: where Z V 4 bb and Z V 4 ll are the renormalization factors for the flavor-diagonal band lightquark temporal vector currents that are calculated nonperturbatively in Ref. [16] and listed in Table IV. The remaining flavor-off-diagonal parameters rho V µ are calculated to one-loop order in perturbation theory, separately from this analysis, and also listed in Table IV. In  order to reduce subjectivity in our analysis, we employed a blinding procedure in the form of a small multiplicative offset applied to the ρ factors and known to only two of the authors.
This blinding factor was subsequently disclosed and removed only after the analysis choices were finalized.

IV. ANALYSIS
With lattice correlation functions in hand, we follow the steps outlined near the end of Sec. III A to determine the form factors defined there, where we make use of the spectral decomposition of the correlation functions to extract the desired parameters. The two-and three-point functions take the form [62]: (4.2c) The (−1) n(t+1) and (−1) n(T −t) terms in Eq. As listed in Table V, there are 4 or 8 time sources for each ensemble. The two-point correlation functions are averaged together and folded around N t /2 before constructing the ensemble-averaged propagators and covariance matrix required for the two-point function fits. We use Bayesian constraints with Gaussian priors to perform fits to the correlation functions which include excited states. We vary the number of states and range of time slices included in the fits to separate excited state contributions from the desired ground state parameters and obtain reliable estimates of the uncertainties. The fit ranges are generally determined according to the following rules: t max is the largest value of t where the fractional error in the correlation function is smaller than 3%; t min is chosen small enough to get a good handle on the excited states and to obtain a good correlated p value as defined in Ref. [66]. The fit ranges [t min , t max ] for different lattice spacings are also adjusted so that the physical distances are similar. Our fit functions include the same number of opposite parity states as regular parity states. The number-of-states parameter N in Eq.  Here C 2 stands for the lattice two-point correlation function for the kaon or B s meson. The Bs , Z   The left panel of Fig. 3 shows an example of the stability plot for the B s meson. Fit intervals are chosen based on these plots and are listed in Table VI. Representative fit results for the kaon are shown in the right panel of Fig. 3.
The fit results for the kaon energies and overlap factors can be compared with the continuum relations to study momentum-dependent discretization errors. As illustrated in Fig. 4 we find that the E K and Z   (t, T ) defined in Refs. [16] and [18]. The prior widths for the above parameters are chosen to be 0.1 or larger in lattice units. The priors for all the other D µ mn are chosen to be 0.1 ± 2.0. The ground-state energies obtained from the combined two-and three-point correlator fits are consistent with those from the two-point fits as described in Sec. IV A. Figure 5 shows that the fitted f lat coming from the combined fit is in slight tension with the constructed ratioR Bs→K 3,0 (t, T ) defined in Refs. [16] and [18]. This small difference comes from excited state contributions still present in the ratio but accounted for in the fit method   [16] and [18]. The black curve is the ratio constructed directly from combined fit results. The small difference between the green band and the ratio comes from excited state contributions still present in the ratio but accounted for in the fit method used here. The three-point correlator fit minimum are fixed to be t 3pt min = 3 and the maximum is varied between 12 and 17. The preferred fit ranges are shown with filled points. used here. We find that they are significant at the present level of precision. Figure 6 shows an example of the stability of the fit result when varying the fit range.
In summary, the form factors f and f ⊥ are obtained from D µ 00 and E K according to Eqs. (3.6) and (4.2c), after adding the renormalization factors as in Eq. (3.6).

C. Heavy bottom quark mass correction
The heavy valence b quark is simulated with the Sheikholeslami-Wohlert (SW) action [47] with the Fermilab interpretation [48]. The b-quark mass is controlled by the hopping parameter κ b . The hopping parameter κ b used in the simulations differs slightly from the physical value κ b as can be seen in Tables II and III. We need to correct the form factors to account for these small shifts. A detailed description of the κ tuning analysis and results is provided in Appendix C of Ref. [14]. We use the method described in Refs. [16] and [18] to adjust the form factors to account for the slightly mistuned values of κ b . The relative change in the form factors under small variations of the b-quark mass can be described as where m 2 is the physical b-quark kinetic mass, m 2 is the b-quark mass used in the production run. The slopes ∂ ln f ∂ ln m 2 were determined in Ref. [16]. The corrections to the form factors are about 0.1-1.8% on different ensembles.

D. Chiral-continuum extrapolation
The lattice form factors extracted from the correlation functions as described in Sec. IV B are obtained at the three finite lattice spacings and unphysical light-quark masses listed in Table I. Here we extrapolate them to the continuum limit and physical light-quark masses using SU(2) hard-kaon heavy-meson rooted staggered chiral perturbation theory (HMrSχPT) [68,69]. Based on previous experience with the analyses of similar processes in Refs. [16] and [18], this best describes the data. Heavy-quark discretization effects are also taken into account in the chiral-continuum extrapolation.
We employ the HMrSχPT expansion at next-to-leading order (NLO) in SU(2), leading order in 1/M B , where M B is the B-meson mass, and include next-to-next-to-leading-order (NNLO) analytic and generic discretization terms. In the SU(2) hard-kaon limit, the valence and sea s-quark masses are taken to be infinitely heavy and hence dropped from the HMrSχPT formula; the large kaon energy is integrated out, and its effects are absorbed into the low-energy constants (LECs). In Ref. [18], the conversion rules for B → K and B → π processes from SU(3) HMrSχPT to SU(2) hard-kaon and hard-pion limits were derived.
Here we follow the same procedure to obtain the corresponding formula for B s → K ν. The details are presented in Appendix A.
The NLO expression for B s → K ν form factors in the SU(2) hard-kaon limit that we obtain is where P = or ⊥, δf SU (2) P,logs are the non-analytic contributions from the light-quark mass and lattice spacing, and the χ variables are dimensionless. They are defined in Eqs. (A3) and (A7). The leading-order factor is where f π is the decay constant involved, and g π is the B * Bπ coupling constant 2 . The ∆ * P term takes the pole contribution into account and is determined by requiring f and f ⊥ to have the same poles as the physical form factors f 0 and f + , respectively. This is reasonable because, by Eq. (3.5), f is dominated by the 0 + contributions of f 0 , and f ⊥ , by the 1 − contributions of f + in the q 2 range considered. Using Eq. (3.2), one obtains the exact expression for ∆ * P : The vector meson (with J P = 1 − ) has been experimentally measured [5] to be M B * = 5324.65 (25) MeV; the scalar B * meson (with J P = 0 + ) has not been observed experimentally, but a lattice calculation [70] estimates the mass difference between 0 + and 0 − states to be around 400 MeV: The vector-meson mass M B * is below the Bπ production threshold that is involved in the NNLO analytic terms are included in the fits to take into account higher-order contributions. The leading heavy b-quark discretization effects are also included. The expressions represents c 0 P , · · · , c 5 P as shown in Eq. (4.6). c NNLO p represents c 6 P , · · · , c 14 P as shown in Eq. (4.10b). h P represents h 1 P , · · · , h 6 P as appears in Eq. (4.10c). for the NNLO fit functions are Here we discuss tests of the robustness of this error estimate to check for the presence of residual truncation effects. We also consider other sources of error not already included in our chiral-continuum fit function and construct a complete systematic error budget over the range of q 2 for which we have lattice data, 19 GeV 2 q 2 24 GeV 2 .

A. Chiral-continuum extrapolation errors
Our central fit uses the NNLO SU(2) hard kaon HMrSχPT fit function described in Eq. (4.10). In order to study truncation effects, we consider variations of the central fit function. We also perform fits which include fewer form factor data.
We estimate chiral truncation effects by comparing our central NNLO fit with fits using either only the NLO function, defined in Eq. (4.6), or a fit function that includes the complete set of next-to-NNLO (NNNLO) terms. The coefficients of the NNNLO terms are constrained with the same priors as the NNLO ones. Figure 8 shows the comparison of results for the f + form factor from the three fits. The corresponding results for f 0 are similar. We see that the results from the three different fits are consistent with each other over the range of q 2 where the simulation data are located. The NNNLO errors at small q 2 are larger since the data points in that region are scarce as can be seen in Fig. 7, and the fits cannot determine the higher order terms accurately. The truncation errors are well saturated in the q 2 19 GeV 2 region, and therefore it is unnecessary to add an additional systematic error.
The SU(2) hard-kaon formula is used for the central fit. To see how other choices of the HMrSχPT formula affect the fits, we performed the fit with soft-kaon HMrSχPT. The resulting difference is small, especially for the f + form factor. This can be seen in Fig. 9.
Since the valence s-quark masses are not equal to the sea ones, the corresponding SU(3) HMrSχPT formula are extremely complicated. We therefore did not perform any trial fits with the SU(3) formula. Nevertheless, from previous experience, [16,18], SU(3) HMrSχPT typically does not provide a good description of the data.
Our results with kaon momentum up to 2π(1, 1, 1)/N s are used in the central chiral fit.
To check how the kaon energy range affects the results, we perform the fit omitting the p K = 2π(1, 1, 1)/N s data. The differences are shown in Fig. 9. Again, the difference is small especially for the f + form factor and also for f 0 at q 2 19 GeV 2 .
Based on the tests discussed above and visually summarized in Figs. 8 and 9, we find that the deviations between the results from the central fit and the alternative fits are smaller than the statistical error of the preferred central fit. We therefore do not assign additional systematic errors due to these sources.

B. Current renormalization uncertainties
The mostly nonperturbative renormalization procedure, described in Eq.  For ρ V 4 and ρ V i , the dominant source of error is the truncation at one-loop order in perturbation theory. As seen in Table IV, the one-loop corrections provided by ρ V 4 and ρ V i are small, with ρ V 4 (ρ V i ) deviating from unity by less than 1% (2.4%). Here, we adopt the estimate of the perturbative truncation error presented in Ref. [16], which yields an uncertainty of 1% on both ρ V 4 and ρ V i . This estimate is consistent with the observed differences between nonperturbative [71] and perturbative [72] calculations of ρ V 4 = ρ A 4 , discussed in Ref. [71]. In particular, the observed differences decrease in the continuum limit, as expected. Note, however, that the nonperturbative result [71] employs the HISQ action for the light quarks, while our one-loop results [72] employ the asqtad action. For this reason, the comparison is suggestive but not definitive.

C. Lattice-scale uncertainties
The dimensionful form factors f ⊥ and f , and meson energies and masses are converted to physical units via the relative scales r 1 /a listed in Table III and the absolute scale r 1 = 0.3117 (22) fm [61]. The statistical errors on r 1 /a are small and their effects on the form factors can be neglected. We estimate the error due to the uncertainty of r 1 , as before, by shifting its value by one standard deviation and repeating the chiral fit. The shifts on the form factors f +,0 are at most 0.8% in the range of simulated momenta.

D. Quark mass uncertainties
The continuum physical form factors are obtained by evaluating the chiral-continuum extrapolated functions, as discussed in Sec. IV D at the physical averaged uand d-quark masses, namely r 1 m phys ud = 0.000965 (33), and the physical s-quark mass r 1 m phys s = 0.0265 (8) as determined by analyzing the light pseudoscalar meson spectrum [22]. The error due to the uncertainties in these masses is obtained by varying their central values by one standard deviation to find the corresponding changes in the form factors. The maximum changes are below 0.15% in the simulated q 2 region.

E. Uncertainties arising from the bottom quark mass correction
As explained in Sec. IV C, the form factors are adjusted to account for the slightly mistuned valence b-quark masses before the chiral-continuum extrapolation. This accounts for the dominant effect from b-quark mass mistuning. The errors on the form factors due to the uncertainties in the κ b -correction factors and the tuned κ b values are taken into account by following the procedure described in Ref. [16]. A q 2 -independent 0.4% error due to tuning κ b is assigned to both f + and f 0 .

F. Finite volume effects
Finite-volume effects, estimated by comparing infinite-volume integrals with finite sums in HMrSχPT, are negligibly small [16,18], so they are omitted from the total error budget.

G. Summary of the statistical and systematic error budgets
The systematic errors discussed in this section are summarized in Fig. 10. We see that the largest source of systematic uncertainty by far comes from the chiral-continuum extrapolation, which includes higher-order discretization effects. This is especially obvious at small q 2 , i.e., large r 1 E K , because the statistical errors of the correlations functions increase with increasing recoil momentum so that the corresponding form factors at large r 1 E K have large errors. This is also due to a lack of data points in the large r 1 E K region as shown in Fig. 7.
Furthermore, the HMrSχPT used to perform the extrapolation is valid only for moderate errors(%) E K . This is a generic feature common to all similar lattice calculations. Our aim, however, is to get the form factor in the whole kinematically allowed region, all the way to q 2 = 0. In the next section, Sec. VI, we will describe how the extrapolation can be done by including physical information to control the error in the small q 2 region.
The sub-dominant errors, excluding the chiral-continuum extrapolation error, have mild q 2 dependence. Following Ref. [16] we therefore treat them as constants in q 2 when propagating them to the z-parametrization fit in Sec. VI B. We conservatively take the maximum estimated error from each source in the simulated q 2 range and add them in quadrature.
Specifically, the overall additional systematic error is 1.4% for both f + and f 0 , which is added to the covariance function of the chiral-continuum fit using the procedure described in Ref. [16] prior to the next step in the analysis described in the following section.

VI. CONTINUUM FORM FACTORS
The continuum form factors obtained from the chiral-continuum extrapolations described in the previous two sections are reliable only in the high momentum transfer q 2 17 GeV 2 region. In this section, we use a model-independent parametrization and expansion, namely the z-parametrization, to extrapolate the form factors to the whole kinematically allowed region. This parametrization and expansion is based on the analyticity of the form factors and angular momentum conservation. The parametrization we used was introduced by Bourrely, Caprini, and Lellouch (BCL) [73] and the fitting procedure and extrapolation technique was first introduced in our previous B → π ν paper [16].
In Sec. VI A, we briefly review the z-parametrization and give the expansion form used in the analysis. In Sec. VI B, we present the extrapolated continuum form factors in the whole kinematically allowed region. The results are shown in Table X, and Figs. 12 and 13.
A comparison with results of other groups is presented in Sec. VI C.

A. z parametrization of form factors
Before discussing the details of the method, let us first consider the properties of the semileptonic form factors. Causality and unitarity [74] imply that the B s → K ν semileptonic form factors are real analytic functions 3 in the complex q 2 -plane with a cut from q 2 > t cut to ∞, except at physical poles below t cut . The parameter t cut is the particle-pair- The pole for the vector form factor is below the cut; while the one for the scalar form factor is above it. The above-threshold pole corresponds to an unstable particle, or resonance, and may appear only on the second Riemann sheet.
From deep-inelastic-scattering experiments and perturbative QCD scaling [75,76], it is known that the semileptonic form factors vanish rapidly as 1/q 2 , up to logarithmic corrections, when q 2 approaches minus infinity.
Near the threshold t cut , the form factors have the following scaling behavior with l = 0 for f 0 and l = 1 for f + , obtained from simple partial wave analysis. Now let us look at the z parametrization. The z parametrization involves a conformal mapping. Conventionally, the variable q 2 is mapped to a new variable z according to where t 0 is a parameter that can be chosen to optimize the mapping. The maximum momentum transfer allowed in the semileptonic B s → K ν decay is defined as for convenience. This conformal mapping was first considered in Ref. [77] and further developed and used to get model-independent constraints, usually called "unitarity bounds", on form factors in Ref. [78]. A stronger constraint based on heavy-quark power counting was derived in Ref. [79]. The conformal transformation Eq. (6.3) maps the physical semileptonic region 0 ≤ q 2 ≤ t − onto a small region on the real z axis, the upper edge of the cut onto the upper edge of the unit circle, the lower edge of the cut onto the lower edge of the unit circle, the limiting points q 2 = ±∞ to z = 1, and q 2 = t cut to z = −1. The complex q 2 cut plane is mapped onto the unit disk in the z plane with the cut mapping onto the unit circle.
The parameter t 0 can be chosen such that the semileptonic region is centered around z = 0 after the conformal mapping. This is obtained by solving the equation The solution for t 0 is This mapping is schematically shown in Fig. 11 with small lepton masses ignored and with the optimized t 0 as defined in Eq. (6.6). Under the above transformation, the form factors are always in the region where |z| < 1, and therefore they can be parametrized as a power series in z. Since the physical semileptonic region in terms of z is usually small, |z| ≤ 0.205 for B s → K ν, this parametrization converges quickly. Table VIII has a list of quantities in terms of r 1 E K , q 2 , and z parameters.
Two commonly used parametrizations are given by Boyd, Grinstein and Lebed (BGL) [80] and by Bourrely, Caprini and Lellouch (BCL) [73]. Here we use the BCL parametrization as given by The factors 1/(1 − q 2 /m 2 B * ) take the poles into account and ensure the asymptotic scaling, f (q 2 ) ∼ 1/q 2 at large q 2 . Moreover, the scaling condition of Eq. (6.2) near t cut is also enforced for f + . Note that Eq. (6.2) in the q 2 -plane imply the following relation The form factors constructed with this BCL z parametrization satisfy all three properties of the semileptonic form factors discussed at the beginning of this section.   Table IX lists the relevant meson masses used in the z-parametrization fit.
The functional method introduced in Ref. [16] is used to perform the z-parametrization fit, where, following Ref. [16], we take as inputs the results from the chiral-continuum extrapolation and systematic error analysis as presented in Sec. V G).
Our preferred (central) fit has K = 4, where K is the number of terms in the expansion in Eq. (6.7). The results of this fit are shown in Table X further improve the f + form-factor fit. The fit parameters also satisfy the unitarity condition [73] and the condition estimated from heavy-quark power counting [79]. Adding the heavy-quark constraint does not affect the fit results. The kinematic constraint is enforced by requiring f + and f 0 to be exactly equal at the q 2 = 0 point. In practice, we set a prior in the z-parametrization fit f + (q 2 = 0) − f 0 (q 2 = 0) = 0, (6.9)  Table VIII. The meson poles are listed in Table IX. with width = 10 −10 . When further increasing the expansion order to K = 4, the central value of the form factors at q 2 = 0 agrees with the results with K = 3, but the error increases. The unitarity and heavy-quark constraints are still satisfied automatically. The results stabilize at K = 4 and do not change with K = 5. We conclude that the K = 4 fit with the kinematic constraint includes the systematic uncertainty due to truncating the z-parametrization series.
The left panel of Fig. 12 shows the preferred K = 4 form-factor results, with poles removed, as functions of z. The q 2 = 0 point is at the right end of the plot. Note that the shape of the form factors as functions of z is parametrization dependent. For convenience, the right panel of Fig. 12 shows the form factors as functions of q 2 . The q 2 dependence of the form factors is parametrization independent and can be used directly to compare with results of other groups.

C. Comparison with existing results
Several other groups have also calculated the same form factors. We note that Refs. [30] and [31] use the B s K threshold instead of Bπ in their implementation of the z parametrization. Since the z-parameter, by definition (see Eq. (6.3)), depends on the threshold (t cut ), TABLE X. The results of the preferred z-parametrization fit from Eqs. (6.7, 6.6, 6.9) and Table IX with K = 4. These values can be used to reconstruct the form factors as explained in Appendix B.
The correlation matrix is listed with only four digits after the decimal point. The correlation matrix has one near zero eigenvalue due to the kinematic constraint used. See Appendix B for details.
Correlation matrix we cannot directly compare the z-dependence of our form factors with those of Refs. [30] and [31]. We therefore compare our form factors with those from other lattice QCD calculations only as functions of q 2 . This is shown in Figure 13.
The results of the HPQCD Collaboration [30] are based on (2+1)-flavor-MILC-asqtad configurations for the sea quarks, and employ the HISQ action for the light valence quarks, and lattice NRQCD for the heavy b-quark. The RBC and UKQCD Collaborations [31] use (2+1)-flavor-domain-wall fermions for the sea quarks and light valence quarks, and a variant [81,82] of the Fermilab action for the heavy b-quark. The ALPHA Collaboration [32] uses leading-order lattice HQET to get the form factors at one point, q 2 = 22.12 GeV 2 .
While our results are consistent with those from Refs. [31] and [32], they are in tension with HPQCD's results [30]. We note that Ref. [30] employs the so-called modified z-expansion, where the chiral-continuum extrapolation is combined with the z-expansion into one fit function by modifying the z-coefficients with lattice-spacing and light-quark-mass dependent terms. This procedure may affect the shape of the form factors. Indeed, in their calculation of the form factors for the B → K + − decay in Ref. [83], the HPQCD Collaboration compared the form factors obtained after the modified z-expansion with the results from a two-step method that is very similar to ours, performing first a chiral-continuum extrapolation, and then a z-expansion fit. While they find only small differences between the two sets of form factors, those obtained from their implementation of the two-step method are in better agreement with the results of Ref. [18]. However, unlike the case at hand, the form factors of Ref. [18] are not in significant tension with HPQCD's results of Ref. [83]. We see that the tension between our B s → K ν form factor results and those of Ref. [30] increases with decreasing q 2 to roughly 2.3σ at q 2 = 0. The RBC and UKQCD Collaborations [31], on the other hand, adopt the same procedure as we do, namely a chiral-continuum extrapolation at high q 2 , followed by a z-expansion extrapolation to q 2 = 0.
A comparison of the form factor at q 2 = 0 is shown in Fig. 14, where we also include results from calculations using light-cone sum rules [34,35], a relativistic quark model [33], and NLO perturbative QCD [36].

VII. PHENOMENOLOGICAL APPLICATIONS
The angular-dependent differential decay rate for B s → K ν is given in Eq. (2.3). One can construct at most three independent observables from there. In the following, we will consider the differential decay rate dΓ/dq 2 in Sec. VII A, the forward-backward asymmetry A F B (q 2 ) in Sec. VII B, and the lepton polarization asymmetry A pol (q 2 ) in Sec. VII C. The latter two quantities are sensitive to the mass of the final-state charged lepton. In Sec. VII D, we also construct the ratios of the scalar and vector form factors between the B s → K ν and B s → D s ν decays.
A. Decay rate The differential decay rate can be obtained from Eq. (2.3) by integrating over the angle θ , which yields In Fig. 15, we plot the Standard Model predictions of the differential decay rate divided by |V ub | 2 over the whole kinematic range of q 2 for B s → Kµν and B s → Kτ ν.
One can also explore the ratio of the differential decay rates Figure 16 shows the prediction for R τ /µ (q 2 ).
The total decay rate is given by  In Appendix C, we also provide partially integrated differential decay rates in evenly  spaced q 2 bins.
The ratio of the total decay rate is which takes the correlations between the form factors into account and is more precise than directly using Eq. (7.4).

B. Forward-backward asymmetry
The forward-backward asymmetry, A F B , which depends on the linear cos θ term in Eq. (2.3), is given by The Standard Model predictions for the forward-backward asymmetry divided by |V ub | 2 are shown in Fig. 17. For the corresponding integrated quantities we find The normalized forward-backward asymmetry is given bȳ The normalized lepton polarization asymmetry is defined as from the differential decay rates with definite lepton helicity [84] dΓ (7.11b) Here the superscripts + (−) imply a right-(left-)handed lepton in the final state. The lepton is produced via the V − A current in the Standard Model, and therefore the electron and muon are mainly left-handed polarized. The A µ pol is close to one in the whole q 2 range. Here we provide the normalized lepton polarization asymmetry A µ pol and A τ pol as functions of q 2 in Fig. 18.   We also calculate the ratios of the scalar and vector form factors between the B s → K ν and B s → D s ν semileptonic decays. The ratios can be used along with future experimental results to determine the ratio of the CKM matrix elements |V ub /V cb |.
First, we reconstruct the B s → D s ν form factors from our previous papers [13,15].  [13] and [15], respectively. They are shown in Fig. 19. The

VIII. SUMMARY AND OUTLOOK
Using six strategically selected ensembles of MILC asqtad 2+1 flavor gauge configurations, we have calculated the form factors f + (q 2 ) and f 0 (q 2 ) needed to understand the semileptonic decay B s → K ν. We present predictions of the differential decay rate (divided by |V ub | 2 ) for both light (e or µ) or heavy (τ ) final-state leptons. Once the experimental data become  construct the form factor ratios as shown in Fig. 21 and Fig. 22 ical interest include the forward-backward asymmetry A FB (q 2 ) and the lepton polarization asymmetry A pol (q 2 ). We also present ratios of the form factors f + and f 0 for B s → K ν and B s → D s ν as functions of both q 2 and w. These may be valuable for determining |V ub /V cb |.
Although there are no published results for the decay B s → K ν, this process is under investigation by the LHCb experiment, and will be studied by the Belle II Collaboration when they run at the Υ(5S) resonance, which is a copious source of B s andB s mesons.
On the theoretical side, we have plans to reduce the contributions from the dominant sources of systematic errors in upcoming calculations, which include chiral extrapolation, light and heavy-quark discretization, and renormalization. The gauge ensembles generated by the MILC collaboration with four flavors of HISQ sea quarks [28,29] are a crucial ingredient in these plans. These ensembles cover a lattice spacing range of approximately 0.15-0.045 fm with physical light quark masses and a dynamical charm quark. The chiral extrapolation becomes a chiral interpolation, and the reduced taste breaking of the HISQ action greatly reduces light quark discretization errors. Using these ensembles, we will be taking two approaches to the b quark. First, we have started a project using Fermilab b quarks (as in this project) and HISQ light valence quarks. Preliminary results were already reported in Refs. [26] and [86]. As a further small improvement compared to this work, we will include the full correlation matrix between form factors for different processes in our final results. In our second approach, we plan to use the HISQ formalism for the b quark to calculate semileptonic B (s) -and D-meson decay form factors again on the HISQ ensembles.
Heavy-quark discretization errors are simpler with the HISQ action than with the Fermilab approach, and can be controlled with high precision by including ensembles with very fine lattice spacings in the range of a ≈ 0.03 − 0.042 fm. The heavy-HISQ approach also allows us to take advantage of Ward identities when renormalizing the currents. Indeed, our recent work [87] where the subscript P stands for or ⊥; c i P are coefficients and the corresponding rescaled quantities in Eq. (A7) will be determined by the chiral fits; δf P,logs contains the one-loop nonanalytic contributions and wave-function renormalizations; m x and m y are the corresponding valence quark masses; m u , m d and m s are sea quark masses; E = p · v is the P xy meson energy in the B x meson rest frame; and a is the lattice spacing. The leading order terms for f and f ⊥ are where f is the decay constant involved; g is the coupling constant; ∆ * xy,P is the mass difference between the quantum number J P = 0 + or 1 − B * y meson and the pseudoscalar B x meson masses at leading order in the chiral expansion, i.e., ∆ * xy = B * y − B x ; and D logs is the nonanalytic self-energy contribution. The scalar pole was not included in Ref. [69] as the 0 + meson is not in leading order HMrSχPT. It is added here phenomenologically as explained in Sec. IV D.
For the B s → K ν analysis considered here, x = s and y = u = d = u = d. Here we use primed quantities to denote the valence quarks and the unprimed for the sea quarks 5 .
All the data generated for B s → K ν analysis are partially quenched points, i.e., m s = m s .
In the SU(2) limit, the s-quark mass is treated as infinitely heavy and all the explicit m s dependent terms are removed from the formula. However, one will still need to keep the mass difference, m s − m s , in the leading-order analytic term to take the partial quenching effects into account. In the hard-kaon limit, the kaon with a large energy, E, is integrated out in the nonanalytic chiral expressions. These two limits greatly simplify the expressions of the chiral logs. Following the recipes presented in the Appendix of Ref. [18], we obtain  We can regroup relevant terms in Eq. (A1), drop the m s dependent term due to the SU (2) limit, and write the formula as the following f P,NLO =f (0) P [c 0 P (1 + δf P,logs ) + (c u P + 2c sea P ) 3 3m u + c sea P (m s − m s ) + c E P E + c E 2 P E 2 + c a 2 P a 2 ].
We can further write all the expansion parameters in terms of dimensionless ones where µ is the leading-order low-energy constant that relates the tree-level mass of a taste-ξ meson composed of quarks of flavor i and j to the corresponding quark masses Here ∆ ξ is the staggered fermion taste splitting. The numerical values of µ and ∆ ξ are determined by the MILC Collaboration and are shown in Table XI   To get the q 2 dependence of the form factors as in the right panel of Fig. 12, one needs the relation between z and q 2 . In this paper, the mapping is defined in Eqs. (6.3), (6.1), (6.4), and (6.6). One can then solve Eq. (6.3) to get q 2 in terms of z: Once we have the form factors as functions of z from Appendix B 1, we can then use Eq. (B1) to change the variable to get the q 2 dependence.

Dealing with the near zero eigenvalue in the covariance matrix
In Table X the fit parameter standard deviations and the correlation matrix are listed.
To get the covariance matrix, one only needs to follow the usual procedure to rescale the correlation matrix. The following is the detailed procedure.
Suppose the standard deviation of the fit parameters is Σ = [σ 1 , σ 2 , · · · , σ n ] , and the matrix D is a diagonal matrix with diagonal elements Σ. The correlation matrix is denoted as R and the covariance matrix is denoted as S. The relations among D, R, and S are Alternatively, one can use the following relation to directly convert the matrix elements where there is no summation over the repeated indices. The covariance matrix, or the inverse of it, is useful when combining form factor results from different sources. It is difficult to calculate the inverse of the covariance matrix from the results listed in Table X.
This is because we imposed the kinematical constraint Eq. (6.9) with = 10 −10 in the zparametrization fit. This results in a near zero eigenvalue in the covariance matrix. The kinematic constraint is equivalent to reducing one parameter in the z parametrization. In principle, one can first reduce one parameter, say b 0 3 , in Eq. (6.7), express it in terms of the other b +,0 TABLE XII. The binned differential decay rates, defined in Eq. (C1), and their correlations for