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

TL;DR

This paper calculates the form factors for the semileptonic decay $B_s o K\,\ell\nu$ using lattice QCD with multiple ensembles, providing theoretical predictions crucial for determining the CKM matrix element $|V_{ub}|$.

Contribution

The study presents the first lattice QCD calculation of $B_s\to K$ form factors across the full kinematic range with multiple lattice spacings and chiral extrapolation techniques.

Findings

Predicted differential decay rates for $B_s\to K\mu\nu$ and $B_s\to K\tau\nu$.

Provided form factor ratios and asymmetries for $B_s\to K\ell\nu$ decays.

Results can be combined with experimental data to extract $|V_{ub}|$.

Abstract

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}|$.

Figures

Tables

Table 1. Table 1: Parameters used for generating the lattice QCD gauge fields. The columns from left to right are approximate lattice spacing a 𝑎 a in fm, the lattice dimensions in lattice units N s 3 × N t superscript subscript 𝑁 𝑠 3 subscript 𝑁 𝑡 N_{s}^{3}\times N_{t} , the sea-quark mass ratios a ​ m l ′ / a ​ m h ′ 𝑎 superscript subscript 𝑚 𝑙 ′ 𝑎 superscript subscript 𝑚 ℎ ′ am_{l}^{\prime}/am_{h}^{\prime} , the gauge coupling β 𝛽 \beta , the tadpole improvement factor u 0 subscript 𝑢 0 u_{0} , the number of gauge-field configurations N conf subscript 𝑁 conf N_{\text{conf}} , and the pion mass times the box linear spatial size M π ​ L subscript 𝑀 𝜋 𝐿 M_{\pi}L ( L = N s ​ a 𝐿 subscript 𝑁 𝑠 𝑎 L=N_{s}a ). The gauge-field configurations can be downloaded using the DOI links provided in Refs.  MILC Collaboration ( 2015a , b , c , d , e , f , g , h , i , j ) .

$\approx$a (fm)	$N_s^3\times N_t$	$a{m}_l^{\prime }/a{m}_h^{\prime }$	$\beta$	$u_0$	$N_{\text{conf}}$	$M_{\pi }L$
0.12 MILC Collaboration (2015a)	${24}^3\times 64$	0.0050/0.050	6.76	0.8678	2099	3.8
0.09 MILC Collaboration (2015b, c, d)	${28}^3\times 96$	0.0062/0.031	7.09	0.8782	1931	4.1
0.09 MILC Collaboration (2015e)	${32}^3\times 96$	0.00465/0.031	7.085	0.8781	1015	4.1
0.09 MILC Collaboration (2015f, g)	${40}^3\times 96$	0.0031/0.031	7.08	0.8779	1015	4.2
0.09 MILC Collaboration (2015h)	${64}^3\times 96$	0.00155/0.031	7.075	0.877805	791	4.8
0.06 MILC Collaboration (2015i, j)	${64}^3\times 144$	0.0018/0.018	7.46	0.88764	827	4.3

Table 2. Table 2: Parameters used for generating the valence quarks. The approximate lattice spacing a 𝑎 a and lattice dimensions N s 3 × N t superscript subscript 𝑁 𝑠 3 subscript 𝑁 𝑡 N_{s}^{3}\times N_{t} in the first two column identify the ensmeble. The light valence quarks m l subscript 𝑚 𝑙 m_{l} are degenerate with the sea quarks m l ′ superscript subscript 𝑚 𝑙 ′ m_{l}^{\prime} . The valence s 𝑠 s quarks m h subscript 𝑚 ℎ m_{h} are better tuned than the sea s 𝑠 s quarks m h ′ superscript subscript 𝑚 ℎ ′ m_{h}^{\prime} . The parameters c sw subscript 𝑐 sw c_{\mathrm{sw}} and κ b ′ superscript subscript 𝜅 𝑏 ′ \kappa_{b}^{\prime} are used in the SW action for b 𝑏 b quarks. The rotation parameter d 1 subscript 𝑑 1 d_{1} is used in the current.

$\approx$a (fm)	$N_s^3\times N_t$	a$m_l$/a$m_h$	$c_{\mathrm{sw}}$	${\kappa }_b^{\prime }$	$d_1$
0.12	${24}^3\times 64$	0.0050/0.0336	1.53	0.0901	0.09332
0.09	${28}^3\times 96$	0.0062/0.0247	1.476	0.0979	0.096765
0.09	${32}^3\times 96$	0.00465/0.0247	1.477	0.0977	0.096708
0.09	${40}^3\times 96$	0.0031/0.0247	1.478	0.0976	0.096688
0.09	${64}^3\times 96$	0.00155/0.0247	1.478	0.0976	0.0967
0.06	${64}^3\times 144$	0.0018/0.0177	1.4298	0.1052	0.0963

Table 3. Table 3: Parameters derived from the simulation. The approximate lattice spacing a 𝑎 a in fm and the lattice dimensions in lattice units N s 3 × N t superscript subscript 𝑁 𝑠 3 subscript 𝑁 𝑡 N_{s}^{3}\times N_{t} are used for identifying the ensemble. Relative scales r 1 / a subscript 𝑟 1 𝑎 r_{1}/a are listed in the third column. The statistical errors on r 1 / a subscript 𝑟 1 𝑎 r_{1}/a are 0.1 to 0.3% and the systematic errors are comparable. The physical κ b subscript 𝜅 𝑏 \kappa_{b}   Bailey  et al. ( 2014 ) 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 subscript 𝜅 crit \kappa_{\mathrm{crit}} for the SW action are listed in the fifth column. The errors of κ crit subscript 𝜅 crit \kappa_{\mathrm{crit}} are in the last digit. We also list the Goldstone pion mass ( M π subscript 𝑀 𝜋 M_{\pi} ) and root-mean-square (RMS) pion mass ( M π RMS superscript subscript 𝑀 𝜋 RMS M_{\pi}^{\textrm{RMS}} ) here.

$\approx$a (fm)	$N_s^3\times N_t$	$r_1/a$	${\kappa }_b$	${\kappa }_{\mathrm{crit}}$	$M_{\pi }$ (MeV)	$M_{\pi }^{\text{RMS}}$ (MeV)
0.12	${24}^3\times 64$	2.73859	0.0868(9)(3)	0.14096	277	456
0.09	${28}^3\times 96$	3.78873	0.0967(7)(3)	0.139119	354	413
0.09	${32}^3\times 96$	3.77163	0.0966(7)(3)	0.139134	307	374
0.09	${40}^3\times 96$	3.75459	0.0965(7)(3)	0.139173	249	329
0.09	${64}^3\times 96$	3.73761	0.0964(7)(3)	0.13919	177	277
0.06	${64}^3\times 144$	5.30734	0.1050(5)(2)	0.137678	224	255

Table 4. Table 4: Parameters for the renormalization of the form factors. The approximate lattice spacing a 𝑎 a and lattice dimensions N s 3 × N t superscript subscript 𝑁 𝑠 3 subscript 𝑁 𝑡 N_{s}^{3}\times N_{t} in the first two column identify the ensmeble. The light-light and heavy-heavy renormalization factors Z V l ​ l 4 subscript 𝑍 subscript superscript 𝑉 4 𝑙 𝑙 Z_{V^{4}_{ll}} and Z V b ​ b 4 subscript 𝑍 subscript superscript 𝑉 4 𝑏 𝑏 Z_{V^{4}_{bb}} are listed in the third and fourth columns. The one-loop estimates of ρ V i subscript 𝜌 superscript 𝑉 𝑖 \rho_{V^{i}} and ρ V 4 subscript 𝜌 superscript 𝑉 4 \rho_{V^{4}} are listed in the fifth and sixth columns. The errors shown are statistical. The complete current renormalization is obtained via Eq. ( 16 ).

$\approx$a (fm)	$N_s^3\times N_t$	$Z_{V_{ll}^4}$	$Z_{V_{bb}^4}$	${\rho }_{V^i}$	${\rho }_{V^4}$
0.12	${24}^3\times 64$	1.7410(30)	0.5015(8)	0.973082	1.006197
0.09	${28}^3\times 96$	1.7770(50)	0.4519(15)	0.975822	0.999308
0.09	${32}^3\times 96$	1.7760(50)	0.4530(15)	0.975775	0.999405
0.09	${40}^3\times 96$	1.7760(50)	0.4536(15)	0.975744	0.999441
0.09	${64}^3\times 96$	1.7760(50)	0.4536(15)	0.975703	0.999416
0.06	${64}^3\times 144$	1.8070(70)	0.4065(21)	0.979176	0.995327

Table 5. Table 5: The number of time sources N src subscript 𝑁 src N_{\mathrm{src}} used in the two- and three-point correlation function generation and the source-sink separations T sink subscript 𝑇 sink T_{\mathrm{sink}} used in the three-point correlation function generation. The approximate lattice spacing a 𝑎 a and lattice dimensions N s 3 × N t superscript subscript 𝑁 𝑠 3 subscript 𝑁 𝑡 N_{s}^{3}\times N_{t} in the first two column identify the ensmeble.

$\approx$a (fm)	$N_s^3\times N_t$	$N_{\mathrm{src}}$	$T_{\mathrm{sink}}$
0.12	${24}^3\times 64$	4	18
0.09	${28}^3\times 96$	4	25
0.09	${32}^3\times 96$	8	25
0.09	${40}^3\times 96$	8	25
0.09	${64}^3\times 96$	4	25
0.06	${64}^3\times 144$	4	36

Table 6. Table 6: Fit ranges [ t min , t max ] subscript 𝑡 min subscript 𝑡 max [t_{\text{min}},t_{\text{max}}] used in the kaon and B s subscript 𝐵 𝑠 B_{s} meson two-point correlator fits.

$\approx$a (fm)	Kaon	$B_s$ meson
0.12	[5,31]	[3,22]
0.09	[7,47]	[4,30]
0.06	[10,71]	[6,44]

Table 7. Table 7: Priors used in the chiral-continuum extrapolation fit. c p NLO superscript subscript 𝑐 𝑝 NLO c_{p}^{\mathrm{NLO}} represents c P 0 , ⋯ , c P 5 superscript subscript 𝑐 𝑃 0 ⋯ superscript subscript 𝑐 𝑃 5 c_{P}^{0},\cdots,c_{P}^{5} as shown in Eq. ( 22 ). c p NNLO superscript subscript 𝑐 𝑝 NNLO c_{p}^{\mathrm{NNLO}} represents c P 6 , ⋯ , c P 14 superscript subscript 𝑐 𝑃 6 ⋯ superscript subscript 𝑐 𝑃 14 c_{P}^{6},\cdots,c_{P}^{14} as shown in Eq. ( 26b ). h P subscript ℎ 𝑃 h_{P} represents h P 1 , ⋯ , h P 6 superscript subscript ℎ 𝑃 1 ⋯ superscript subscript ℎ 𝑃 6 h_{P}^{1},\cdots,h_{P}^{6} as appears in Eq. ( 26c ).

$f_{\pi }$	$g_{\pi }$	$c_P^{\mathrm{NLO}}$	$c_P^{\mathrm{NNLO}}$	$h_P$
130.4 MeV	0.45(8)	0(1.0)	0(0.6)	0(1.0)

Table 8. Table 8: Quantities in terms of different parameters.

	$r_1{E}_K$	$q^2\left({\mathrm{GeV}}^2\right)$	$z$
Lattice data range	$[0.846,1.71]$	$[17.4,23.3]$	$[-0.186,-0.0174]$
Physical range	$[0.780,4.28]$	$[0,23.7]$	$[-0.205,0.205]$
$t_{\_}={\left(M_{B_s}-M_K\right)}^2$	0.780	23.7	$-0.205$
$t_{\mathrm{cut}}={\left(M_B+M_{\pi }\right)}^2$	$-0.0395$	29.3	$-1.0$
$t_0=t_{\mathrm{cut}}-\sqrt{t_{\mathrm{cut}}\left(t_{\mathrm{cut}}-t_{-}\right)}$	1.84	16.5	0.0
$M_{B^{\ast }}^2\left(1^{-}\right)$	0.102	28.4	$-0.569$
$M_{B^{\ast }}^2\left(0^{+}\right)$	$-0.473$	32.3	$-0.625+0.781i$

Table 9. Table 9: Input meson masses used in the z 𝑧 z -parametrization fit.

	$M_{B_s}$	$M_K$	$M_B$	$M_{\pi }$	$M_{B^{\ast }}\left(1^{-}\right)$	$M_{B^{\ast }}\left(0^{+}\right)$
Value $\left(\mathrm{GeV}\right)$	5.36682	0.493677	5.27931	0.1349766	5.32465	5.68

Table 10. Table 10: The results of the preferred z 𝑧 z -parametrization fit from Eqs. ( 33 ,  32 ,  35 ) and Table  9 with K = 4 𝐾 4 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
	Value	$b_0^{+}$	$b_1^{+}$	$b_2^{+}$	$b_3^{+}$	$b_0^0$	$b_1^0$	$b_2^0$	$b_3^0$
$b_0^{+}$	0.3623(0.0178)	1.0000	0.6023	0.0326	-0.1288	0.7122	0.6035	0.5659	0.5516
$b_1^{+}$	-0.9559(0.1307)		1.0000	0.4735	0.2677	0.7518	0.9086	0.9009	0.8903
$b_2^{+}$	-0.8525(0.4783)			1.0000	0.9187	0.5833	0.7367	0.7340	0.7005
$b_3^{+}$	0.2785(0.6892)				1.0000	0.4355	0.5553	0.5633	0.5461
$b_0^0$	0.1981(0.0101)					1.0000	0.8667	0.7742	0.7337
$b_1^0$	-0.1661(0.1130)						1.0000	0.9687	0.9359
$b_2^0$	-0.6430(0.4385)							1.0000	0.9899
$b_3^0$	-0.3754(0.4535)								1.0000

Table 11. Table 11: Fixed parameters used in the chiral-continuum extrapolation fit function. The r 1 2 ​ a 2 ​ Δ ξ superscript subscript 𝑟 1 2 superscript 𝑎 2 subscript Δ 𝜉 r_{1}^{2}a^{2}\Delta_{\xi} with ξ = P , A , T , V , I 𝜉 𝑃 𝐴 𝑇 𝑉 𝐼 \xi=P,A,T,V,I and r 1 2 ​ a 2 ​ δ V / A ′ superscript subscript 𝑟 1 2 superscript 𝑎 2 subscript superscript 𝛿 ′ 𝑉 𝐴 r_{1}^{2}a^{2}\delta^{\prime}_{V/A} are taste splittings and hairpin parameters.

$\approx$a (fm)	0.12	0.09	0.06	0
$r_1\mu$	6.831904	6.638563	6.486649	6.015349
$r_1^2{a}^2{\Delta }_P$	0	0	0	0
$r_1^2{a}^2{\Delta }_A$	0.22705	0.07469	0.02635	0
$r_1^2{a}^2{\Delta }_T$	0.36616	0.12378	0.04298	0
$r_1^2{a}^2{\Delta }_V$	0.48026	0.15932	0.05744	0
$r_1^2{a}^2{\Delta }_I$	0.60082	0.22065	0.07039	0
$r_1^2{a}^2{\delta }_V^{\prime }$	0.0	0.0	0.0	0
$r_1^2{a}^2{\delta }_A^{\prime }$	$-0.28$	$-0.09$	$-0.03$	0

Table 12. Table 12: The binned differential decay rates, defined in Eq. ( 63 ), and their correlations for B s → K ​ μ ​ ν → subscript 𝐵 𝑠 𝐾 𝜇 𝜈 B_{s}\to K\mu\nu in twelve evenly spaced q 2 superscript 𝑞 2 q^{2} bins.

Correlation matrix
$\Delta q^2\left({\mathrm{GeV}}^2\right)$	Value	0-2	2-4	4-6	6-8	8-10	10-12	12-14	14-16	16-18	18-20	20-22	22-24
0-2	0.20707(0.14609)	$1.0000$	$0.9981$	$0.9911$	$0.9766$	$0.9510$	$0.9096$	$0.8459$	$0.7520$	$0.6171$	$0.4280$	$0.1869$	$-0.0067$
2-4	0.25741(0.14256)		$1.0000$	$0.9974$	$0.9880$	$0.9681$	$0.9332$	$0.8765$	$0.7893$	$0.6599$	$0.4715$	$0.2200$	$0.0060$
4-6	0.30678(0.13329)			$1.0000$	$0.9965$	$0.9836$	$0.9564$	$0.9082$	$0.8297$	$0.7076$	$0.5217$	$0.2600$	$0.0237$
6-8	0.35666(0.12140)				$1.0000$	$0.9952$	$0.9773$	$0.9396$	$0.8720$	$0.7599$	$0.5790$	$0.3085$	$0.0484$
8-10	0.40477(0.10755)					$1.0000$	$0.9933$	$0.9683$	$0.9147$	$0.8160$	$0.6440$	$0.3674$	$0.0830$
10-12	0.44785(0.09239)						$1.0000$	$0.9906$	$0.9545$	$0.8742$	$0.7169$	$0.4396$	$0.1323$
12-14	0.48121(0.07653)							$1.0000$	$0.9862$	$0.9306$	$0.7967$	$0.5285$	$0.2033$
14-16	0.49783(0.06056)								$1.0000$	$0.9778$	$0.8804$	$0.6382$	$0.3069$
16-18	0.48688(0.04502)									$1.0000$	$0.9585$	$0.7714$	$0.4589$
18-20	0.43098(0.03061)										$1.0000$	$0.9170$	$0.6736$
20-22	0.30246(0.01787)											$1.0000$	$0.9090$
22-24	0.08453(0.00502)												$1.0000$

Table 13. Table 13: The binned differential decay rates, defined in Eq. ( 63 ), and their correlations for B s → K ​ τ ​ ν → subscript 𝐵 𝑠 𝐾 𝜏 𝜈 B_{s}\to K\tau\nu in eleven evenly spaced q 2 superscript 𝑞 2 q^{2} bins.

Correlation matrix
$\Delta q^2\left({\mathrm{GeV}}^2\right)$	Value	2-4	4-6	6-8	8-10	10-12	12-14	14-16	16-18	18-20	20-22	22-24
2-4	0.00500(0.00269)	$1.0000$	$0.9993$	$0.9953$	$0.9851$	$0.9648$	$0.9281$	$0.8640$	$0.7522$	$0.5581$	$0.2683$	$0.0393$
4-6	0.09085(0.04127)		$1.0000$	$0.9982$	$0.9909$	$0.9740$	$0.9414$	$0.8819$	$0.7747$	$0.5837$	$0.2913$	$0.0544$
6-8	0.19913(0.07215)			$1.0000$	$0.9971$	$0.9857$	$0.9595$	$0.9074$	$0.8081$	$0.6229$	$0.3277$	$0.0797$
8-10	0.28718(0.08104)				$1.0000$	$0.9956$	$0.9779$	$0.9359$	$0.8477$	$0.6718$	$0.3761$	$0.1156$
10-12	0.36067(0.07799)					$1.0000$	$0.9931$	$0.9643$	$0.8915$	$0.7302$	$0.4384$	$0.1660$
12-14	0.42097(0.06852)						$1.0000$	$0.9885$	$0.9370$	$0.7986$	$0.5192$	$0.2379$
14-16	0.46455(0.05576)							$1.0000$	$0.9786$	$0.8763$	$0.6252$	$0.3429$
16-18	0.48267(0.04190)								$1.0000$	$0.9554$	$0.7621$	$0.4978$
18-20	0.45879(0.02887)									$1.0000$	$0.9162$	$0.7121$
20-22	0.36241(0.01781)										$1.0000$	$0.9256$
22-24	0.14148(0.00609)											$1.0000$

Equations

$$⟨K(p_K)|{\mathcal{V}}^{\mu }|B_s(p_{B_s})⟩=\left(p_K^{\mu }+p_{B_s}^{\mu }-q^{\mu }\frac{M_{B_s}^2-M_K^2}{q^2}\right)f_{+}(q^2)+q^{\mu }\frac{M_{B_s}^2-M_K^2}{q^2}{f}_0(q^2),$$

$$f_{+}(0)=f_0(0),$$

$$\frac{d^2\Gamma }{d{q}^{2}d\cos {\theta }_{\ell }}=$$

$$+\frac{4m_{\ell }^2}{q^2}(M_{B_s}^2-M_K^2)M_{B_s}|{\mathrm{p}}_K|\cos {\theta }_{\ell }\mathrm{Re}\left[f_{+}(q^2)f_0^{\ast }(q^2)\right]$$

$$+\frac{m_{\ell }^2}{q^2}(M_{B_s}^2-M_K^2{)}^2|f_0(q^2){|}^2]$$

$$v^{\mu }=p_{B_s}^{\mu }/M_{B_s}.$$

$$q^2$$

$$p_{\perp }^{\mu }\equiv p_K^{\mu }-(p_K\cdot v)v^{\mu }$$

$$⟨K(p_K)|{\mathcal{V}}^{\mu }|B_s(p_{B_s})⟩=\sqrt{2M_{B_s}}\left[v^{\mu }{f}_{\parallel }(E_K)+p_{\perp }^{\mu }{f}_{\perp }(E_K)\right].$$

$$f_{+}(q^2)$$

$$f_0(q^2)$$

$$f_{\parallel }(E_K)$$

$$f_{\perp }(E_K)$$

$${\mathcal{V}}^{\mu }(x)=Z_{V_{\mu }}{V}^{\mu }(x).$$

$$C_2^{B_s}(t;{\mathbf{p}}_{B_s}=0)$$

$$C_2^K(t;{\mathbf{p}}_K)$$

$$C_{3,\mu }^{B_s\to K}(t,T;{\mathbf{p}}_K)$$

$${\mathcal{O}}_K(t,\mathbf{x})=\bar{\chi }(t,\mathbf{x})(-1{)}^{t+x_1+x_2+x_3}\chi (t,\mathbf{x}),$$

$${\mathcal{O}}_{B_s}(t,\mathbf{x})$$

$$\Omega (t,\mathbf{x})$$

$$V^{\mu }(x)$$

$$\Psi =(1+a{d}_1\mathbf{\gamma }\cdot {\mathbf{D}}_{\mathrm{lat}})\psi ,$$

$$Z_{V_{bl}^{\mu }}={\rho }_{V^{\mu }}\sqrt{Z_{V_{bb}^4}{Z}_{V_{ll}^4}},$$

$$C_2^{B_s}(t;0)$$

$$C_2^K(t;{\mathbf{p}}_K)$$

$$C_{3,\mu }^{B_s\to K}(t,T;{\mathbf{p}}_K)$$

$$Z_{B_s}^{(n)}$$

$$Z_K^{(n)}({\mathbf{p}}_K)$$

$$D_{mn}^{\mu }$$

$$m_{\mathrm{eff}}\equiv -\log \left[C_2(t+1)/C_2(t)\right],Z_{\mathrm{eff}}^2\equiv e^{+m_{\mathrm{eff}}t}{C}_2(t).$$

$$E_K^2=M_K^2+{\mathbf{p}}_K^2,Z_K^{(0)}({\mathbf{p}}_K)=Z_K^{(0)}(0)\sqrt{\frac{M_K}{E_K}}$$

$$f(m_2^{\prime })=f(m_2)\left[1-\frac{\partial \ln f}{\partial \ln m_2}\left(\frac{m_2}{m_2^{\prime }}-1\right)\right],$$

$$f_{P,\mathrm{NLO}}=f_P^{(0)}\left[c_P^0(1+\delta f_{P,\mathrm{logs}}^{\mathrm{SU}(2)})+c_P^{\text{1}}{\chi }_{\text{l}}+c_P^{\text{2}}{\chi }_{\text{h}}+c_P^3{\chi }_E+c_P^4{\chi }_E^2+c_P^5{\chi }_a^2\right],$$

$$f_P^{(0)}=\frac{1}{f_{\pi }}\frac{g_{\pi }}{E_K+{\Delta }_P^{\ast }},$$

$${\Delta }_P^{\ast }=\frac{M_{B^{\ast }}^2-M_{B_s}^2-M_K^2}{2M_{B_s}}.$$

$$M_{B^{\ast }}(0^{+})-M_B\approx 400\mathrm{MeV}.$$

Full text

††preprint: FERMILAB-PUB-19-005-T

1. I Introduction

2. II Matrix elements and form factors

3. III Lattice-QCD calculation

4. III.1 Definitions

5. III.2 Actions and parameters

6. III.3 Interpolating operators, currents, and correlation functions

7. IV Analysis

8. IV.1 Analysis of the two-point correlation functions

9. IV.2 Extracting form factors from two- and three-point correlation functions

10. IV.3 Heavy bottom quark mass correction

11. IV.4 Chiral-continuum extrapolation

12. V Systematic error estimations

13. V.1 Chiral-continuum extrapolation errors

14. V.2 Current renormalization uncertainties

15. V.3 Lattice-scale uncertainties

16. V.4 Quark mass uncertainties

17. V.5 Uncertainties arising from the bottom quark mass correction

18. V.6 Finite volume effects

19. V.7 Summary of the statistical and systematic error budgets

20. VI Continuum form factors

21. VI.1 $z$ parametrization of form factors

22. VI.2 $z$-parametrization fit and extrapolation

23. VI.3 Comparison with existing results

24. VII Phenomenological applications

25. VII.1 Decay rate

26. VII.2 Forward-backward asymmetry

27. VII.3 Lepton polarization asymmetry

28. VII.4 Ratio of the $B_{s}\to K\ell\nu$ and $B_{s}\to D_{s}\ell\nu$ form factors

29. VIII Summary and Outlook

30. A $B_{s}\to K\ell\nu$ form factors in SU(2) chiral perturbation theory

31. B Reconstructing the $B_{s}\to K\ell\nu$ form factors

32. B.1 Reconstructing the form factors as functions of $z$

33. B.2 Reconstructing the form factors as functions of $q^{2}$

34. B.3 Dealing with the near zero eigenvalue in the covariance matrix

35. C $B_{s}\to K\ell\nu$ differential decay rate bin tables

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}\to K\ell\nu$ depends on the same matrix element $V_{ub}$ as the decay $B\to\pi\ell\nu$. 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}\to K\ell\nu$ decay can enable a precise $|V_{ub}|$ determination. This, in turn, can provide a useful test of $|V_{ub}|$ determinations from the exclusive $B\to\pi\ell\nu$ and $\Lambda_{b}\to p\ell\nu$ processes, and, if consistent, a reduced error on $|V_{ub}|$ (exclusive) after combination.

On the experimental side, however, while BaBar del Amo Sanchez et al. (2011); Lees et al. (2012) and Belle Ha et al. (2011); Sibidanov et al. (2013) have published precise measurements of the differential decay rate for $B\to\pi\ell\nu$, no such measurements exist yet for $B_{s}\to K\ell\nu$. The branching fraction of the former decay is $(7.80\pm 0.27)\times 10^{-5}$ Tanabashi et al. (2018), as it is Cabibbo suppressed compared to final states with charm. As BaBar and Belle observed many more $B\bar{B}$ than $B_{s}\bar{B}_{s}$ events, it is not surprising that experimental measurements of the latter decay have not yet been reported. In contrast, the LHCb experiment at the CERN LHC collider observes decays of all $b$-flavored hadrons, including $B_{s}$ mesons. They are expected to publish the results of their ongoing $B_{s}\to K\ell\nu$ decay study within the coming year Ciezarek et al. (2017). The Belle II experiment Urquijo (2015), where the $e^{+}e^{-}$ collisions provide a cleaner environment than at the LHC, also expects to study this decay. The current plans are that Belle II will collect about $50~{}\text{ab}^{-1}$ at the $\Upsilon(4S)$ resonance (which decays predominantly into $B$-meson pairs), and $5~{}\text{ab}^{-1}$ at the $\Upsilon(5S)$, a rich source of $B_{s}$-meson pairs Urquijo (2015). Thus, we do not expect the experimental accuracy for Belle II’s future measurement of $B_{s}\to K\ell\nu$ decay rates to rival that of their expected results for $B\to\pi\ell\nu$, but we do expect this decay to be studied by Belle II.

This work is part of a broad study of flavor physics by the Fermilab Lattice and MILC Collaborations to determine a number of CKM matrix elements from semileptonic $K$ Bazavov et al. (2013a), $D_{(s)}$ Aubin et al. (2005), and $B_{(s)}$ Bernard et al. (2009); Bailey et al. (2009, 2012a, 2012b, 2014, 2015a, 2015b, 2015c, 2016); Du et al. (2016) decays using the asqtad $2+1$ flavor ensembles generated by the MILC Collaboration Bernard et al. (2001); Aubin et al. (2004); Bazavov et al. (2010a). These studies are currently being extended Bazavov et al. (2014); Gámiz et al. (2016); Primer et al. (2017); Gelzer et al. (2018a); Bazavov et al. (2018a) to use HISQ $2+1+1$ flavor ensembles Bazavov et al. (2010b, 2013b). These newer ensembles include ones with physical-mass Goldstone pions at several lattice spacings that significantly improve our control of the chiral limit. In order to provide a systematic mode by mode comparison of results obtained with the two sets of configurations, it is important to complete this analysis of $B_{s}\to K\ell\nu$.

The techniques used here are very similar to those employed in Ref. Bailey et al. (2015b), where the functional $z$ expansion was introduced. However, in this work, we use a subset of six MILC ensembles covering a range of lattice spacing $a$ between approximately 0.12 and 0.06 fm. Prior work used 12 ensembles including one with $a\approx 0.045$ fm.

The decay $B_{s}\to K\ell\nu$ has been studied by three other lattice-QCD groups, the HPQCD Collaboration Bouchard et al. (2014), the RBC and UKQCD Collaborations Flynn et al. (2015), and the ALPHA Collaboration Bahr et al. (2016), each choosing different actions for the $b$-quark and for the light sea and valence quarks. Other previous calculations of the $B_{s}\to K\ell\nu$ decay form factors are based on the relativistic quark model Faustov and Galkin (2013), light-cone sum rules Duplancic and Melic (2008); Khodjamirian and Rusov (2017), and next-to-leading-order (NLO) perturbative QCD Wang and Xiao (2012). In Sec. VI.3, we compare our results with the prior results. Preliminary reports on this study can be found in Refs. Liu et al. (2014) and Liu et al. (2018), 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.2. 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}\to K\ell\nu$ decay can be described via the Feynman diagram shown in Fig. 1.

The relevant hadronic matrix element can be written as

[TABLE]

where $\mathcal{V}^{\mu}\equiv\bar{u}\gamma^{\mu}b$ is the vector current, $p^{\mu}_{B_{s}}$ and $p^{\mu}_{K}$ are the $B_{s}$ and $K$ four-momenta, respectively, $M_{B_{s}}$ and $M_{K}$ are the corresponding meson masses, $q^{\mu}=p_{B_{s}}^{\mu}-p_{K}^{\mu}$ is the momentum transferred to the lepton pair, and $f_{+}(q^{2})$ and $f_{0}(q^{2})$ are the vector and scalar form factors corresponding to the exchange of $1^{-}$ and $0^{+}$ particles. These two form factors are subject to a kinematic constraint:

[TABLE]

which eliminates the spurious pole at $q^{2}=0$ in Eq. (1). The tensor form factor $f_{T}$ parametrizes the hadronic matrix element of the tensor current $T^{\mu\nu}=i\bar{u}\sigma^{\mu\nu}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}\to K\ell\nu$ can be written as

[TABLE]

in the $B_{s}$ meson rest frame. Here $G_{F}$ is the Fermi constant, $V_{ub}$ is an element of the CKM matrix, $m_{\ell}$ is the lepton mass, and $\theta_{\ell}$ 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 $\tau$, 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 definitions of form factors and correlation functions are given in Sec. III.1. The lattice actions and simulation parameters are described in Sec. III.2. The lattice interpolating operators, currents, and correlation functions are presented in Sec. III.3.

III.1 Definitions

For lattice calculations and Heavy Quark Effective Theory (HQET), it is convenient to work in the $B_{s}$ rest frame and introduce the $B_{s}$ four-velocity

[TABLE]

The square of the lepton momentum transfer $q^{2}$ can then be expressed as

[TABLE]

where $E_{K}=p_{K}\cdot v$ is the kaon energy. Defining

[TABLE]

as the projection of the kaon momentum in the direction perpendicular to $v^{\mu}$ and using Eq. (5), one can rewrite the matrix element Eq. (1) in terms of the form factors $f_{\parallel}(E_{K})$ and $f_{\perp}(E_{K})$ as

[TABLE]

The relations to the original form factors $f_{+}$ and $f_{0}$ are given by

[TABLE]

The kinematic constraint, Eq. (2), is automatically satisfied in Eq. (8).

In the $B_{s}$ rest frame, which we use throughout the lattice-QCD calculation, the form factors $f_{\parallel}$ and $f_{\perp}$ are related to the temporal and spatial components of the matrix element of the vector current $\mathcal{V}^{\mu}$ via

[TABLE]

Note that there is no summation over the superscript $i$ in Eq. (9b). The continuum-QCD current is related to the lattice current operator $V^{\mu}$ by a multiplicative renormalization factor, i.e.,

[TABLE]

The lattice current $V^{\mu}$ is defined in Sec. III.3, below. We use a mostly nonperturbative method to compute $Z_{V_{\mu}}$. The details are explained in Sec. III.3.

The desired matrix elements (and hence form factors) can be calculated from suitably defined two- and three-point correlation functions:

[TABLE]

where $\mathcal{O}_{B_{s}}$ and $\mathcal{O}_{K}$ are lattice interpolating operators, which are defined in Sec. III.3, below. Further, $\bm{p}_{K}$ is the kaon spatial momentum, whose components in a finite volume are integer multiples of $2\pi/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. (1) in lattice QCD is the following:

For each ensemble:

1. (i)

Determine the lattice $B_{s}$ meson masses, kaon masses and energies from the lattice two-point correlation functions. 2. (ii)

Determine the lattice form factors $f_{\parallel}^{\mathrm{lat}}$ and $f_{\perp}^{\mathrm{lat}}$ at several discrete kaon momenta $\bm{p}_{K}$ from the two- and three-point correlation functions. 3. (iii)

Obtain the renormalized form factors by matching the lattice current to the continuum as in Eq. (10). 2. 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_{\parallel}$ and $f_{\perp}$ as functions of the kaon recoil energy $E_{K}$ in the interval covered by the simulation, roughly $0.5$ GeV $\lesssim E_{K}\lesssim 1$ GeV. 3. 3.

Construct the continuum form factors $f_{+}$ and $f_{0}$ from $f_{\parallel}$ and $f_{\perp}$ via Eq. (8) 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$.

III.2 Actions and parameters

We use lattice gauge configurations with $N_{f}=2+1$ flavors generated by the MILC Collaboration Bazavov et al. (2010a); Bernard et al. (2001); Aubin et al. (2004). These configurations include two degenerate dynamical light quarks, acting as $u$ and $d$ quarks, and one heavier, $s$, quark. The gluon fields are simulated with the one-loop improved Lüscher-Weisz action Lüscher and Weisz (1985). The $a^{2}$ tadpole-improved staggered action (asqtad) Blum et al. (1997); Lepage (1998); Lagaë and Sinclair (1999); Lepage (1999); Orginos and Toussaint (1999); Orginos et al. (1999); Bernard et al. (2000a) is used for generating dynamical light quarks ($u$, $d$, and $s$). Reference Bazavov et al. (2010a) is a review of simulations and formalism of improved staggered quarks.

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 Sheikholeslami and Wohlert (1985) with the Fermilab interpretation El-Khadra et al. (1997).

Some of the parameters used to generate the configurations are listed in Table 1. Six ensembles with three different lattice spacings, $a\approx 0.12$, 0.09, and 0.06 fm, are used. For each lattice spacing, we have dynamical sea quarks with light-to-strange quark mass ratio $m^{\prime}_{l}/m^{\prime}_{h}=0.1$111In this paper, we use primed quantities to denote the sea quarks and the unprimed for the valence quark.. For the intermediate lattice spacing $a\approx 0.09~{}\textrm{fm}$, we have three additional values of $m^{\prime}_{l}/m^{\prime}_{h}=0.05,0.15$, and 0.2 to provide results for the chiral extrapolation. The subset of ensembles used for the analysis is based on experience from previous semileptonic form factor analyses Bailey et al. (2015b) and Bailey et al. (2016). 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 2. The valence light quarks are degenerate with the sea quarks, i.e., $am_{l}=am_{l}^{\prime}$; 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}^{\prime}$. The heavy $b$ quark Wilson fermions with SW lattice action are controlled by the hopping parameter $\kappa$ and the clover coefficient of the SW action $c_{\mathrm{sw}}$. We use $\kappa_{b}^{\prime}$ to denote the values used in the computation. We use the tadpole-improved tree-level value for $c_{\mathrm{sw}}=u_{0}^{-3}$, with $u_{0}$ listed in Table 1. The parameter $d_{1}$ is used for the correlation function generation and will be explained later in Sec. III.3.

Table 3 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_{1}^{2}F(r_{1})=1.0$ Sommer (1994); Bernard et al. (2000b). 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_{\pi}$: $r_{1}=0.3117(22)~{}\mathrm{fm}$ Bazavov et al. (2010a, 2012). The physical value $\kappa_{b}$ Bailey et al. (2014), corresponding to the physical $b$-quark mass, and the critical value $\kappa_{\mathrm{crit}}$, corresponding to the zero quark masses in the SW action on each ensemble, are also listed in Table 3. They will be used only for correcting the $b$-quark masses as will be discussed in Sec. IV.3. The Goldstone pion mass $M_{\pi}$ and the root-mean-square (RMS) pion mass $M_{\pi}^{\textrm{RMS}}$ are listed in the last two columns of Table 3.

III.3 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. (11). For the kaon, the local pseudoscalar interpolating operator is used

[TABLE]

where $\chi(t,\bm{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 Bailey et al. (2009); Wingate et al. (2003); Kawamoto and Smit (1981)

[TABLE]

where $\psi(t,\bm{y})$ is the four-component $b$-quark field, and $S(\bm{x},\bm{y})$ is a spatial smearing function. We use two smearing functions for the $B_{s}$ meson. One is the local $S(\bm{x},\bm{y})=\delta(\bm{x}-\bm{y})$. The other one is the ground-state 1S wave function of the Richardson potential Bazavov et al. (2012).

The lattice vector current operator in Eqs. (10) and (11c) is defined as in Refs. Bailey et al. (2009) and Wingate et al. (2003)

[TABLE]

where the rotated $b$-quark field $\Psi$, defined by

[TABLE]

removes $O(a)$ discretization effects from the current El-Khadra et al. (1997). Here $\bm{D}_{\mathrm{lat}}$ is a symmetric nearest-neighbor covariant difference operator. The coefficient $d_{1}$, shown in Table 2, 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_{\mu}}$, needed to match the lattice vector current to its continuum counterpart (see Eq. (10)), is determined using a mostly nonperturbative renormalization procedure Harada et al. (2002); El-Khadra et al. (2001):

[TABLE]

where $Z_{V^{4}_{bb}}$ and $Z_{V^{4}_{ll}}$ are the renormalization factors for the flavor-diagonal $b$- and light-quark temporal vector currents that are calculated nonperturbatively in Ref. Bailey et al. (2015b) and listed in Table 4. The remaining flavor-off-diagonal parameters $rho_{V^{\mu}}$ are calculated to one-loop order in perturbation theory, separately from this analysis, and also listed in Table 4. In order to reduce subjectivity in our analysis, we employed a blinding procedure in the form of a small multiplicative offset applied to the $\rho$ factors and known to only two of the authors. This blinding factor was subsequently disclosed and removed only after the analysis choices were finalized.

In the generation of the correlation functions defined in Eqs. (11), (12), (13), and (14), we increase statistics by repeating the calculation at $N_{\mathrm{src}}$ source times evenly distributed in the $N_{t}$ direction. The three-point correlation functions are generated with two adjacent temporal source-sink separations: $T=T_{\mathrm{sink}}$ and $T=T_{\mathrm{sink}}+1$. Both $N_{\mathrm{src}}$ and $T_{\mathrm{sink}}$ are listed in Table 5. For the kaon recoil momenta we include the following lowest possible values: $\bm{p}_{K}/(2\pi/N_{s})=(0,0,0),(1,0,0),(1,1,0),(1,1,1)$ and $(2,0,0)$. In practice, the largest momentum $\bm{p}_{K}=2\pi(2,0,0)/N_{s}$ is too noisy and is excluded from the analysis.

IV Analysis

With lattice correlation functions in hand, we follow the steps outlined near the end of Sec. III.1 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 Wingate et al. (2003):

[TABLE]

with

[TABLE]

The $(-1)^{n(t+1)}$ and $(-1)^{n(T-t)}$ terms in Eq. (17) arise because with our choice for the light-quark valence action the interpolating operators also generate opposite-parity (scalar) states. The overlap factors $Z_{B_{s}}^{(n)}$ and $Z_{K}^{(n)}(\bm{p}_{K})$ describe the overlap of the interpolating operators with the states $|{B_{s}}^{(n)}\rangle$ and $|K^{(n)}(\bm{p}_{K})\rangle$, respectively, while the $D_{mn}^{\mu}$ contain the desired matrix element.

In Sec. IV.1, we extract the meson masses, and the overlap factors $Z_{B_{s}}^{(n)}$ and $Z_{K}^{(n)}(\bm{p}_{K})$ from the two-point correlation functions. We explain how the lattice form factors are extracted from the two- and three-point correlation functions in Sec. IV.2. We briefly describe the heavy $b$-quark mass corrections in Sec. IV.3. The chiral-continuum fit function and extrapolation are described in Sec. IV.4.

IV.1 Analysis of the two-point correlation functions

The $B_{s}$-meson masses, kaon masses, kaon energies, and $B_{s}$ and kaon overlap factors are obtained from fitting the two-point correlation functions to the functional forms in Eqs. (17a) and (17b).

As listed in Table 5, 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_{\mathrm{max}}$ is the largest value of $t$ where the fractional error in the correlation function is smaller than $3\%$; $t_{\mathrm{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. Bazavov et al. (2016). The fit ranges $[t_{\mathrm{min}},t_{\mathrm{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. (17) therefore refers to a fit function with the pseudoscalar ground state plus $N-1$ of its radial excitations and $N$ scalar states. Our central value fits have $N=3$. The prior central values for the ground state energies and overlap factors are guided by the effective mass and effective amplitude evaluated at large times $t$. The effective mass $m_{\mathrm{eff}}$ and effective amplitude $Z_{\mathrm{eff}}$ are constructed from the two-point correlation functions via

[TABLE]

Here $C_{2}$ stands for the lattice two-point correlation function for the kaon or $B_{s}$ meson. The prior central values for $M_{B_{s}}^{(0)}$, $M_{K}^{(0)}\equiv E_{K}^{(0)}$, $Z_{B_{s}}^{(0)}$, $Z_{K}^{(0)}(0)$, and $Z_{K}^{(0)}(\bm{p}_{K})$ are set according to Eq. (19) and the widths are set to be 0.1 or larger in lattice units. The prior central values for $M_{B_{s}}^{(n\neq 0)}$ and $E_{K}^{(n\neq 0)}$ are set using the energy difference between ground states and the corresponding excited states from the PDG Patrignani et al. (2016) values as a guide wherever available and the widths are set to be 0.1 or larger in lattice units. The prior central values for $Z_{B_{s}}^{(1)}$, $Z_{K}^{(1)}(0)$, and $Z_{K}^{(1)}(\bm{p}_{K})$ are set using the $N=1$ fit results as a guide and the widths are set to be 0.1 or larger in lattice units. The prior central values for $Z_{B_{s}}^{(2)}$, $Z_{K}^{(2)}(0)$, and $Z_{K}^{(2)}(\bm{p}_{K})$ are set using $Z_{B_{s}}^{(0)}$, $Z_{K}^{(0)}(0)$, and $Z_{K}^{(0)}(\bm{p}_{K})$ as a guide and the widths are set to be 0.1 or larger in lattice units. Finally the prior central values for $Z_{B_{s}}^{(3,4,5)}$, $Z_{K}^{(3,4,5)}(0)$, and $Z_{K}^{(3,4,5)}(\bm{p}_{K})$ are set to be 0.1 and the widths are set to be 1.0 or larger in lattice units. The prior widths in general are set to be large enough so that no bias is introduced in the fits. An example of the $B_{s}$ effective mass, prior, and fit result is shown in the left panel of Fig. 2. The corresponding kaon effective mass has smaller oscillations and much smaller errors as shown in the right panel of Fig. 2. Our fit results are stable over a range of $t_{\mathrm{min}}$ choices and consistent with results from $N=2$ fits. We find that the lattice correlation functions are precise enough to determine the first excited and opposite-parity, $N=2$, states. Including extra $N=3$ excited states better stabilizes the errors of fit posteriors.

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 6. 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

[TABLE]

to study momentum-dependent discretization errors. As illustrated in Fig. 4 we find that the $E_{K}$ and $Z_{K}^{(0)}(\bm{p}_{K})$ satisfy Eq. (20) albeit with increasing statistical errors at higher momenta. We therefore use the continuum relations for the kaon energies and $Z_{K}$ factors whenever possible.

IV.2 Extracting form factors from two- and three-point correlation

functions

The form factors are related to the semileptonic matrix elements via Eq. (9), and the lattice matrix elements are contained in the three-point correlation function as in Eqs. (17c) and (18). To get the lattice form factors $f_{\parallel,\perp}^{\textrm{lat}}$, we fit the two- and three-point correlation functions together. In particular, we perform combined two- and three-point correlation-function fits according to Eqs. (17) and (18) with $N=3$. The three-point fit ranges are chosen to be $[t_{\text{min}}^{K},T-t_{\text{min}}^{B_{s}}]$ with $T=T_{\mathrm{sink}}$ or $T_{\mathrm{sink}}+1$. The parameters to be fitted are $M_{B_{s}}^{(n)}$, $E_{K}^{(n)}$, $Z_{B_{s}}^{(n)}$, $Z_{K}^{(n)}(\bm{p}_{K})$, and $D_{mn}^{\mu}$. The prior central values for the $B_{s}$-meson and kaon masses, $Z_{B_{s}}^{(0)}$, and $Z_{K}^{(0)}(\bm{p}_{K})$ are chosen as the posteriors of the two-point correlator fits. The kaon energies and $Z_{K}^{(0)}(\bm{p}_{K})$ are constrained according to Eq. (20). The $Z_{B_{s}}^{(n\neq 0)}$ and $Z_{K}^{(n\neq 0)}(\bm{p}_{K})$ central values are taken to be the same size as $Z_{B_{s}}^{(0)}$ and $Z_{K}^{(n\neq 0)}(0)$. The priors for $D_{00}^{\mu}$ are guided by the constructed ratio $\bar{R}_{3,0}^{B_{s}\to K}(t,T)$ defined in Refs. Bailey et al. (2015b) and Bailey et al. (2016). 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}^{\mu}$ are chosen to be $0.1\pm 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.1.

Figure 5 shows that the fitted ${f_{\parallel}}^{\textrm{lat}}$ coming from the combined fit is in slight tension with the constructed ratio $\bar{R}_{3,0}^{B_{s}\to K}(t,T)$ defined in Refs. Bailey et al. (2015b) and Bailey et al. (2016). This small difference comes from excited state contributions still present in the ratio but accounted for in the fit method 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_{\parallel}$ and $f_{\perp}$ are obtained from $D_{00}^{\mu}$ and $E_{K}$ according to Eqs. (9) and (18c), after adding the renormalization factors as in Eq. (9).

IV.3 Heavy bottom quark mass correction

The heavy valence $b$ quark is simulated with the Sheikholeslami-Wohlert (SW) action Sheikholeslami and Wohlert (1985) with the Fermilab interpretation El-Khadra et al. (1997). The $b$-quark mass is controlled by the hopping parameter $\kappa_{b}$. The hopping parameter $\kappa_{b}^{\prime}$ used in the simulations differs slightly from the physical value $\kappa_{b}$ as can be seen in Tables 2 and 3. We need to correct the form factors to account for these small shifts. A detailed description of the $\kappa$ tuning analysis and results is provided in Appendix C of Ref. Bailey et al. (2014). We use the method described in Refs. Bailey et al. (2015b) and Bailey et al. (2016) to adjust the form factors to account for the slightly mistuned values of $\kappa_{b}$. The relative change in the form factors under small variations of the $b$-quark mass can be described as

[TABLE]

where $m_{2}$ is the physical $b$-quark kinetic mass, $m_{2}^{\prime}$ is the $b$-quark mass used in the production run. The slopes $\frac{\partial\ln f}{\partial\ln m_{2}}$ were determined in Ref. Bailey et al. (2015b). The corrections to the form factors are about 0.1–1.8% on different ensembles.

IV.4 Chiral-continuum extrapolation

The lattice form factors extracted from the correlation functions as described in Sec. IV.2 are obtained at the three finite lattice spacings and unphysical light-quark masses listed in Table 1. 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$\chi$PT) Aubin and Bernard (2006, 2007). Based on previous experience with the analyses of similar processes in Refs. Bailey et al. (2015b) and Bailey et al. (2016), this best describes the data. Heavy-quark discretization effects are also taken into account in the chiral-continuum extrapolation.

We employ the HMrS$\chi$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$\chi$PT formula; the large kaon energy is integrated out, and its effects are absorbed into the low-energy constants (LECs). In Ref. Bailey et al. (2016), the conversion rules for $B\rightarrow K$ and $B\rightarrow\pi$ processes from SU(3) HMrS$\chi$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}\to K\ell\nu$. The details are presented in Appendix A.

The NLO expression for $B_{s}\to K\ell\nu$ form factors in the SU(2) hard-kaon limit that we obtain is

[TABLE]

where $P=~{}\parallel$ or $\perp$, $\delta f_{P,\mathrm{logs}}^{\mathrm{SU(2)}}$ are the non-analytic contributions from the light-quark mass and lattice spacing, and the $\chi$ variables are dimensionless. They are defined in Eqs. (51) and (55). The leading-order factor is

[TABLE]

where $f_{\pi}$ is the decay constant involved, and $g_{\pi}$ is the $B^{*}B\pi$ coupling constant222$SU(3)$ breaking effects renormalize the $g_{\pi}/f_{\pi}$ ratio; however, since it results in a overall multiplicative factor, it has been reabsorbed in the fitting coefficients.. The $\Delta_{P}^{*}$ term takes the pole contribution into account and is determined by requiring $f_{\parallel}$ and $f_{\perp}$ to have the same poles as the physical form factors $f_{0}$ and $f_{+}$, respectively. This is reasonable because, by Eq. (8), $f_{\parallel}$ is dominated by the $0^{+}$ contributions of $f_{0}$, and $f_{\perp}$, by the $1^{-}$ contributions of $f_{+}$ in the $q^{2}$ range considered. Using Eq. (5), one obtains the exact expression for $\Delta_{P}^{*}$:

[TABLE]

The vector meson (with $J^{P}=1^{-}$) has been experimentally measured Tanabashi et al. (2018) to be $M_{B^{*}}=5324.65(25)~{}\mathrm{MeV}$; the scalar $B^{*}$ meson (with $J^{P}=0^{+}$) has not been observed experimentally, but a lattice calculation Gregory et al. (2011) estimates the mass difference between $0^{+}$ and $0^{-}$ states to be around $400~{}\mathrm{MeV}$:

[TABLE]

The vector-meson mass $M_{B^{*}}$ is below the $B\pi$ production threshold that is involved in the $B_{s}\to K\ell\nu$ decay, and the scalar-meson mass $M_{B^{*}}(0^{+})$ is above the threshold. The inclusion of the scalar pole and its exact location have little impact on the chiral fit results but stabilizes the form factor extrapolations.

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 for the NNLO fit functions are

[TABLE]

where the heavy-quark discretization effects are modeled with $\delta f_{\mathrm{HQ}}$. The mismatch functions $f_{E,X,Y,B,3}$ are defined in the appendix of Ref. Bazavov et al. (2012). The next-to-leading-order (NLO) analytic term $f_{P,\mathrm{NLO}}$ was defined previously in Eq. (22). A Bayesian method is used in the chiral-continuum fit. The priors are listed in Table 7. The fit results using the NNLO fit function in Eq. (26) are used as the central fit and are shown in Fig. 7.

V Systematic error estimations

The chiral-continuum extrapolated form factors are given in Sec. IV. The statistical-fit errors, which are propagated through each step of the analysis already include the effects of NNLO terms in the chiral expansion as well as light- and heavy-quark discretization. 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~{}{\rm GeV}^{2}\lesssim q^{2}\lesssim 24~{}{\rm GeV}^{2}$.

V.1 Chiral-continuum extrapolation errors

Our central fit uses the NNLO SU(2) hard kaon HMrS$\chi$PT fit function described in Eq. (26). 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. (22), 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}\gtrsim 19~{}\mathrm{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$\chi$PT formula affect the fits, we performed the fit with soft-kaon HMrS$\chi$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$\chi$PT formula are extremely complicated. We therefore did not perform any trial fits with the SU(3) formula. Nevertheless, from previous experience, Bailey et al. (2015b, 2016), SU(3) HMrS$\chi$PT typically does not provide a good description of the data.

Our results with kaon momentum up to $2\pi(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 $\bm{p}_{K}=2\pi(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}\gtrsim 19~{}\mathrm{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.

V.2 Current renormalization uncertainties

The mostly nonperturbative renormalization procedure, described in Eq. (10), used to renormalize the matrix elements (and hence the form factors) requires, as inputs, the factors $Z_{V^{4}_{bb}}$, $Z_{V^{4}_{ll}}$, $\rho_{V^{4}}$ and $\rho_{V^{i}}$. We estimate the error on $f_{\parallel},f_{\perp}$ due to the uncertainties of the nonperturbatively determined $Z_{V^{4}_{bb}}$ and $Z_{V^{4}_{ll}}$ by varying their central values by one standard deviation in each direction. As expected, the resulting changes in the form factors are small, yielding errors on $f_{\parallel},f_{\perp}$ in the range of 0.2–0.3%.

For $\rho_{V^{4}}$ and $\rho_{V^{i}}$, the dominant source of error is the truncation at one-loop order in perturbation theory. As seen in Table 4, the one-loop corrections provided by $\rho_{V^{4}}$ and $\rho_{V^{i}}$ are small, with $\rho_{V^{4}}$ ($\rho_{V^{i}}$) deviating from unity by less than $1\%$ ($2.4\%$). Here, we adopt the estimate of the perturbative truncation error presented in Ref. Bailey et al. (2015b), which yields an uncertainty of $1\%$ on both $\rho_{V^{4}}$ and $\rho_{V^{i}}$. This estimate is consistent with the observed differences between nonperturbative Chakraborty et al. (2014) and perturbative El-Khadra et al. (2007) calculations of $\rho_{V^{4}}=\rho_{A^{4}}$, discussed in Ref. Chakraborty et al. (2014). In particular, the observed differences decrease in the continuum limit, as expected. Note, however, that the nonperturbative result Chakraborty et al. (2014) employs the HISQ action for the light quarks, while our one-loop results El-Khadra et al. (2007) employ the asqtad action. For this reason, the comparison is suggestive but not definitive.

V.3 Lattice-scale uncertainties

The dimensionful form factors $f_{\perp}$ and $f_{\parallel}$, and meson energies and masses are converted to physical units via the relative scales $r_{1}/a$ listed in Table 3 and the absolute scale $r_{1}=0.3117(22)$ fm Bazavov et al. (2012). 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.

V.4 Quark mass uncertainties

The continuum physical form factors are obtained by evaluating the chiral-continuum extrapolated functions, as discussed in Sec. IV.4 at the physical averaged $u$- and $d$-quark masses, namely $r_{1}m_{ud}^{\text{phys}}=0.000965(33)$, and the physical $s$-quark mass $r_{1}m_{s}^{\text{phys}}=0.0265(8)$ as determined by analyzing the light pseudoscalar meson spectrum Bazavov et al. (2010a). 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.

V.5 Uncertainties arising from the bottom quark mass correction

As explained in Sec. IV.3, 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 $\kappa_{b}$-correction factors and the tuned $\kappa_{b}$ values are taken into account by following the procedure described in Ref. Bailey et al. (2015b). A $q^{2}$-independent $0.4\%$ error due to tuning $\kappa_{b}$ is assigned to both $f_{+}$ and $f_{0}$.

V.6 Finite volume effects

Finite-volume effects, estimated by comparing infinite-volume integrals with finite sums in HMrS$\chi$PT, are negligibly small Bailey et al. (2015b, 2016), so they are omitted from the total error budget.

V.7 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$\chi$PT used to perform the extrapolation is valid only for moderate $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. Bailey et al. (2015b) we therefore treat them as constants in $q^{2}$ when propagating them to the $z$-parametrization fit in Sec. VI.2. 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. Bailey et al. (2015b) 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}\gtrsim 17~{}\text{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) Bourrely et al. (2009) and the fitting procedure and extrapolation technique was first introduced in our previous $B\to\pi\ell\nu$ paper Bailey et al. (2015b).

In Sec. VI.1, we briefly review the $z$-parametrization and give the expansion form used in the analysis. In Sec. VI.2, we present the extrapolated continuum form factors in the whole kinematically allowed region. The results are shown in Table 10, and Figs. 12 and 13. A comparison with results of other groups is presented in Sec. VI.3.

VI.1 $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 Eden et al. (1966) imply that the $B_{s}\to K\ell\nu$ semileptonic form factors are real analytic functions333 An analytic function $f(x)$ is real analytic if it satisfies $f(x^{*})=(f(x))^{*}$. If $f(x)$ is a real analytic function with a branch point at $x_{0}$, then $f(x)$ is real for $x<x_{0}$ and its discontinuity across the cut is purely imaginary: $f(x+i\epsilon)-f(x-i\epsilon)=2i\mathrm{Im}f(x+i\epsilon)$.

in the complex $q^{2}$-plane with a cut from $q^{2}>t_{\mathrm{cut}}$ to $\infty$, except at physical poles below $t_{\mathrm{cut}}$. The parameter $t_{\mathrm{cut}}$ is the particle-pair-production threshold. For $B_{s}\to K\ell\nu$, this is

[TABLE]

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 Lepage and Brodsky (1980); Akhoury et al. (1994), 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_{\mathrm{cut}}$, the form factors have the following scaling behavior

[TABLE]

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

[TABLE]

where $t_{0}$ is a parameter that can be chosen to optimize the mapping. The maximum momentum transfer allowed in the semileptonic $B_{s}\to K\ell\nu$ decay is defined as

[TABLE]

for convenience. This conformal mapping was first considered in Ref. Meiman (1963) and further developed and used to get model-independent constraints, usually called “unitarity bounds”, on form factors in Ref. Okubo (1971). A stronger constraint based on heavy-quark power counting was derived in Ref. Becher and Hill (2006). The conformal transformation Eq. (29) maps the physical semileptonic region $0\leq q^{2}\leq 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}=\pm\infty$ to $z=1$, and $q^{2}=t_{\mathrm{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

[TABLE]

The solution for $t_{0}$ is

[TABLE]

This mapping is schematically shown in Fig. 11 with small lepton masses ignored and with the optimized $t_{0}$ as defined in Eq. (32).

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|\leq 0.205$ for $B_{s}\to K\ell\nu$, this parametrization converges quickly. Table 8 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) Boyd et al. (1995) and by Bourrely, Caprini and Lellouch (BCL) Bourrely et al. (2009). Here we use the BCL parametrization as given by

[TABLE]

The factors $1/({1-q^{2}/m_{B^{*}}^{2}})$ take the poles into account and ensure the asymptotic scaling, $f(q^{2})\sim 1/q^{2}$ at large $q^{2}$. Moreover, the scaling condition of Eq. (28) near $t_{\mathrm{cut}}$ is also enforced for $f_{+}$. Note that Eq. (28) in the $q^{2}$-plane imply the following relation

[TABLE]

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.

VI.2 $z$-parametrization fit and extrapolation

We use Eq. (33) to perform the $z$-parametrization fit to our chiral-continuum-extrapolated form factor results obtained in Secs. IV.4 and V. The vector pole $M_{B^{*}(1^{-})}$ is taken to be $M_{B^{*}(1^{-})}=5.32465(25)\mathrm{GeV}$ Tanabashi et al. (2018), and the above threshold scalar pole $M_{B^{*}(0^{+})}$ is taken to be the theoretically predicted value $M_{B^{*}(0^{+})}=5.68\mathrm{GeV}$ Gregory et al. (2011). The parameter $t_{0}$ is chosen as in Eq. (32), and the corresponding value for the $B_{s}\to K\ell\nu$ process is $16.5~{}\mathrm{GeV}^{2}$. Table 9 lists the relevant meson masses used in the $z$-parametrization fit.

The functional method introduced in Ref. Bailey et al. (2015b) is used to perform the $z$-parametrization fit, where, following Ref. Bailey et al. (2015b), we take as inputs the results from the chiral-continuum extrapolation and systematic error analysis as presented in Sec. V.7).

Our preferred (central) fit has $K=4$, where $K$ is the number of terms in the expansion in Eq. (33). The results of this fit are shown in Table 10. These can be used to reconstruct the final form factors as described in Appendix B. We arrive at this preferred fit choice by first simultaneously fitting the form factors $f_{+}$ and $f_{0}$ with $K=2$ and without constraining the $z$-parametrization parameters $b_{i}^{+,0}$ in Eq. (33). The coefficients $b_{0}^{+}$ and $b_{0}^{0}$ are well determined, but the quality of this fit is poor. When increasing $K$ from 2 to 3, the quality of the fit improves, and all the $b_{i}^{+,0}$ coefficients can be determined well. The kinematic constraint Eq. (2) is satisfied within errors444Note that the kinematic constraint is automatically satisfied in Eq. (8) before taking the extrapolation as is being done in this section. After the extrapolation, this constraint is not guaranteed if not imposed in the fit.. Enforcing this kinematic constraint, as explained below, further improve the $f_{+}$ form-factor fit. The fit parameters also satisfy the unitarity condition Bourrely et al. (2009) and the condition estimated from heavy-quark power counting Becher and Hill (2006). 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

[TABLE]

with width $\epsilon=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.

VI.3 Comparison with existing results

Several other groups have also calculated the same form factors. We note that Refs. Bouchard et al. (2014) and Flynn et al. (2015) use the $B_{s}K$ threshold instead of $B\pi$ in their implementation of the $z$ parametrization. Since the $z$-parameter, by definition (see Eq. (29)), depends on the threshold ($t_{\rm cut}$), we cannot directly compare the $z$-dependence of our form factors with those of Refs. Bouchard et al. (2014) and Flynn et al. (2015). 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 Bouchard et al. (2014) 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 Flynn et al. (2015) use (2+1)-flavor-domain-wall fermions for the sea quarks and light valence quarks, and a variant Lin and Christ (2007); Christ et al. (2007) of the Fermilab action for the heavy $b$-quark. The ALPHA Collaboration Bahr et al. (2016) uses leading-order lattice HQET to get the form factors at one point, $q^{2}=22.12~{}\mathrm{GeV}^{2}$. While our results are consistent with those from Refs. Flynn et al. (2015) and Bahr et al. (2016), they are in tension with HPQCD’s results Bouchard et al. (2014). We note that Ref. Bouchard et al. (2014) 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\to K\ell^{+}\ell^{-}$ decay in Ref. Bouchard et al. (2013), 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. Bailey et al. (2016). However, unlike the case at hand, the form factors of Ref. Bailey et al. (2016) are not in significant tension with HPQCD’s results of Ref. Bouchard et al. (2013). We see that the tension between our $B_{s}\to K\ell\nu$ form factor results and those of Ref. Bouchard et al. (2014) increases with decreasing $q^{2}$ to roughly 2.3$\sigma$ at $q^{2}=0$. The RBC and UKQCD Collaborations Flynn et al. (2015), 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 Duplancic and Melic (2008); Khodjamirian and Rusov (2017), a relativistic quark model Faustov and Galkin (2013), and NLO perturbative QCD Wang and Xiao (2012).

VII Phenomenological applications

The angular-dependent differential decay rate for $B_{s}\to K\ell\nu$ is given in Eq. (3). One can construct at most three independent observables from there. In the following, we will consider the differential decay rate $d\Gamma/dq^{2}$ in Sec. VII.1, the forward-backward asymmetry $A_{FB}^{\ell}(q^{2})$ in Sec. VII.2, and the lepton polarization asymmetry $A_{pol}^{\ell}(q^{2})$ in Sec. VII.3. The latter two quantities are sensitive to the mass of the final-state charged lepton. In Sec. VII.4, we also construct the ratios of the scalar and vector form factors between the $B_{s}\to K\ell\nu$ and $B_{s}\to D_{s}\ell\nu$ decays.

VII.1 Decay rate

The differential decay rate can be obtained from Eq. (3) by integrating over the angle $\theta_{\ell}$, which yields

[TABLE]

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}\to K\mu\nu$ and $B_{s}\to K\tau\nu$.

One can also explore the ratio of the differential decay rates

[TABLE]

Figure 16 shows the prediction for $R^{\tau/\mu}(q^{2})$.

The total decay rate is given by

[TABLE]

with $q_{\mathrm{max}}^{2}=t_{-}=(M_{B_{s}}-M_{K})^{2}$, as in Eq. (30). The numerical results for $\Gamma/|V_{ub}|^{2}$ are

[TABLE]

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

[TABLE]

which takes the correlations between the form factors into account and is more precise than directly using Eq. (39).

VII.2 Forward-backward asymmetry

The forward-backward asymmetry, $A_{FB}$, which depends on the linear $\cos\theta_{\ell}$ term in Eq. (3), is given by

[TABLE]

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

[TABLE]

The normalized forward-backward asymmetry is given by

[TABLE]

and the corresponding numerical values are

[TABLE]

VII.3 Lepton polarization asymmetry

The normalized lepton polarization asymmetry is defined as

[TABLE]

from the differential decay rates with definite lepton helicity Meiß ner and Wang (2014)

[TABLE]

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_{\text{pol}}^{\mu}$ is close to one in the whole $q^{2}$ range. Here we provide the normalized lepton polarization asymmetry $A_{\text{pol}}^{\mu}$ and $A_{\text{pol}}^{\tau}$ as functions of $q^{2}$ in Fig. 18.

VII.4 Ratio of the $B_{s}\to K\ell\nu$ and

$B_{s}\to D_{s}\ell\nu$ form factors

We also calculate the ratios of the scalar and vector form factors between the $B_{s}\to K\ell\nu$ and $B_{s}\to D_{s}\ell\nu$ 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}\to D_{s}\ell\nu$ form factors from our previous papers Bailey et al. (2012b, 2015a). Form factor ratios, $f_{+,0}^{2012}(B_{s}\to D_{s})/f_{+,0}^{2012}(B\to D)$, and the $B\to D\ell\nu$ form factors, $f_{+,0}^{2015}(B\to D)$, are calculated in Refs. Bailey et al. (2012b) and Bailey et al. (2015a), respectively. They are shown in Fig. 19.

The $B_{s}\to D_{s}\ell\nu$ form factor can be reconstructed via

[TABLE]

With the reconstructed $B_{s}\to D_{s}\ell\nu$ form factors $f_{+,0}^{\mathrm{reco}}(B_{s}\to D_{s})$ shown in Fig. 20, we obtain the form-factor ratios, $f_{+,0}(B_{s}\to K)/f_{+,0}^{\mathrm{reco}}(B_{s}\to D_{s})$, shown in Fig. 21 as functions of $q^{2}$ and in Fig. 22 as functions of $w$. Although the 2012 analysis was carried out on a subset of the ensembles used in the 2015 analysis, we neglect any correlations in the two form factors in Eq. (47). Here $q^{2}$ is the usual square of the lepton momentum transfer as defined in Eq. (5). The recoil parameter $w$ for $B_{s}\to D_{s}\ell\nu$ is defined as

[TABLE]

and the corresponding one for the $B_{s}\to K\ell\nu$ is defined by replacing $M_{D_{s}}$ with $M_{K}$. The relation between $w$ and $q^{2}$ in Eq. (48), and the kinematically allowed regions for the two types of processes are shown in Fig. 23. The ratios constructed with different parameters $q^{2}$ and $w$ as shown in Figs. 21 and 22 allow us to probe the different $B_{s}\to K\ell\nu$ form factor regions.

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}\to K\ell\nu$. We present predictions of the differential decay rate (divided by $|V_{ub}|^{2}$) for both light ($e$ or $\mu$) or heavy ($\tau$) final-state leptons. Once the experimental data become available, our form factors can be used to determine $|V_{ub}|$, which can then be compared to and, if consistent, combined with the $|V_{ub}|$ determinations from other exclusive decay processes. Hence they may help shed light on the discrepancy with $|V_{ub}|$ from inclusive decays $B\to X_{u}\ell\nu$ and, perhaps, contribute to evidence for new physics beyond the Standard Model by enabling more stringent tests of the CKM paradigm. Other quantities of phenomenological interest include the forward-backward asymmetry $A_{\text{FB}}^{\ell}(q^{2})$ and the lepton polarization asymmetry $A_{\text{pol}}^{\ell}(q^{2})$. We also present ratios of the form factors $f_{+}$ and $f_{0}$ for $B_{s}\to K\ell\nu$ and $B_{s}\to D_{s}\ell\nu$ 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}\to K\ell\nu$, this process is under investigation by the LHCb experiment, and will be studied by the Belle II Collaboration when they run at the $\Upsilon(5\text{S})$ resonance, which is a copious source of $B_{s}$ and $\bar{B}_{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 Bazavov et al. (2010b, 2013b) 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. Gelzer et al. (2018a) and Gelzer et al. (2018b). 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\approx 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 Bazavov et al. (2018b) employing the heavy HISQ method for the $B$- and $D$-meson decay constants has reached unprecedented precision. We have recently started to generate the correlation functions for this project. In summary, with the improvements outlined above, we expect, in the coming years, to obtain the form factors for $B_{s}\to K\ell\nu$ (and related decays) with percent level precision, at least in the low recoil region of the phase space.

Acknowledgements.

We thank Jon A. Bailey for participating in the early stage of the project. We thank Chris M. Bouchard for discussions and comments on the form factor comparison section. We thank members of the LHCb Collaboration for discussions and in particular Svende Braun, Marta Calvi, and Mika A. Vesterinen for sharing with us their analysis status and the preferred $q^{2}$ bins. Computations for this work were carried out with resources provided by the USQCD Collaboration, the National Energy Research Scientific Computing Center, the Argonne Leadership Computing Facility, the Blue Waters sustained-petascale computing project, the National Institute for Computational Science, the National Center for Atmospheric Research, the Texas Advanced Computing Center, and Big Red II+ at Indiana University. USQCD resources are acquired and operated thanks to funding from the Office of Science of the U.S. Department of Energy. The National Energy Research Scientific Computing Center is a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. The Blue Waters sustained-petascale computing project is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the State of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. This work is also part of the “Lattice QCD on Blue Waters” and “High Energy Physics on Blue Waters” PRAC allocations supported by the National Science Foundation (award numbers 0832315 and 1615006). This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562 Towns et al. (2014). Allocations under the Teragrid and XSEDE programs included resources at the National Institute for Computational Sciences (NICS) at the Oak Ridge National Laboratory Computer Center, The Texas Advanced Computing Center and the National Center for Atmospheric Research, all under NSF teragrid allocation TG-MCA93S002. Computer time at the National Center for Atmospheric Research was provided by NSF MRI Grant CNS-0421498, NSF MRI Grant CNS-0420873, NSF MRI Grant CNS-0420985, NSF sponsorship of the National Center for Atmospheric Research, the University of Colorado, and a grant from the IBM Shared University Research (SUR) program. Computing at Indiana University is supported by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute. This project was supported in part by the URA Visitig Scholar Award 12-S-15 (Y.L.); by the U.S. Department of Energy under grants No. DE-FG02-91ER40628 (C.B.), No. DE-FC02-12ER41879 (C.D.), No. DE-FG02-13ER42001 (A.X.K.), No. DE-SC0015655 (A.X.K., Z.G.), No. DE-SC0010120 (S.G.), No. DE-FG02-91ER40661 (S.G.), No. DE-SC0010113 (Y.M.), No. DE-SC0010005 (E.T.N.), No. DE-FG02-13ER41976 (D.T.); by the U.S. National Science Foundation under grants PHY14-14614 and PHY17-19626 (C.D.), and PHY14-17805 (J.L.); by the MINECO (Spain) under grants FPA2013-47836-C-1-P and FPA2016-78220-C3-3-P (E.G.); by the Junta de Andalucía (Spain) under grant No. FQM-101 (E.G.); by the Fermilab Distinguished Scholars program (A.X.K.); by the German Excellence Initiative and the European Union Seventh Framework Program under grant agreement No. 291763 as well as the European Union’s Marie Curie COFUND program (A.S.K.). Brookhaven National Laboratory is supported by the United States Department of Energy, Office of Science, Office of High Energy Physics, under Contract No. DE-SC0012704. This document was prepared by the Fermilab Lattice and MILC Collaborations using the resources of the Fermi National Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office of Science, HEP User Facility. Fermilab is managed by Fermi Research Alliance, LLC (FRA), acting under Contract No. DE-AC02-07CH11359.

Appendix A $B_{s}\to K\ell\nu$ form factors in SU(2) chiral

perturbation theory

In this appendix, we derive Eq. (22), the SU(2) chiral formula for $B_{s}\to K\ell\nu$.

We start from the next-to-leading order (NLO) SU(3) HMrS$\chi$PT expression for $B_{x}\to P_{xy}$ semileptonic decay. It is expressed as Aubin and Bernard (2007)

[TABLE]

where the subscript $P$ stands for $\parallel$ or $\perp$; $c^{i}_{P}$ are coefficients and the corresponding rescaled quantities in Eq. (55) will be determined by the chiral fits; $\delta f_{P,\mathrm{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\cdot 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_{\parallel}$ and $f_{\perp}$ are

[TABLE]

where $f$ is the decay constant involved; $g$ is the coupling constant; $\Delta_{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., $\Delta_{xy}^{*}=B_{y}^{*}-B_{x}$; and $D_{\text{logs}}$ is the nonanalytic self-energy contribution. The scalar pole was not included in Ref. Aubin and Bernard (2007) as the $0^{+}$ meson is not in leading order HMrS$\chi$PT. It is added here phenomenologically as explained in Sec. IV.4.

For the $B_{s}\to K\ell\nu$ analysis considered here, $x=s^{\prime}$ and $y=u^{\prime}=d^{\prime}=u=d$. Here we use primed quantities to denote the valence quarks and the unprimed for the sea quarks555This is different from the convention used in the main text. For example, in Tables 1 and 2 the prime quantities denote the sea quarks and the unprimed for the valence ones.. All the data generated for $B_{s}\to K\ell\nu$ analysis are partially quenched points, i.e., $m^{\prime}_{s}\neq 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^{\prime}_{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. Bailey et al. (2016), we obtain the SU(2) hard-kaon chiral log terms in $f_{\parallel}$ and $f_{\perp}$ for the $B_{s}\to K\ell\nu$ as

[TABLE]

The summation $\xi$ is over 16 staggered fermion tastes (P, V, T, A, or I); $I_{1}(m)$ is the chiral logarithm defined as

[TABLE]

Meson masses for the 2+1 case in the SU(2) limit ($m_{u}=m_{d}$ and $m_{s}\rightarrow\infty$) are Aubin and Bernard (2003); Bailey et al. (2016)

[TABLE]

The $[V\to A]$ in Eq. (51) stands for terms with subscripts changed from V to A. The hairpin parameters $\delta_{V(A)}^{\prime}$ in Eq. (53b) are listed in Table 11. The $m^{2}_{ij,\xi}$ are defined later in Eq. (56).

We can regroup relevant terms in Eq. (49), drop the $m_{s}^{\prime}$ dependent term due to the SU(2) limit, and write the formula as the following

[TABLE]

We can further write all the expansion parameters in terms of dimensionless ones

[TABLE]

where $\mu$ is the leading-order low-energy constant that relates the tree-level mass of a taste-$\xi$ meson composed of quarks of flavor $i$ and $j$ to the corresponding quark masses

[TABLE]

Here $\Delta_{\xi}$ is the staggered fermion taste splitting. The numerical values of $\mu$ and $\Delta_{\xi}$ are determined by the MILC Collaboration and are shown in Table 11. The average taste splitting in Eq. (55d) is $\bar{\Delta}=\frac{1}{16}\sum_{\xi}\Delta_{\xi}$.

Combining the above information, one arrives at the final NLO form used in the chiral-continuum extrapolation in this work, Eq. (22).

Appendix B Reconstructing the $B_{s}\to K\ell\nu$ form factors

In this appendix, we document the procedure of reconstructing the form factors from the fitting results obtained in Sec. VI.2.

B.1 Reconstructing the form factors as functions of

$z$

The form factors are parametrized in a BCL Bourrely et al. (2009) form with coefficients $b_{i}^{+,0}$ as shown in Eq. (33). The meson masses used in the $z$-parametrization fit are listed in Table 9. The fitted coefficients $b_{i}^{+,0}$ are listed in Table 10. To get the form factors as functions of $z$ and reproduce the left panel result of Fig. 12, one should use Eq. (33) with ${M_{B^{*}}(1^{-})}$ and ${M_{B^{*}}(0^{+})}$ meson mass values in Table 9, and the $b_{i}^{+,0}$ values and the correlation matrix in Table 10.

B.2 Reconstructing the form factors as functions of

$q^{2}$

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. (29), (27), (30), and (32). One can then solve Eq. (29) to get $q^{2}$ in terms of $z$:

[TABLE]

Once we have the form factors as functions of $z$ from Appendix B.1, we can then use Eq. (57) to change the variable to get the $q^{2}$ dependence.

B.3 Dealing with the near zero eigenvalue in the covariance matrix

In Table 10 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

[TABLE]

and the matrix $D$ is a diagonal matrix with diagonal elements $\Sigma$. The correlation matrix is denoted as $R$ and the covariance matrix is denoted as $S$. The relations among $D$, $R$, and $S$ are

[TABLE]

Alternatively, one can use the following relation to directly convert the matrix elements

[TABLE]

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 10. This is because we imposed the kinematical constraint Eq. (35) with $\epsilon=10^{-10}$ in the $z$-parametrization 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_{3}^{0}$, in Eq. (33), express it in terms of the other $b_{i}^{+,0}$ parameters, and then perform the $z$-parametrization fit. This, however, will make the expression Eq. (33b) cumbersomely complicated to handle when performing the fit. In practice, we use the expressions Eq. (33) and perform the fits as described in Sec. VI.2. Whenever one needs to invert the covariance matrix, one simply needs to reduce the size of the matrix by removing one column and one row corresponding to one parameter $b_{r}$. The parameter $b_{r}$ can be any one of the $b_{i}^{+,0}$ parameters. Without loss of generality, let us pick $b_{r}$ to be $b_{K-1}^{0}=b_{3}^{0}$ for our $K=4$ preferred fit. From Eqs. (33 and 35), we can get

[TABLE]

with

[TABLE]

as derived from Eq. (29).

Appendix C $B_{s}\to K\ell\nu$ differential decay rate bin

tables

In this appendix, we present the quantity

[TABLE]

for the $B_{s}\to K\mu\nu$ and $B_{s}\to K\tau\nu$ decays, in bins of $q^{2}$, in Tables 12 and 13. Since we also include the correlations between $q^{2}$ bins in these tables, the results therein can be combined with the corresponding experimental measurements to determine $|V_{ub}|$.

[FIGURE:]

[FIGURE:]