Extraction of bare Form Factors for $\mathrm B_\mathrm s \to \mathrm K \ell \nu$ Decays in non-perturbative HQET

TL;DR

This paper presents a non-perturbative lattice QCD method to extract form factors for B_s to K decays, emphasizing excited state removal and achieving about 2% precision, with future work needed on multi-hadron effects.

Contribution

It introduces a lattice QCD approach using HQET with excited state suppression techniques to accurately determine B_s to K form factors.

Findings

Achieved ~2% precision in form factor extraction.

Controlled excited state contributions effectively.

Identified the need to study multi-hadron contaminations.

Abstract

We discuss the extraction of the ground state $\langle \mathrm{K} ({\bf p})|V_\mu(0)|\mathrm{B} ({\bf 0})\rangle$ matrix elements from Euclidean lattice correlation functions. The emphasis is on the elimination of excited state contributions. Two typical gauge-field ensembles with lattice spacings $0.075, \; 0.05$ fm and pion masses $330,\;270$ MeV are used from the O($a$)- improved CLS 2-flavour simulations and the final state momentum is $|{\bf p}|=0.5\,{\rm GeV}$. The b-quark is treated in HQET including the $1/m_\mathrm{b}$ corrections. Fits to two-point and three-point correlation functions and suitable ratios including summed ratios are used, yielding consistent results with precision of around 2% which is $not$ limited by the $1/m_\mathrm{b}$ corrections but by the dominating static form factors. Excited state contributions are under reasonable control but are the bottleneck towards precision. We do not yet include a specific investigation of multi-hadron contaminations, a gap in the literature which ought to be filled soon.

Tables

Table 1. Table 1 : Overview of the subset of N f = 2 subscript 𝑁 f 2 N_{\rm f}=2 CLS ensembles used. They have T = 2 × L 𝑇 2 𝐿 T=2\times L . Lattice spacings, a 𝑎 a , are taken from Ref.  [ 19 ] and Ref.  [ 20 ] . The pion mass is denoted by m π subscript 𝑚 𝜋 m_{\pi} , and we have m π ​ L ≥ 4 subscript 𝑚 𝜋 𝐿 4 m_{\pi}L\geq 4 such that finite-volume effects are sufficiently (exponentially) suppressed at the level of accuracy we are aiming for.

ensemble	$\beta$	$L/a$	$a$ [fm]	$m_{\pi }$ [MeV]	$m_{\pi }L$
A5	$5.2$	$32$	$0.0749\left(8\right)$	$330$	$4.0$
O7	$5.5$	$64$	$0.0483\left(4\right)$	$270$	$4.2$

Table 2. Table 2 : Details of the measurements. The hopping parameter of the valence strange quark, κ s subscript 𝜅 𝑠 \kappa_{s} , has been determined as described in the text. The number of configurations on which we computed the correlations functions is denoted by N cfg subscript 𝑁 cfg N_{\rm cfg} . The number of smearing iterations applied to the light quark(s) in the K K {\rm K} and B s subscript B s {\textnormal{B}_{\textnormal{s}}} meson is given by I K superscript 𝐼 K I^{\rm K} and I r B s subscript superscript 𝐼 subscript B s 𝑟 I^{\textnormal{B}_{\textnormal{s}}}_{r} , respectively. The integer vector 𝐧 𝐧 {\bf n} and the scalar twist angle ϑ italic-ϑ \vartheta determine the momentum of the K K {\rm K} meson by setting θ s = 𝐧 ​ ϑ superscript 𝜃 s 𝐧 italic-ϑ \theta^{\mathrm{s}}={\bf n}\vartheta in ( 2.19 ), i.e. 𝐩 K = ( 2 ​ π − ϑ ) ​ 𝐧 / L subscript 𝐩 K 2 𝜋 italic-ϑ 𝐧 𝐿 {\bf p}_{\rm K}=(2\pi-\vartheta){\bf n}/L .

ensemble	${\kappa }_s$	$N_{\mathrm{cfg}}$	$I^K$	$(I_1^{{\text{B}}_{\text{s}}},I_2^{{\text{B}}_{\text{s}}},I_3^{{\text{B}}_{\text{s}}})$	$𝐧$	$\vartheta$
A5	$0.13535\left(2\right)$	$1000$	35	$(30,60,155)$	${(1,0,0)}^T$	$+0.2141$
O7	$0.13627\left(1\right)$	$490$	75	$(68,135,338)$	${(1,1,0)}^T$	$-0.3593$

Table 3. Table 3 : Results for the ground-state form factors φ μ ( 0 , 0 ) superscript subscript 𝜑 𝜇 0 0 \varphi_{\mu}^{(0,0)} . The second error in the result for ℛ μ = 1 II subscript superscript ℛ II 𝜇 1 \mathcal{R}^{\,\textrm{II}}_{\mu=1} for ensemble A5 is the systematic error due to the wrapper effect. The last column shows the results of the linear fit to 𝒮 μ , 3 I ​ ( τ ) subscript superscript 𝒮 I 𝜇 3 𝜏 \mathcal{S}^{\textrm{I}}_{\mu,3}(\tau) , see the shaded band in Fig. ( 14 ).

id	$\mu$	Fit		Ratio	Summed ratio
		$t_{\mathrm{B2}}^{\min }\approx 0.32\mathrm{fm}$	$t_{\mathrm{B2}}^{\min }\approx 0.38\mathrm{fm}$	${ℛ}^{\text{II}}$	${ℳ}^{\text{I}}$
A5	0	1.103(12)	1.104(10)	1.129(18)	1.097(17)
A5	1	0.579(9)	0.578(8)	0.579(13)(05)	0.587(11)
O7	0	1.118(16)	1.114(13)	1.113(19)	1.109(13)
O7	1	0.610(13)	0.606(11)	0.602(15)	0.601(10)

Table 4. Table 4 : Matrix elements at 𝒪 ​ ( 1 / m ) 𝒪 1 𝑚 \mathcal{O}\left(1/m\right) : kin and spin contributions.

$\mu$	$k$	${\rho }_{\mu }^k$ A5	${\rho }_{\mu }^k$ O7
0	kin	-0.554(13)	-0.493(12)
1	kin	-0.316(24)	-0.335(21)
0	spin	-0.363(03)	-0.348(02)
1	spin	-0.223(05)	-0.177(03)

Table 5. Table 5 : Overview of the 1 / m 1 𝑚 1/m vector current insertions, their tree-level matching coefficients, and the results extracted from the highest light-quark smearing, r = 3 𝑟 3 r=3 . We use symmetric covariant derivatives ∇ i S subscript superscript ∇ S 𝑖 \nabla^{\textnormal{S}}_{i} .

$\mu$	$j$	${\Gamma }_{\mu }^j$	${\omega }_{\mu ,j}^{\text{tree}}\cdot m_h$	$a{\rho }_{\mu }^j$ A5	$a{\rho }_{\mu }^j$ O7
0	1	${\sum }_l{\gamma }_l\frac{1}{2}\left({\nabla }_l^{\text{S}}-{\overleftarrow{\nabla }}_l^{\text{S}}\right)$	1/2	-0.0182(12)	-0.0241(07)
0	2	${\sum }_l{\gamma }_l\frac{1}{2}\left({\nabla }_l^{\text{S}}+{\overleftarrow{\nabla }}_l^{\text{S}}\right)$	1/2	-0.0434$\left(4\right)$	-0.0289(02)
$i$	1	${\sum }_l\frac{1}{2}\left({\nabla }_l^{\text{S}}-{\overleftarrow{\nabla }}_l^{\text{S}}\right){\gamma }_l{\gamma }_i$	1/2	0.3508(27)	0.2852(15)
$i$	2	$\frac{1}{2}\left({\nabla }_i^{\text{S}}-{\overleftarrow{\nabla }}_i^{\text{S}}\right)$	-1	0.0759(17)	0.0709(11)
$i$	3	${\sum }_l\frac{1}{2}\left({\nabla }_l^{\text{S}}+{\overleftarrow{\nabla }}_l^{\text{S}}\right){\gamma }_l{\gamma }_i$	1/2	0.1737(16)	0.1132(08)
$i$	4	$\frac{1}{2}\left({\nabla }_i^{\text{S}}+{\overleftarrow{\nabla }}_i^{\text{S}}\right)$	-1	0.1737(16)	0.1133(08)

Table 6. Table 6 : Results for the ground-state form factors φ μ ( 0 ) superscript subscript 𝜑 𝜇 0 \varphi_{\mu}^{(0)} using the fit with N B s = 2 subscript 𝑁 subscript B s 2 N_{\textnormal{B}_{\textnormal{s}}}=2 .

id	$\mu$	Fit value
A5	0	1.099(12)
A5	1	0.580(9)
O7	0	1.100(13)
O7	1	0.606(9)

Equations

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

$$f_0(q^2)$$

$$(2m_{{\text{B}}_{\text{s}}}{)}^{-1/2}⟨\mathrm{K}({\mathbf{p}}_{\mathrm{K}})|{\hat{V}}_{\mathrm{M}}^0(0)|{\text{B}}_{\text{s}}(0)⟩$$

$$(2m_{{\text{B}}_{\text{s}}}{)}^{-1/2}⟨\mathrm{K}({\mathbf{p}}_{\mathrm{K}})|{\hat{V}}_{\mathrm{M}}^k(0)|{\text{B}}_{\text{s}}(0)⟩$$

$$\begin{array}{rl}{\mathcal{C}}_{\mu }^{{\text{B}}_{\text{s}}\to \mathrm{K}}(t_{\mathrm{K}},t_{{\text{B}}_{\text{s}}};{\mathbf{p}}_{\mathrm{K}})=& \frac{a^9}{L^3}\sum _{{\mathbf{x}}_{\mathrm{f}},{\mathbf{x}}_{\mathrm{v}},{\mathbf{x}}_{\mathrm{i}}}{e}^{-i{\mathbf{p}}_{\mathrm{K}}({\mathbf{x}}_{\mathrm{f}}-{\mathbf{x}}_{\mathrm{v}})}\times ⟨{\overline{\mathcal{O}}}_{\mathrm{us}}(x_{\mathrm{f}})V_{\mu }(x_{\mathrm{v}}){\mathcal{O}}_{\mathrm{bs}}(x_{\mathrm{i}})⟩,\\ & t_{\mathrm{K}}=(x_{\mathrm{f}}-x_{\mathrm{v}}{)}_0>0,t_{{\text{B}}_{\text{s}}}=(x_{\mathrm{v}}-x_{\mathrm{i}}{)}_0>0,\end{array}$$

$$V_{\mu }(x)={\overline{\psi }}_{\mathrm{u}}(x){\gamma }_{\mu }{\psi }_{\mathrm{b}}(x).$$

$${\mathcal{C}}^{\mathrm{K}}(t;{\mathbf{p}}_{\mathrm{K}})$$

$${\mathcal{C}}^{{\text{B}}_{\text{s}}}(t)$$

$$⟨\mathrm{K}(\mathbf{p})|{\hat{V}}_{\mathrm{M}}^0(0)|{\text{B}}_{\text{s}}(0)⟩$$

$$⟨\mathrm{K}(\mathbf{p})|{\hat{V}}_{\mathrm{M}}^k(0)|{\text{B}}_{\text{s}}(0)⟩$$

$$\log ({\mathcal{C}}_{\mu }^{{\text{B}}_{\text{s}}\to \mathrm{K}})=\log Z_{V_{\mu }}^{\mathrm{HQET}}+\log ([{\mathcal{C}}_{\mu }^{{\text{B}}_{\text{s}}\to \mathrm{K}}{]}_{\mathrm{stat}})+\sum _k{\omega }_k{\rho }_{\mu }^k+{\delta }_{\mathcal{O}}+\mathcal{O}\left(1/m_{\mathrm{b}}^2\right),$$

$${\rho }_{\mu }^k=\frac{[{\mathcal{C}}_{\mu }^{{\text{B}}_{\text{s}}\to \mathrm{K}}{]}_k}{[{\mathcal{C}}_{\mu }^{{\text{B}}_{\text{s}}\to \mathrm{K}}{]}_{\mathrm{stat}}}.$$

$$V_0^{\text{stat}}$$

$$V_k^{\text{stat}}$$

$$c_{{\text{V}}_k}=0.0518g_0^2,c_{{\text{V}}_0}=0.0380g_0^2$$

$${\mathcal{O}}_{\mathrm{bs},r}(x)$$

$${\overline{\mathcal{O}}}_{\mathrm{bs},r}(x)$$

$${\psi }_{\mathrm{s},r}^{\mathrm{sm}}=(1-{\kappa }_{\mathrm{G}}{a}^2\Delta {)}^{I_r^{{\text{B}}_{\text{s}}}}{\psi }_s.$$

$${\mathcal{O}}_{\mathrm{su}}(x_{\mathrm{f}})$$

$${\overline{\mathcal{O}}}_{\mathrm{su}}(x_{\mathrm{f}})$$

$$q^2=(m_{{\text{B}}_{\text{s}}}^{\mathrm{phys}}{)}^2+(m_{\mathrm{K}}^{\mathrm{phys}}{)}^2-2m_{{\text{B}}_{\text{s}}}^{\mathrm{phys}}\sqrt{(m_{\mathrm{K}}^{\mathrm{phys}}{)}^2+{\mathbf{p}}_{\mathrm{K}}^2}$$

$${\psi }_{\mathrm{s}}(x_0,\mathbf{x}+\hat{\mathbf{k}}L)$$

$${\psi }_{\mathrm{b}}(x_0,\mathbf{x}+\hat{\mathbf{k}}L)$$

$$(p_{\mathrm{K}}{)}_k=(2\pi n_k-{\theta }_k^{\mathrm{s}})/L,$$

$$\delta q^2\approx 10\mathrm{GeV}\cdot \delta m_B+5GeV\cdot \delta m_K+2Ge{V}^2\cdot \delta \vartheta +5Ge{V}^2\cdot \frac{\delta a(\beta )}{a(\beta )}.$$

$${\mathcal{C}}^{\mathrm{K}}(t_{\mathrm{K}})=$$

$${\mathcal{C}}_{rs}^{{\text{B}}_{\text{s}}}(t_{{\text{B}}_{\text{s}}})=$$

$$\begin{array}{rl}{\mathcal{C}}_{\mu ,r}^{{\text{B}}_{\text{s}}\to \mathrm{K}}(t_{\mathrm{K}},t_{{\text{B}}_{\text{s}}})=& \sum _{n=0}^{\infty }\sum _{m=0}^{\infty }{\kappa }^{(m)}{\phi }_{\mu }^{(m,n)}{\beta }_r^{(n)}{e}^{-E_{\mathrm{K}}^{(m)}{t}_{\mathrm{K}}}{e}^{-E_{{\text{B}}_{\text{s}}}^{(n)}{t}_{{\text{B}}_{\text{s}}}}\cong \\ & \sum _{n=0}^{N_{\mathrm{K}}-1}\sum _{m=0}^{N_{{\text{B}}_{\text{s}}}-1}{\kappa }^{(m)}{\phi }_{\mu }^{(m,n)}{\beta }_r^{(n)}{e}^{-E_{\mathrm{K}}^{(m)}{t}_{\mathrm{K}}}{e}^{-E_{{\text{B}}_{\text{s}}}^{(n)}{t}_{{\text{B}}_{\text{s}}}},\end{array}$$

$$⟨\mathrm{K}(\mathbf{p})|{\hat{V}}_{\mathrm{M}}^0(0)|{\text{B}}_{\text{s}}(0)⟩=$$

$${\mathcal{C}}^{\mathrm{K}}(t)\cong ({\kappa }^{(0)}{)}^2(e^{-E_{\mathrm{K}}^{(0)}t}+e^{-E_{\mathrm{K}}^{(0)}(T-t)}).$$

$${\mathcal{C}}^{{\text{B}}_{\text{s}}}(t)v^{(n)}(t,t_0)={\lambda }^{(n)}(t,t_0){\mathcal{C}}^{{\text{B}}_{\text{s}}}(t_0)v^{(n)}(t,t_0),n=1,\dots ,N_{{\text{B}}_{\text{s}}},$$

$$E_{{\textnormal{B}_{\textnormal{s}}},{\rm eff}}^{(n)}(t)=-\tilde{\partial}_{t}\log\big{(}\lambda^{(n)}(t,t_{0}))\big{)}=\frac{1}{2a}\log\big{(}\lambda^{(n)}(t-a,t_{0})/\lambda^{(n)}(t+a,t_{0})\big{)}.$$

$${\mathcal{C}}_{rr}^{{\text{B}}_{\text{s}}}(t)\cong \sum _{n=0}^{N_{{\text{B}}_{\text{s}}}-1}({\beta }_r^{(n)}{)}^2{e}^{-E_{{\text{B}}_{\text{s}}}^{(n)}t}$$

$${\mathcal{C}}_{\mu ,r}^{{\text{B}}_{\text{s}}\to \mathrm{K}}(t_{\mathrm{K}},t_{{\text{B}}_{\text{s}}})\cong$$

$${\kappa }^{(0)}{\xi }_{\mu ,r}{e}^{-E_{{\mathrm{B}}^{\ast }}{t}_{{\text{B}}_{\text{s}}}}{e}^{-E_{\mathrm{K}}^{(0)}(T-t_{{\text{B}}_{\text{s}}}-t_{\mathrm{K}})},$$

$${\mathcal{C}}_{\mu ,r}^{{\text{B}}_{\text{s}}\to \mathrm{K}}(t_{\mathrm{K}})\cong B_{\mu ,r}{e}^{-E{t}_{\mathrm{K}}}+C_{\mu ,r}{e}^{+E{t}_{\mathrm{K}}},$$

Full text

1. 1 Introduction

2. 2 Setup for the non-perturbative evaluation

3. 2.1 Form factors from Euclidean correlation functions

4. 2.2 Lattice Setup

5. 3 Static correlation functions

6. 3.1 Kaon correlation functions

7. 3.2 ${\textnormal{B}_{\textnormal{s}}}$ correlation functions

8. 3.2.1 Energies from the GEVP

9. 3.2.2 Amplitudes

10. 3.3 Three-point static correlation functions

11. 4 Matrix elements at static order from a combined fit

12. 4.1 Fit ranges

13. 4.2 Fit quality and results

14. 5 Matrix elements at static order from ratios

15. 5.1 Ordinary ratios

16. 5.2 Summed ratios

17. 5.3 GEVP method

18. 6 Matrix elements at $1/m$ order

19. 6.1 Numerical results for the kin and spin insertions at $1/m$ order

20. 6.2 Current insertions at $1/m$ order

21. 7 Discussion

22. A Computation of the correlation functions

23. A.1 Improvement and $1/m$ terms

24. B Alternative static fit with $N_{\rm K}=1,N_{\textnormal{B}_{\textnormal{s}}}=2$

25. C Note on the convergence of the static summed ratio

1 Introduction

Decays of b-quarks in the form $\mathrm{B}\to\pi\ell\nu$ and $\mathrm{B}_{\mathrm{s}}\to\mathrm{K}\ell\nu$ are very relevant in constraining the Standard Model of particle physics and hence also in the search for deviations from it. Lattice QCD is the method of choice for predicting the necessary form factors [1, 2, 3, 4, 5, 6, 7, 8, 9]. While computations are progresssing significantly, there are still few groups carrying them out. Crosschecks by independent (lattice) methods are largely absent. In Ref. [10] we have started an investigation which takes significantly different avenues in almost every choice that can be made. We use (improved) Wilson fermions, treat the b-quark in non-perturbative HQET and compute three-point functions for all time separations rather than fixing its total time-span (called $\tau={t_{\rm K}}+{t_{{\textnormal{B}_{\textnormal{s}}}}}$ later) to one or few values. We opt for non-perturbative HQET because a complete and practical method is known to non-perturbatively renormalize the theory and match it to QCD [11, 12, 13, 14, 15]. When completed, this will provide an independent crosscheck.

Once these choices are made, one needs to carry out the basic steps as written in Ref. [10]:

• a)

obtain the ground state matrix elements $\langle{\rm K}|V^{\mu}(0)|{\textnormal{B}_{\textnormal{s}}}\rangle$ that mediate the transition,

• b)

renormalize the currents (and thus matrix elements) and, if an effective theory is used, relate them to QCD (“matching”),

• c)

take the continuum limit of the matrix elements,

• d)

extrapolate to the quark masses realized in Nature,

• e)

map out the $q^{2}$ dependence.

In the previous work we saw that the continuum limit is smooth in the static approximation, which is expected to yield the by far dominating part of the full result. In this paper we discuss the step a) in detail, including the $1/m$ terms. It is not at all an easy enterprise. As discussed at length below, the issue is that ground state matrix elements are obtained at large Euclidean time separations. At time separations of around $\frac{1}{2}$ fm — a typical QCD length scale — the low-lying states start dominating the two- and three-point functions. Unfortunately the statistical quality of the Monte Carlo estimates is typically deteriorating fast once 2 fm in total (or 1 fm for $K$ and ${\textnormal{B}_{\textnormal{s}}}$ each) are reached. It is thus a matter of the details of methods and QCD dynamics whether there is a window to determine the desired form factors at the percent level.

In other words, control over the extraction of the bare form factors is crucial in a proper determination of the decay rates. In this publication we thus explain our methods and their limitations in detail.

2 Setup for the non-perturbative evaluation

2.1 Form factors from Euclidean correlation functions

To leading order in the weak interactions, the ${\textnormal{B}_{\textnormal{s}}}\to{\rm K}\ell\nu$ decay rate depends on two form factors, $h_{\parallel}(E_{\rm K})$ and $h_{\perp}(E_{\rm K})$. which are related to the commonly used

[TABLE]

The kinematics is explained in Ref. [10], and in the rest-frame of the ${\textnormal{B}_{\textnormal{s}}}$, we have

[TABLE]

where $(E_{\rm K},{\bf p}_{\rm K})$ is the four-momentum of the final state and $|{\textnormal{B}_{\textnormal{s}}}({\bf 0})\rangle$ is the initial state with the ${\textnormal{B}_{\textnormal{s}}}$ at rest, ${\bf p_{\textnormal{B}_{\textnormal{s}}}}={\bf 0}$. The above hadron-to-hadron matrix elements are in Minkowski space, as indicated by the subscript M. They can be obtained from the Euclidean three-point function (fig. 1), where the spatial volume is a $L^{3}$ torus,

[TABLE]

with suitable interpolating fields ${\cal O}_{qq^{\prime}},\overline{{\cal O}}_{qq^{\prime}}$ (see (2.14),(2.16)) and

[TABLE]

Removing the overlaps $\langle{\textnormal{B}_{\textnormal{s}}}({\bf 0})|{\cal O}_{\rm bs}|0\rangle$, $\langle 0|\overline{{\cal O}}_{\rm us}|{\rm K}({\bf p}_{\rm K})\rangle$ (amputation), requires also the two-point functions

[TABLE]

The explicit relation between the Euclidean correlation functions and the desired matrix elements is (see e.g. [16])

[TABLE]

A few comments are in order

We chose to define the correlation functions on a lattice. Replacing $a^{3}\sum_{\bf x}\to\int\mathrm{d}^{3}{\bf x}$ yields the continuum expressions.

• -

The derivation of Eq. (LABEL:eqn:me1) (on a lattice) is based just on the existence of a transfer matrix with standard properties and on interpolating fields which are local in time; non-locality (smearing) in space is allowed.

• -

Of course, on a finite lattice the limit $t\to\infty$ has to be replaced by a proper procedure. Controlling it is the main topic of this paper. As a first, straightforward, procedure we note that it would be sufficient to have $t$ (as well as the spatial box length $L$) large compared to the typical QCD scales $\frac{1}{2}{\rm fm},\;1/m_{\pi}$.

• -

In order to obtain the physical form factors, the current $V_{\mu}$ has to be properly renormalized and the continuum limit has to be taken.

So far we have spoken of standard relativistic (lattice) QCD. The non-perturbative HQET expressions including $\mathcal{O}\left(1/m\right)$ terms are obtained by the following generic expansion of correlation functions [11, 17, 18]

[TABLE]

where $[{\cal C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}]_{\mathrm{stat}}$ is the correlation function in the static approximation and

[TABLE]

The "stat" terms refer to lowest order (static) HQET, while $k\in\{\rm kin,spin\}$ are the terms due to the $1/m$ corrections to the action, and the remaining terms in the sum over $k$ are $1/m$ corrections to the currents, see Sect. 6.2. In $\delta_{\cal O}$, we collect all terms which arise from the renormalization and $1/m$ corrections to the interpolating fields. Their form is irrelevant since they cancel in the ratios of interest Eq. (LABEL:eqn:me1). Apart from $\delta_{\cal O}$, all divergences strictly cancel within Eq. (2.10) iff the parameters $\log Z^{\mathrm{HQET}}_{V_{\mu}}$ and $\omega_{k}$ are determined non-perturbatively. The expansion of $\log({\cal C}^{{\textnormal{B}_{\textnormal{s}}}})$ reads the same apart from the current-terms which are absent. Finally, functions of correlation functions such as the ratio Eq. (LABEL:eqn:me1) are trivially obtained from the expansion of the correlation functions.

For completeness we remind the reader of the rule of the game: the expansion is “derived” by considering the parameters $\omega_{k}$ as $\mathcal{O}\left(1/m\right)$ and expanding in $1/m$, despite the fact that individual pieces may be highly divergent. The logarithm is taken in Eq. (2.10) since then all expressions are automatically linear in the parameters and higher order terms in the parameters (and thus $1/m$) are directly avoided. Those would otherwise have to be dropped to preserve renormalizability.

2.2 Lattice Setup

For the lattice computation of the correlation functions we use gauge field configurations generated by the Coordinated Lattice Simulations (CLS) effort [21] with two degenerate flavors of O$(a)$-improved Wilson fermions [22] and Wilson gauge action[23]. In the following, we present only results from the two representative ensembles, labeled A5 and O7, see Table 1. The lattice spacings were determined from the Kaon decay constant $f_{\rm K}$ at physical pion mass in Ref. [21] and updated in Ref. [24]. They do not depend on the quark mass, i.e. we use a mass-independent renormalisation scheme.

The natural choice for the mass of the valence up quark is equal to the mass of the degenerate up and down sea quarks, $m_{\rm l}$, and thus fixed by the gauge field ensembles. For the mass of the spectator (valence) strange quark we are free, however, to choose any smooth function $m_{{\rm s}}(m_{\rm l})$ passing through the physical point. We follow strategy 1 of Ref. [21] and define $m_{{\rm s}}(m_{\rm l})$ by fixing the squared ratio between the Kaon mass and the Kaon decay constant to its physical value $R^{\rm phys}_{\mathrm{K}}\equiv(m^{\rm phys}_{\rm K}/f^{\rm phys}_{\rm K})^{2}$ with $m^{\rm phys}_{\rm K}=494.2\,{\rm MeV}$ and $f^{\rm phys}_{\rm K}=155\,{\rm MeV}$. We expect that this will lead to a flat extrapolation to the physical value of $m_{\rm l}$, the “physical point”. At the simulated (unphysical) values of $m_{\rm l}$, however, the resulting Kaon mass is $m_{\rm K}(m_{\rm l})=f_{\rm K}(m_{\rm l})\cdot(R^{\rm phys}_{\mathrm{K}})^{1/2}$ and only approximately equal to its physical value (see e.g. Fig. 3).

In the static correlation functions on the r.h.s. of Eq. (2.10), the heavy quark field $\psi_{\mathrm{b}}$ is replaced by a static quark $\psi_{\mathrm{h}}=\gamma_{0}\psi_{\mathrm{h}}$ with the HYP2 discretizations [25, 26]. At NLO also the $k\in\{\rm kin,spin\}$ correlation functions are needed, see Appendix A.1. They are static correlation functions with an extra insertion of ${\cal O}^{\mathrm{kin}}$ or ${\cal O}^{\mathrm{spin}}$ from the HQET action [11]. Finally there are the extra $1/m$ terms in the current $V_{\mu}$, see Sec. 6.2 and Appendix A.1. We implement O($a$) improvement perturbatively by writing the vector current as

[TABLE]

with [27]

[TABLE]

when we evaluate the static matrix elements, but we use $c_{\textnormal{V}_{k}}=c_{\textnormal{V}_{0}}=0$ when $V_{\mu}^{\rm stat}$ enters ratios at order $1/m$ such as Eq. (2.11).

For the interpolating field of the ${\textnormal{B}_{\textnormal{s}}}$ meson we use three different quark bilinears ($r=1,2,3$)

[TABLE]

The smeared light quark fields are constructed by $I^{\textnormal{B}_{\textnormal{s}}}_{r}$ iterations of Gaussian smearing [28, 29]

[TABLE]

where the gauge links in the covariant Laplacian $\Delta$ are defined by 3 iterations of (spatial) APE smearing [30] and $\kappa_{\mathrm{G}}=0.1$.

The interpolating field for the Kaon,

[TABLE]

is constructed from two smeared quark fields, $\overline{\psi}^{\mathrm{sm}}_{\mathrm{s}}$ and $\psi^{\mathrm{sm}}_{\mathrm{u}}$, each with $I^{\rm K}$ iterations of Gaussian smearing.

Integrating out the Grassmann-valued quark fields yields the correlation functions in terms of traces of products of quark propagators and other factors such as smearing operators. In order to make full use of translation invariance, we evaluate the traces by a stochastic estimator which we represent by a single random $U(1)$ vector with support on time slice $(x_{\mathrm{f}})_{0}$. All values of $(x_{\mathrm{f}})_{0}$ are averaged over (“full time dilution”). More details can be found in Appendix A.

The momentum transfer $q^{2}$ is computed from the continuum relation

[TABLE]

with ${\bf p}_{\rm K}$ in the rest frame of the ${\textnormal{B}_{\textnormal{s}}}$ and physical values of the masses (using $m^{\rm phys}_{\textnormal{B}_{\textnormal{s}}}=5366.77\,{\rm MeV}$). In order to have $q^{2}$ fixed to the same value on all ensembles, we need to adjust ${\bf p}_{\rm K}$ accordingly. This is achieved by using flavor-twisted boundary conditions for the (valence) $\mathrm{s}$ and $\mathrm{b}$ quarks, i.e. instead of periodic boundary conditions, we impose

[TABLE]

for the unit vectors $\hat{{\bf k}}$ in $k$-direction. In this way, the momentum of the ${\rm K}$ meson becomes

[TABLE]

while the ${\textnormal{B}_{\textnormal{s}}}$ meson is kept at rest (${\bf p}_{\textnormal{B}_{\textnormal{s}}}\sim\theta^{\mathrm{b}}-\theta^{\mathrm{s}}=0$) by choosing $\theta^{\mathrm{b}}=\theta^{\mathrm{s}}$. We fix the value of $|{\bf p}_{\rm K}|=0.535\,{\rm GeV}$, and hence $q^{2}=21.23\,{\rm GeV}^{2}$. 111 Note that small changes of the input parameters lead to

$\delta q^{2}\approx 10\,{\rm GeV}\cdot\delta m_{B}+5\,GeV\cdot\delta m_{K}+2\,GeV^{2}\cdot\delta\vartheta+5\,GeV^{2}\cdot\displaystyle\frac{\delta a(\beta)}{a(\beta)}\,.$

The particular value for ${\bf p}_{\rm K}$ was chosen on the “N6” lattice which we do not discuss here, but which will be used in the physics analysis. On N6 we have $L/a=48$ and $\beta=5.5$ such that ${\bf n}=(1,0,0)^{T}$ and $\theta^{\mathrm{s}}={\bf 0}$ lead to $|{\bf p}_{\rm K}|=0.535\,{\rm GeV}$.

3 Static correlation functions

For ${t_{\rm K}}$, ${t_{{\textnormal{B}_{\textnormal{s}}}}}\gg a$, where any possible violations of positivity are exponentially damped222In the pure Wilson regularization and with static quarks, positivity of the lattice transfer matrix is exact. With the O$(a)$-improvement term this may be different, but universality means that non-positive contributions may appear only at distances close to the lattice spacing. For a concrete example of such contributions we refer to Ref. [31]., and assuming infinite ${\rm T}$ for the moment, we can decompose the Euclidean correlation functions as (in the following we keep the ${\bf p}$ dependence of $\mathcal{C}^{\rm K}$ and $\mathcal{C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}$ implicit):

[TABLE]

where the indices $r,s$ label the smearing levels used for the ${\textnormal{B}_{\textnormal{s}}}$ meson, while the indices $m,n$ label the Kaon and ${\textnormal{B}_{\textnormal{s}}}$ meson energy levels respectively. A rather small number of states leads to a good precision if the time separations $t_{\rm K},\,t_{\textnormal{B}_{\textnormal{s}}}$ are sufficiently large. In our fits we will use $N_{\rm K}=1$ and $N_{\textnormal{B}_{\textnormal{s}}}\leq 3$.

By comparing the Euclidean representations of the correlation functions with Eq. (LABEL:eqn:me1) we see that the bare form factors of interest are given by $\varphi^{(0,0)}_{\mu}$ from Eq. (3.22), e.g.

[TABLE]

In the following we describe two methods to determine them:

By directly fitting333All the fits in this paper are uncorrelated, i.e. we use a diagonal covariance matrix. the correlation functions in a window of time extents where one can reliably limit the influence of excited states, statistical noise and other systematic effects (such as the finite-$T$ extent of the available periodic lattices). This will be the focus of Section 4. 2. 2.

By forming suitable ratios of correlation functions (like the one in Eq. (LABEL:eqn:me1)) where at large enough $t$ the dependence on all or most of the additional parameters cancels out. This will be the topic of Section 5.

As we shall see, for both of these methods it is beneficial to first obtain a good determination of the two-point function parameters, as well as to determine regions in the three-point functions where one is safe from the finite-$T$ effects. These prerequisite steps are described in the following subsections. First we describe our methods to extract $\kappa^{(0)}$ and $E_{\rm K}^{(0)}$ from Eq. (3.20).

3.1 Kaon correlation functions

While Eq. (3.20) is true in the limit of infinite $T$, in practice we have to take into account that $T$ is finite. With our periodic boundary conditions in time, the ground-state contribution reads

[TABLE]

Since the signal-to-noise problem in the Kaon sector is very mild, we average $\mathcal{C}^{\rm K}(t)$ and $\mathcal{C}^{\rm K}(T-t)$ and always fit until the middle of the lattice, $t_{\rm K2}^{\rm max}=T/2$. For the start of the fit range, we follow Refs. [19, 32] and use the (rather conservative) criterion that $t^{{\rm K}2}_{\rm min}$ is the first $t$ for which the excited-state contribution, estimated by a two-state fit, is less than 25 % of the statistical error of the effective mass $E_{{\rm K},{\rm eff}}(t)=-\frac{1}{a}\log[\mathcal{C}^{\rm K}(t+a)/\mathcal{C}^{\rm K}(t)]$. In practice, we obtain $t_{\rm K2}^{\rm min}\approx 1.3-1.4\,\textrm{fm}$. In this range, we fit the correlator to Eq. (3.24) to get the ground-state energy and amplitude. Note that the gap between the ground state and excited states in the Kaon sector is large, approximately $1\,{\rm GeV}\approx 5\,{\rm fm}^{-1}$.

To check the behaviour of the signal-to-noise ratio and the dispersion relation, we calculated the Kaon two-point function at two values of momentum. The effective mass plots for both momenta are presented together with the corresponding fits in Fig. 2. At higher momentum the signal-to-noise problem becomes significantly more severe. The dispersion relation, also including a measurement for ${\bf p}=0$, is presented in Fig. 3. The continuum dispersion relation describes the data very well even at the coarsest lattice spacing, showing that the cutoff effects in $E_{\rm K}/m_{\rm K}$ are below our statistical precision.

3.2 ${\textnormal{B}_{\textnormal{s}}}$ correlation functions

For the two-point heavy-light (${\textnormal{B}_{\textnormal{s}}}$) correlator, we have three different smearings (in addition to the unsmeared correlator, which is not used due to its large contamination from highly excited states). With the diagonal and the off-diagonal terms, we have six independent correlators in a symmetric $\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}_{rs}$ matrix. The determination of the parameters of the ${\textnormal{B}_{\textnormal{s}}}$ correlation functions is divided into several steps, which are described in the following subsections.

3.2.1 Energies from the GEVP

We first determine the energies by solving the generalized eigenvalue problem (GEVP) on $\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}_{rs}$ [33]. We use all $N_{{\textnormal{B}_{\textnormal{s}}}}=3$ interpolating fields. The GEVP is defined as

[TABLE]

and yields approximations to the lowest $N_{{\textnormal{B}_{\textnormal{s}}}}$ energy levels via

[TABLE]

We use $t_{0}=\lceil t/2\rceil$ since this asymptotically accelerates the convergence (with $t$) to the energies $E_{\textnormal{B}_{\textnormal{s}}}^{(n)}$ [34]. For each $n$ separately we then find the GEVP estimate of the energy $E_{{\textnormal{B}_{\textnormal{s}}}}^{(n)}$ by a plateau fit (a weighted average) of $E_{{\textnormal{B}_{\textnormal{s}}},{\rm eff}}^{(n)}(t)$. The upper end, $t_{\textrm{max}}$, of the plateau is chosen with a noise criterion: only $E_{{\textnormal{B}_{\textnormal{s}}},{\rm eff}}^{(n)}(t)$ with a relative error less than 20% (50% for the excited states) are included in the average. The start of the plateau, $t_{\rm min}$, is chosen by requiring $r(t_{\rm min})=\frac{|E^{(n)}_{\textnormal{plat}}(t_{\rm min})-E^{(n)}_{\textnormal{plat}}(t_{\rm min}-\delta t)|}{\sqrt{\sigma^{2}(t_{\rm min})+\sigma^{2}(t_{\rm min}-\delta t)}}\leq 3,$ where $E^{(n)}_{\textnormal{plat}}(t_{\rm min})$ is the result of the plateau fit starting at $t_{\rm min}$ and $\delta t\approx 2/(E_{\textnormal{B}_{\textnormal{s}}}^{(N)}-E_{\textnormal{B}_{\textnormal{s}}}^{(0)})$, which we approximate by choosing $\delta t=0.3$ fm. For the ground state this criterion coincides with that of Ref. [35], ensuring that the statistical error dominates over the systematic one. Examples of $E_{{\textnormal{B}_{\textnormal{s}}},{\rm eff}}^{(0)}(t)$ are presented in Fig. 4 (together with $E_{{\rm K},{\rm eff}}(t)$ for comparison). Note that a good estimate of the ground-state energy will be important in Sec. 5.

For the excited states one should in principle use a higher value of $\delta t\approx 2/(E_{\textnormal{B}_{\textnormal{s}}}^{(N)}-E_{\textnormal{B}_{\textnormal{s}}}^{(n)})$ but this is not possible with the current precision of the data, as one ends up in regions dominated by noise. However, as opposed to the ground state, for the excited-state energies we only need reasonable rough estimates, as any residual systematic error will be eliminated in the subsequent fitting steps as described in Sec. 4.

3.2.2 Amplitudes

Having obtained the energies, we determine the amplitudes $\beta_{r}^{(n)}$. We start with the diagonal elements of ${\cal C}^{{\textnormal{B}_{\textnormal{s}}}}$,

[TABLE]

and estimate the squared amplitudes through a linear fit. These values are used as initial guess for a non-linear fit for the amplitudes $\beta^{(n)}_{r}$ to all the correlators, including the off-diagonal ones. In this way also the relative signs of the amplitudes are obtained.

Then, we try the full fit for both amplitudes and energies whenever possible. On certain ensembles the full fit becomes unstable with $N_{\textnormal{B}_{\textnormal{s}}}=3$, in this case we fall back to the fit with fixed energies.

The upper end of the fit range, $t_{\rm B2}^{\rm max}$, is common for all the six correlators and is determined by the relative noise criterion, i.e. as the last point at which the relative error is smaller than 2.5% for all of the $\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}_{rs}$. For our ensembles this gives $2.1-2.4\;\textrm{fm}$. Also the lower end of the fit range, $t_{\rm B2}^{\rm min}$, is the same for all $\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}_{rs}$. It is a tunable parameter that is common for the determination of the amplitudes and for the combined fit (cf. Sec. 4).

3.3 Three-point static correlation functions

Due to the finite time extent of the lattice, the three-point function also receives contributions from particles propagating (“wrapping”) around the torus in time direction and is the sum of the contributions from the two diagrams shown in Fig. 5. At large enough time separations, the contributions of the corresponding lowest states can be written as

[TABLE]

where $\xi_{\mu,r}=\langle 0|V_{\mu}|{\rm B}^{*}\rangle\langle{\rm B}^{*}|P_{\rm hl}|{\rm K}\rangle$ contains unknown matrix elements of the state ${\rm B}^{*}$, the lightest heavy-light state contributing to the “wrapper” diagram. At static order the energy $E_{{\rm B}^{*}}$ is equal to $E_{\rm{\textnormal{B}_{\textnormal{s}}}}^{(0)}$.

Instead of introducing extra fit parameters $\xi_{\mu,r}$ (and possibly others, if excited states need to be included), we restrict ourselves to time separations where the wrapper contributions are negligible.

For this purpose, at every fixed $t_{\textnormal{B}_{\textnormal{s}}}$ we fit the three-point function to the form

[TABLE]

with $B_{\mu,r}$ and $C_{\mu,r}$ being linear fit parameters (which one can express in terms of the amplitudes and matrix elements from Eq. (3.28)) and $E$ being a non-linear fit parameter (although one could in principle set it to $E_{\rm K}$ extracted from the two-point function). The fit is done in a region where ground-state dominance in the Kaon sector is expected (starting at approximately $t_{\rm K3}^{\rm min}=0.8$ fm). Fig. 6 shows an example of such a fit at fixed $t_{\textnormal{B}_{\textnormal{s}}}$.

For every $t_{\textnormal{B}_{\textnormal{s}}}$, $\mu$, and $r$ separately, we then find $t_{\rm K3}^{\rm max,wr}$ as the last $t_{\rm K}$ that fulfills the condition:

[TABLE]

with $c_{\rm wr}=0.25$.

The times $t_{\rm K3,\mu}^{\rm max,wr}(t_{\textnormal{B}_{\textnormal{s}}})$ are shown by the red lines in Fig. 7 for the two representative ensembles A5 and O7 (together with other constraints for the fits as discussed in Sect. 4). We observe that the wrapper criterion yields considerably different $t_{\rm K3}^{\rm max,wr}$ for $\mu=0$ and $\mu=1$ (with lower $t_{\rm K3}^{\rm max,wr}$ for $\mu=1$). On the other hand, the different smearing levels give very similar limits, but are still all treated separately.

4 Matrix elements at static order from a combined fit

From the two-point function fits, we have obtained (estimates of) the energies and amplitudes, $\kappa^{(m)}$, $\beta_{r}^{(n)}$. Taking them as input values, we can now fit Eq. (3.22) for each $\mu=0,1$ with $\varphi^{(m,n)}_{\mu}$, i.e. the form factors, as $N_{\rm K}\times N_{\textnormal{B}_{\textnormal{s}}}$ linear parameters. These could in principle serve as the final result. However, for better stability, accuracy, and control of the systematic errors, we take all these results only as initial conditions for a combined simultaneous fit to the two-point and the three-point functions. In this way, we can use the information on the energies and amplitudes also from the three-point functions, which contain many data points. To estimate the errors we use [36, 37] but with the derivatives with respect to the fit parameters calculated analytically.

The numerical results presented in this section are for the fit with $N_{\rm K}=1$ and $N_{\textnormal{B}_{\textnormal{s}}}=3$, corresponding to a total of 20 fit parameters444$N_{\rm K}+N_{\textnormal{B}_{\textnormal{s}}}$ energies, $N_{\rm K}^{2}+N_{\textnormal{B}_{\textnormal{s}}}^{2}$ amplitudes, and $2\times N_{\rm K}\times N_{\textnormal{B}_{\textnormal{s}}}$ form factors. We find that this choice gives the best stability with respect to the changes of the fit ranges, and therefore we use it to extract the values of $\varphi^{(0)}_{\mu}$. An alternative fit with $N_{\rm K}=1$, $N_{\textnormal{B}_{\textnormal{s}}}=2$ is discussed in Appendix B. The fits with both $N_{\rm K}>1$ and $N_{\textnormal{B}_{\textnormal{s}}}>1$ can easily become very unstable without some form of prior assumptions about the values of the excited-state form factors.

To keep statistical and systematic errors under control, only a small and carefully chosen subset of all the available time separations $0\leq t_{\rm K}+t_{\textnormal{B}_{\textnormal{s}}}<T$ of the three-point function in the $t_{\textnormal{B}_{\textnormal{s}}}$-$t_{\rm K}$-plane can be used in the fit.

Moreover, it is important to have a criterion for good determination of the ground-state form factors by the combined fit. However, it is not easy to find a simple strict criterion, like the one from Sec. 3.1 for $E_{\rm K}$, and we monitor several criteria.

Details of both issues, the choice of the fit ranges and the quality of the fit, are discussed in the following.

4.1 Fit ranges

For the three-point functions we select the data points to be included in the combined fit by the constraints $t_{\rm B3}^{\rm min}\leq t_{\textnormal{B}_{\textnormal{s}}}$ and $t_{\rm K3}^{\rm min}(t_{\textnormal{B}_{\textnormal{s}}})\leq t_{\rm K}\leq t_{\rm K3}^{\rm max}(t_{\textnormal{B}_{\textnormal{s}}})$.

The function $t_{\rm K3}^{\rm max}(t_{\textnormal{B}_{\textnormal{s}}})$ is chosen, for each $\mu$ and smearing level separately555 In [32, 10] we used one and the same value of $t_{\rm K3,\mu}^{\rm max}$ for all $t_{\textnormal{B}_{\textnormal{s}}}$, corresponding roughly to the one that produced a rectangle of largest area in the plane of the allowed points (for both values of $\mu$ separately). , as

[TABLE]

were $t_{\rm K3}^{\rm max,wr}(t_{\textnormal{B}_{\textnormal{s}}})$ is given by the wrapper criterion of (3.30) and $t_{\rm K3}^{\rm max,noise}(t_{\textnormal{B}_{\textnormal{s}}})$ is determined from a relative noise criterion, i.e. as the largest $t_{\rm K}$ which fulfills

[TABLE]

with $c_{\rm noise}=2.5\%$. The corresponding curves are shown in Fig. 7 and illustrate that A5 is an example of a wrapper-limited lattice (due to a small physical time extent). On the other hand, O7 is noise limited and by reducing statistical errors one may hope to improve the determination of the form factors.

Also the function $t_{\rm K3}^{\rm min}(t_{\textnormal{B}_{\textnormal{s}}})$ is selected by an automatic criterion666In Ref. [10] $t_{\rm K3}^{\rm min}$ was a manually tunable parameter of the fit proceduce, chosen in common for all smearings and $\mu$.

analogous to the one described in Sec. 3.1: for each $t_{\textnormal{B}_{\textnormal{s}}}$, $\mu$, and smearing, we do a two-exponential fit to $t_{\rm K}\in[0.4\;{\rm fm},t_{\rm K3}^{\rm max}(t_{\textnormal{B}_{\textnormal{s}}})]$ (with the energies fixed to the ones extracted from the two-point function). We then find the minimum time $t_{\rm K3*}^{\rm min}(t_{\textnormal{B}_{\textnormal{s}}})$ at which the excited-state contribution is smaller than 25% of the statistical error. The final value is set to

[TABLE]

to avoid values getting too small in the regions with large noise, where the two-state fit does not work well. The resulting $t_{\rm K3}^{\rm min}$ is shown by the blue curves in Fig. 7.

The remaining tunable parameters of the fit are $t_{\rm B2}^{\rm min}$ and $t_{\rm B3}^{\rm min}$. In Ref. [10] we used $t_{\rm B2}^{\rm min}\approx 0.45$ fm, which seems sufficient when $N_{{\textnormal{B}_{\textnormal{s}}}}=3$. With better smearings and improved analysis methods, we can use slightly lower values, with two choices: $t_{\rm B2}^{\rm min}\approx 0.32\,{\rm fm}$ (“aggressive”) and $0.38\,{\rm fm}$ (“conservative”).

Concerning $t_{\rm B3}^{\rm min}$ it was observed in Ref. [10] that the fit performs better if one takes a positive $t_{\rm B\Delta}^{\rm min}\equiv t_{\rm B3}^{\rm min}-t_{\rm B2}^{\rm min}$. This can be understood by noting that the amplitudes in the two-point function are proportional to $|\beta_{i}^{(n)}|^{2}$, while in the three-point functions they are only proportional to $\beta_{i}^{(n)}$. Thus, the suppression of excited states is expected to be stronger in the two-point functions if $\beta_{i}^{(0)}>\beta_{i}^{(1)}$. We find that this is in fact the case, except for the lowest smearing.777For the ensembles we analysed, typical values are: $\beta_{3}^{(1)}/\beta_{3}^{(0)}\approx 0.25-0.4$, $\beta_{2}^{(1)}/\beta_{2}^{(0)}\approx 0.67$, $\beta_{1}^{(1)}/\beta_{1}^{(0)}\approx 1$.

With increasing time separation $t_{\textnormal{B}_{\textnormal{s}}}$ the suppression of the first excited state is proportional to $e^{-(E_{\textnormal{B}_{\textnormal{s}}}^{(1)}-E_{\textnormal{B}_{\textnormal{s}}}^{(0)})t}$, which approximately amounts to 0.67 at $t=0.2\;{\rm fm}$, 0.4 at $t=0.4\;{\rm fm}$, and 0.25 at $t=0.5\;{\rm fm}$. Although it is not immediately obvious which value one should take, this justifies a non-zero value of $t_{\rm B\Delta}^{\rm min}$. We keep $t_{\rm B\Delta}^{\rm min}\approx 0.2\,{\rm fm}$, which gives a good compromise in suppressing the higher excited-state contributions and keeping the number of data points large enough for good fit stability.

4.2 Fit quality and results

The combined fit is uncorrelated, i.e. we use a diagonal covariance matrix. For such a big matrix the full covariance matrix is very badly determined and conditioned (in fact, if the number of data points is bigger than the number of measurements, it is not possible to invert [38, 39]).

As a way to monitor the quality of the uncorrelated fit, we use $\chi^{2}_{\rm exp}$, the expectation value of $\chi^{2}$ given normally distributed data with the measured covariance matrix [40]. It can be estimated efficiently including autocorrelations [41], but the uncertainty on its determination can be large and it is not easy to exactly quantify what is a good fit. In general we consider fits with $\chi^{2}\lesssim 2\chi^{2}_{\rm exp}$ as acceptable.

Another way of monitoring the fit quality is to look at the contributions to the $\chi^{2}$ coming from different correlation functions. In particular, an unusually large contribution from $\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}$ should be monitored as it can indicate that $t_{\rm B2}^{\rm min}$ was chosen too small (even though the overall $\chi^{2}$ looks acceptable because of the many points in the three-point functions).

Apart from monitoring the $\chi^{2}$, we analyze the stability of the extracted fit values and errors with respect to changes of $t_{\rm B2}^{\rm min}$ and $t_{\rm B3}^{\rm min}$.

The stability plots are always organized in the following way: the data is divided in groups with different $t_{\rm B2}^{\rm min}$, which is also plotted on the $x$-axis in the middle of each group. Inside every group $t_{\rm B\Delta}^{\rm min}$ varies from 0 to approx. 0.5 fm. Additionally, we highlight in different colors the values for selected $t_{\rm B2/3}^{\rm min}$, which are described in the following.

In choosing $t_{\rm B2}^{\rm min}$ and $t_{\rm B3}^{\rm min}$ one must find a good window between too small values, where the fit is plagued by contamination from the higher excited states (this can be seen by high $\chi^{2}$ as well as by lack of stability of the extracted fit parameters with respect to small changes in $t_{\rm B2/3}^{\rm min}$) and too high values, where it can no longer resolve the three ${\textnormal{B}_{\textnormal{s}}}$ states. The latter is demonstrated in Fig. 8. At large $t_{\rm B2/3}^{\rm min}$ we see that for both ensembles at least one of the energies is no longer resolved. Therefore, we refrain from showing the largest values of $t_{\rm B2}^{\rm min}$ in Figs. 9 – 11.

In general, the ground-state energies (cf. Fig. 9) and the ground-state form factors (cf. Fig. 10) are reasonably stable with respect to changes of $t_{\rm B2/3}^{\rm min}$ (within quickly growing errors). We also show representative examples for the excited-state form factors in Fig. 11. Their precision, especially of $\varphi^{(2)}_{\mu}$, is much worse than of their ground-state counterparts.

The values for the ground-state form factors extracted with both choices of $t_{\rm B2}^{\rm min}$ are collected in Table 3, together with the results from ratio methods described in the next section.

5 Matrix elements at static order from ratios

5.1 Ordinary ratios

An alternative to extracting the ground-state matrix elements from fitting the correlation functions directly is to construct an appropriate ratio of the correlation functions such that the dependence on all or most of the other parameters cancels in the limit of large time.

One example of such a ratio was actually used to define the form factor, cf. Eq. (LABEL:eqn:me1). In the following it will be beneficial to generalize this definition to the case where ${t_{\rm K}}\neq{t_{{\textnormal{B}_{\textnormal{s}}}}}$ in the three-point function ($\tau\equiv{t_{\rm K}}+{t_{{\textnormal{B}_{\textnormal{s}}}}}$):

[TABLE]

In the limit ${t_{\rm K}},{t_{{\textnormal{B}_{\textnormal{s}}}}}\to\infty$, this ratio converges to the desired bare form factor $\varphi^{(0,0)}_{\mu}$ if additional parameters $\tilde{E}_{\rm K}(\tau),\tilde{E}^{(0)}_{\textnormal{B}_{\textnormal{s}}}(\tau)$ satisfy $\tilde{E}(\tau)=E+\mathcal{O}\left(\exp(-\Delta E\tau)\right)$. Possible choices will be discussed below.

In Sec. 5.2 we will see the advantage of this definition, while its obvious disadvantage is that in the denominator one needs the two-point correlation functions at time $\tau={t_{\rm K}}+{t_{{\textnormal{B}_{\textnormal{s}}}}}$. This means that e.g. for a simple estimate of a plateau in the three point function at ${t_{\rm K}}\approx{t_{{\textnormal{B}_{\textnormal{s}}}}}\approx 1\,{\rm fm}$ one needs the two-point correlation functions at a large time separation $\tau\approx 2\,{\rm fm}$, which is particularly problematic in the case of ${\cal C}^{\textnormal{B}_{\textnormal{s}}}$ where the signal-to-noise problem is much more severe than for $\mathcal{C}^{\rm K}$.

For lattices with short time extent, like A5, $\tau$ may come close to the middle of the lattice and we need to take into account the influence of the second term in Eq. (3.24). We do that by multiplying the ratio by an additional factor \big{(}1+\operatorname{e}^{E_{\rm K}(T-4t)}\big{)}^{1/2} which cancels the unwanted contribution up to excited states which are negligible in this region.

In addition, one can consider another ratio

[TABLE]

which has a more favourable signal-to-noise behaviour of the denominator but requires energy estimates even for the choice ${t_{\rm K}}={t_{{\textnormal{B}_{\textnormal{s}}}}}$.

A similar ratio is

[TABLE]

where the additional normalization factors $\mathcal{N}^{\rm K}$, $\mathcal{N}_{r}^{\textnormal{B}_{\textnormal{s}}}$ can be expressed in terms of the correlation function parameters $\kappa^{(0)}$ and $\beta^{(0)}_{r}$ respectively. In practice, our data for $\mathcal{R}^{\,\textrm{III}}$ shows very similar or slightly inferior behaviour compared to $\mathcal{R}^{\,\textrm{II}}$, therefore we do not present the numerical results for this ratio.

Let us now discuss the choice of the ground-state energy estimates. One can set them to the effective masses of the corresponding two-point functions, avoiding any fit procedure. This, however, results in large statistical fluctuations in the large-time region. For our data, it is beneficial to instead use the time-independent estimates $E_{\rm K}^{(0)},E_{\rm{\textnormal{B}_{\textnormal{s}}}}^{(0)}$ extracted in Sections 3.1 and 3.2 respectively. The plateaux in both cases were chosen in a conservative manner, and as a crosscheck we calculated the ratios for different choices of the plateau ranges of both energies. From that exercise we find that the systematic error associated with using the fitted ground-state energies is negligible with respect to our statistical uncertainties.

The resulting ratios $\mathcal{R}^{\,\textrm{I}}$ and $\mathcal{R}^{\,\textrm{II}}$ are presented in Figs. 12 and 13, respectively. In addition to the data we also plot the fit results from Sec. 4, including both the fitted values of $\varphi^{(0,0)}_{\mu}$ and the respective ratios of correlation functions reconstructed by inserting the fit parameters ($\varphi^{(0,n)}_{\mu}$, $\kappa^{(0)}$, $\beta_{r}^{(n)}$, $E_{\rm K}^{(0)}$ and $E_{\textnormal{B}_{\textnormal{s}}}^{(n)}$) into Eqs. (3.20)–(3.22). There is a good agreement of the actual ratios with their fitted counterparts at large values of $t$, where excited states neglected in the fits are irrelevant.

On ensemble A5, due to the short time extent of the lattice, the wrapper criterion Eq. (3.30) becomes relevant. For $\mu=1$ it restricts the available times to approx. 1.1 fm, cf. Fig. 7. At larger times we observe that the ratio $\mathcal{R}^{\,\textrm{II}}_{1}$ starts to grow rapidly, while for $\mathcal{R}^{\,\textrm{I}}_{1}$ this is most likely masked by the large uncertainties. Furthermore, for $\mu=0$ the data is close to violating the wrapper criterion (see Fig. 7) at the maximal time separations used.

On both lattices, A5 and O7, we see that at large times, $t\gtrsim 1$ fm, $\mathcal{R}^{\,\textrm{II}}$ is clearly superior to $\mathcal{R}^{\,\textrm{I}}$ in terms of the signal to noise and has a comparable precision to the results of the combined fit. To use $\mathcal{R}^{\,\textrm{II}}$ for the extraction of $\varphi^{(0,0)}_{\mu}$, we need to select a suitable plateau: looking at the bands from the fit which give us an estimate of the excited-state contamination from the ${\textnormal{B}_{\textnormal{s}}}$ sector, we start the plateau at $t^{\rm min}=1.1$ fm and fit until the loss of precision of the signal below 5% or until we hit the wrapper criterion. On A5 for $\mu=1$ there is no valid plateau, because the wrapper is hit before $t^{\rm min}$. In this case, we quote the first data point above 1.1 fm as the final result and add to it a systematic error associated with the wrapper contribution estimated by the fit.

5.2 Summed ratios

A way to get improved convergence to the ground state is to sum888We write this in the form used in our analysis, but note that often it is advantageous and natural to sum the position of the operator that causes the transition over all spacetime instead [42, 43]. the ratio [42, 44, 43] $\mathcal{R}^{\,\textrm{I}}$ and determine the matrix element from

[TABLE]

The asymptotic excited-state contaminations are then $\mathcal{O}(\tau\Delta{\rm e}^{-\tau\Delta})$, where $\Delta=\min(E_{\rm K}^{(1)}\!\!-E_{\rm K}^{(0)},E_{\rm{\textnormal{B}_{\textnormal{s}}}}^{(1)}\!\!-E_{\rm{\textnormal{B}_{\textnormal{s}}}}^{(0)})$, as opposed to $\mathcal{O}({\rm e}^{-\tau\Delta/2})$ in ordinary ratios [44, 43].999One could in principle attempt to sum the other ratios as well. However, the improvement in convergence is only proven to work (at least asymptotically) for $\mathcal{R}^{\,\textrm{I}}$. Empirically we see that it does not work for $\mathcal{R}^{\,\textrm{II}}$ and $\mathcal{R}^{\,\textrm{III}}$. The accelerated convergence of the summed ratios can be a decisive advantage in the case of ensembles with a limited extent in time such as A5: contaminations by wrappers are less relevant. Example results are shown in Fig. 14. We observe that the convergence is indeed improved.

In fact, especially for $\mu=0$ the plateaux seem to start very early. This early onset of the plateaux is discussed in Appendix C where we show that it is caused by an accidental cancellation of the excited-state contributions in the ${\textnormal{B}_{\textnormal{s}}}$ and ${\rm K}$ sector and therefore should be treated with caution.

Note, however, that the ratios reconstructed from the fits of Sec. 4 only include the influence of the excited states in the ${\textnormal{B}_{\textnormal{s}}}$ sector. As in the previous subsection, these fit bands can be used to estimate the start of the plateaux. We choose $\tau^{\rm min}=1$ fm and $\tau^{\rm max}$ by the 5% relative noise criterion.101010One could also choose 10%, the plateaux would be a bit longer but that does not help in reducing the statistical error. In this range we extract the form factor from the slope of the linear fit directly to $\mathcal{S}^{\textrm{I}}(\tau)$. Note that on both lattices we are not affected by the wrapper criterion as opposed to ordinary ratios.

The results from the different methods are gathered in Table 3 and are in good agreement. One slight difference is seen in $\mathcal{R}^{\,\textrm{II}}_{0}$ on ensemble A5. However, in this case if we start the plateau one data point earlier (which, judging from Fig. 13, seems legitimate) the result goes down to 1.117(15) which removes most of the difference. In general the results from the fits and the summed ratio have similar precision, in which case we prefer to quote the result from the summed ratio as the final results, as the latter method is simpler.

5.3 GEVP method

We have also implemented the GEVP method of Ref. [43], both the ratio one, eq. (2.16) of [43], and the summed one, eq. (3.7) of [43]. In both cases we specialized to using the GEVP only in the ${\textnormal{B}_{\textnormal{s}}}$ channel. The results are very similar to the ones of the previous two sections (with the largest wave function), but with errors which are a little bit larger. The fact that we do not see a significant improvement might be due to our interpolating fields, which likely do not distinguish between single hadron and (excited) multi-hadron states. The GEVP is then not able to significantly reduce the multi-hadron contributions. We comment further on this issue in the conclusions.

6 Matrix elements at $1/m$ order

The HQET expansion of B-meson observables becomes a precision tool only when $1/m$ terms are included. We now discuss the determination of these crucial terms for the matrix elements. As in the static approximation we have the option to perform fits to the (two-point and three-point) correlation functions, or to consider ratios or summed ratios. Fits have been discussed and applied in Ref. [18] for the somewhat simpler case of the B-to-vacuum matrix element $f_{\textnormal{B}_{\textnormal{s}}}$. As there, in our present case, all parameters $\beta_{s}^{(n)},\,E^{(n)}$ get triplicated with static and $1/m$ pieces corresponding to kin and spin insertions. Furthermore pure exponentials turn into exponentials plus terms of the form $E^{\mathrm{kin}}_{n}t\exp(-E^{{\rm stat}}_{n}t)$. This proliferation of terms in addition to the $1/m$ corrections to the matrix elements themselves makes an analysis in terms of fits cumbersome and difficult. We therefore concentrate on the analysis of ratios, which we already have seen to be just as good as fits in the static approximation. We will also see that there are some simplifications in the $1/m$ expansion of ratios which make them rather accessible.

As a preparation for expanding the ratios, we start with the underlying two- and three-point functions. First note that we choose the arbitrary interpolating fields ${\cal O}_{\mathrm{bs},r}$, eq. (2.14), not to contain a $1/m$ piece. Therefore, the $1/m$ expansion of the two-point function reads

[TABLE]

Only the (dimensionful) HQET parameters of the action, $m_{\textnormal{bare}}$ as well as $\omega_{\mathrm{k}}\sim 1/m$, appear. Here and everywhere below all O$(1/m^{2})$ terms are dropped without notice. Note also that $m_{\textnormal{bare}}$ drops out in the expressions for matrix elements.

The energy $E^{\textnormal{B}_{\textnormal{s}}}=\lim_{t\to\infty}-\partial_{t}\log(\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}(t))$ is $1/m$-expanded as $E^{\textnormal{B}_{\textnormal{s}}}=m_{\textnormal{bare}}+E^{{\rm stat}}+\omega_{\mathrm{kin}}E^{\mathrm{kin}}+\omega_{\mathrm{spin}}E^{\mathrm{spin}}$ with $E^{x}$ appearing in the large time behaviour of [45]

[TABLE]

All $E^{x}$ refer to the ground state; the first excited state contribution leads to the term with $\Delta E^{{\rm stat}}=E^{{\rm stat}}_{2}-E^{{\rm stat}}_{1}$. Because one always first expands in $1/m$ and then takes a limit of large time, terms such as $\exp(-\Delta E^{\mathrm{kin}}t)$ do not appear.

In the numerical applications we will take $\partial_{t}$ to be the forward derivative, $a\partial_{t}f(t)=f(t+a)-f(t)$. On integrating Eq. (6.40) , we get:

[TABLE]

The integration constants ${\rm A}^{{\textnormal{B}_{\textnormal{s}}},\textnormal{k}}_{r}$ do depend on the $\mathrm{kin}$ and $\mathrm{spin}$ insertions as well as the smearing level used.

In complete analogy we have for the three-point functions ($t=t_{\textnormal{B}_{\textnormal{s}}}=t_{\rm K}$):

[TABLE]

Now we turn to the the ratios and insert the $1/m$ expansion of the quantities which enter their definition in (5.34) and (5.35). To understand the structure, consider the expansion

[TABLE]

with

[TABLE]

The sums over $k$ run as indicated and for $j$ the range is seen in Table 5. It follows that the desired $1/m$ corrections to the ground state matrix elements are given by the large time limits,

[TABLE]

These equalities are what we alluded to before as simplification in the $1/m$ expansion. In fact, since the ground state matrix elements in the last expression are just given by the $A$ terms, it should not come as a surprise that the same final formulae for $\rho^{j}_{\mu}$ and $\rho^{k}_{\mu}$ are obtained if one starts from ratios $\mathcal{R}^{\,\textrm{II}}$ or $\mathcal{R}^{\,\textrm{III}}$.

The strategy to obtain the $1/m$ terms of the matrix elements is then to extract ${\rm A}^{{\textnormal{B}_{\textnormal{s}}},k}_{r}$ and ${\rm A}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K},k}_{\mu,r}$ from fits to Eq. (6.41) and Eq. (6.42) and $\rho^{j}_{\mu}$ from Eq. (6.46).

6.1 Numerical results for the kin and spin insertions at $1/m$ order

As a first analysis we performed combined fits with parameters ${\rm A}^{{\textnormal{B}_{\textnormal{s}}},k}_{r}$, ${\rm A}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K},k}_{\mu,r}$ and $E^{k}$ to Eq. (6.41) and Eq. (6.42) in the time region where a linear behavior in $t$ is observed. This was done for fixed insertion $k\in\{\mathrm{kin},\mathrm{spin}\}$ and fixed smearing level, typically the largest smearing. The resulting errors on the matrix elements were rather large because the data did not constrain the fitted energies $E^{k}$ so well. An improvement could be achieved by determining $E^{k}$ from the GEVP at order $1/m$, exactly as described in Ref. [46] and then using that as a constraint in Eq. (6.41) and Eq. (6.42). The GEVP takes into account information from all smearing levels of the two point functions. In a little more detail, we expand the ground state GEVP eigenvalues in $1/m$ (see Ref. [46] for explicit formulae) as $\lambda^{(0)}(t,t_{0})=\lambda^{{\rm stat}}(t,t_{0})+\omega_{\mathrm{kin}}\lambda^{\mathrm{kin}}(t,t_{0})+\omega_{\mathrm{spin}}\lambda^{\mathrm{spin}}(t,t_{0})$ and then form

[TABLE]

We consider just $t_{0}\geq t/2$ as the asymptotic convergence is proven to be much better [46] under that condition. For the present case, it turns out that there is rather little dependence on $t_{0}$ in practice but errors of course grow as it is increased. We then form a weighted average of the first (up to) three values with $t_{0}$ at and above $t/2$. These averages are shown as $E^{k}_{\mathrm{eff}}(t)$ in Fig. 15 and 16 together with the effective masses of the best smearing level. The GEVP estimates behave significantly better than just the best smearing which we show for comparison. Our final numbers come from plateaux fits of the GEVP estimates starting at $t=0.5{\rm fm}$, and are shown as bands in Figs. 15 and 16.

Effective ${\rm A}^{{\textnormal{B}_{\textnormal{s}}},k}_{r}(t)$ are determined by inserting the energies into ${\rm A}^{{\textnormal{B}_{\textnormal{s}}},\textnormal{k}}_{r}(t)={{\cal C}}^{{\textnormal{B}_{\textnormal{s}}},\textnormal{k}}_{rr}(t)/{{\cal C}}^{{\textnormal{B}_{\textnormal{s}}},{\rm stat}}_{rr}(t)+E^{\textnormal{k}}t$ and similarly for ${\rm A}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K},k}_{\mu,r}$. These effective A-estimates are shown in Figs. 17 and 18. Precise, early and long plateaux are present for the two-point function ${\rm A}^{{\textnormal{B}_{\textnormal{s}}},k}_{r}$. In contrast, for the three-point function ${\rm A}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K},k}_{\mu,r}$ an agreement with a plateau is only seen starting at $t=0.8\,{\rm fm}$ and in fact we would like the plateaux to be more convincing in one or two cases. Nevertheless, taking weighted averages starting at $t=0.8\,{\rm fm}$ is reasonable and we collect their results combined to $\rho^{k}_{\mu}$ in Table 4. The bands in Figs. 17 and 18 show the chosen fit values.

6.2 Current insertions at $1/m$ order

The $1/m$ vector current contributions are

[TABLE]

with the operators $\Gamma_{\mu,j}$ detailed in Table 5. Note that with a momentum purely along the x-axis (A5 ensemble), only $i=1$ contributes and $\rho_{1}^{3}=\rho_{1}^{4}$ is exact.

The results obtained for Eq. (6.44) are presented in Fig. 19. We see plateaux starting at roughly 0.8-0.9 fm and average from 0.85 fm on. The precision is better than that of the kin and spin terms. A full quantitative error budget for the $1/m$ contributions to the form factors has to wait until the corresponding non-perturbative matching coefficients $\omega$ are available.

7 Discussion

Flavor changing transitions are important channels for learning about possible limitations of the standard model. Exclusive semi-leptonic decays of B or ${\textnormal{B}_{\textnormal{s}}}$ mesons are particularly clean theoretically. However, lattice computations of the relevant form factors are needed and they are non-trivial in practice. A major reason is the infamous signal-to-noise problem. At large Euclidean time separations the noise in Monte Carlo evaluation of the correlation functions is too large to determine the form factors (matrix elements) with interesting precision.

While the issue is not new, we have exposed the problem in a few graphs more clearly than often done, see Figs 4, 14. For the three-point functions and derived (summed) ratios, this was possible because we have evaluated correlation functions at all time separations. We are only considering the pseudoscalar sector, where the signal-to-noise problem is very mild for a relativistic formulation and at zero momentum. However, we need finite momentum (see Fig. 2 for the momentum dependence) and for reasons explained in the beginning of the paper we use HQET for the b-quark (see Fig. 4 for the difference of Kaon and static ${\textnormal{B}_{\textnormal{s}}}$ meson). In HQET, the signal-to-noise problem becomes worse as one decreases the lattice spacing. Nevertheless, we are determining the matrix elements at total time separations $\tau$ of the three-point functions of around 2 fm and for $a\approx 0.05\,{\rm fm}$. Only for the summed ratio we use about half that time-separation (Fig. 14), in agreement with the predicted better suppression of excited states after summation.

Despite our use of HQET, these separations are larger or similar to the ones typically used. E.g. most recently Ref. [9] used a fixed $\tau\approx 2.2\,{\rm fm}$.

We have presented good evidence that the chosen plateaux or fit-windows are reasonably safe, but nevertheless it would be better to have larger times accessible. Maybe multilevel strategies  [47, 48] will help to reach those separations in the future. At present we derive our confidence from the good agreement of

• •

fits with {t_{{\textnormal{B}_{\textnormal{s}}}}}\raisebox{-2.58334pt}{\stackrel{{\scriptstyle\textstyle{>}}}{{\sim}}}0.3-0.6 and 3 states in the ${\textnormal{B}_{\textnormal{s}}}$ sector (and also with 2 ${\textnormal{B}_{\textnormal{s}}}$ states and {t_{{\textnormal{B}_{\textnormal{s}}}}}\raisebox{-2.58334pt}{\stackrel{{\scriptstyle\textstyle{>}}}{{\sim}}}0.5-0.7, cf. Appendix B)

• •

ratios Eq. (5.35) with \tau={t_{{\textnormal{B}_{\textnormal{s}}}}}+{t_{\rm K}}\,\raisebox{-2.58334pt}{\stackrel{{\scriptstyle\textstyle{>}}}{{\sim}}}\,2\,{\rm fm},

• •

and summed ratios Eq. (5.37) at total separation \tau\,\raisebox{-2.58334pt}{\stackrel{{\scriptstyle\textstyle{>}}}{{\sim}}}\,1\,{\rm fm} .

Note that also the precision of the different estimates is quite comparable. It is therefore preferable to use the technically easier ratio methods.

At the lowest (static) order in $1/m$ the most relevant form factors for $\mu=1$ have an accuracy around 2%, which is good for precision physics. The first order corrections in $1/m$ are actually more precise than that. The relative errors induced into the form factors are given by the absolute ones of the quantities $\rho$ in Tables 4,5 multiplied with the appropriate $\omega$-coefficients. From the errors of $\rho^{\mathrm{kin}}$ and $\rho^{\mathrm{spin}}$, using the non-perturbative $\omega_{\mathrm{kin}}\,,\omega_{\mathrm{spin}}$ [13], we get error contributions of 1% and 0.5% to the form factors, respectively. The uncertainties of $\rho^{j}$, translate into below 0.5% errors assuming coefficients $\omega_{j}$ which do not exceed the tree-level values by more than a factor of two. Of course, the uncertainties in the coefficients $\omega$ are to be added separately.

A positive result of our detailed analysis is thus that $1/m$ corrections can be determined precisely, when the coefficients (HQET parameters) are known with reasonable accuracy. Once they are available it will be possible to provide further crosschecks on the existing analysis, see e.g. Fig. 13 of [9].

We would finally like to point out a possible danger in current and future semi-leptonic form factor computations on the lattice. It consists in the contribution of multi-hadron states, such as $|{\textnormal{B}_{\textnormal{s}}}^{(\ast)},\pi\rangle$ or $|{\textnormal{B}_{\textnormal{s}}}^{(\ast)},{\rm K}\rangle$, to our Euclidean correlation functions. Formally, these contributions are just particular ones in the sums over excited states in Eqs. (3.20-3.22) and are thus covered by our analysis. However, there are two properties which make such states special. First, when the spatial volume becomes large and the light-quark masses small, there are several low-lying states with small gaps $E^{(n)}-E^{(0)}$. All used methods may have difficulties in separating those. This has also been pointed out recently in [49]. Second, normalized overlaps $\beta_{r}^{(n)}/\beta_{r}^{(0)}$ may be very similar for different smearing levels $r$ when $n$ corresponds to a multi-hadron state.111111See the discussion of the field in the chiral effective theory representing the action density at finite flow time, $E(x,t)$ in Ref. [50]. As mentioned in Sect. 5.3 the GEVP-method is of little help in such a situation. For B-mesons these contaminations at finite ${t_{{\textnormal{B}_{\textnormal{s}}}}}$ have not been investigated at all, while there is considerable discussion in the analogous determination of nucleon matrix elements (see Refs. [51, 52] and references therein). It appears likely that form factors ${\textnormal{B}_{\textnormal{s}}}\to K$, $\rm B\to\pi$ are in somewhat better shape than nucleon matrix elements, since larger Euclidean time separations are reached. Nevertheless systematic studies, especially theoretical ones such as the ones carried out in chiral perturbation theory [53, 54] for nucleon matrix elements are urgently needed to make quantitative statements.

Appendix A Computation of the correlation functions

Integration over the fermion fields yields the correlation functions as

[TABLE]

averaged over the gauge fields, $U$, (with effective action including the log of the quark determinant). The dependence of the functions

[TABLE]

on the times $(x_{0})_{s}$ ($s=\mathrm{f},\mathrm{v},\mathrm{i}$) as well as the gauge fields is suppressed. In the above expressions $S_{q}$ are the quark propagators, $W_{\rm K}=G^{I^{\rm K}},\;G=1-\kappa_{\mathrm{G}}\,a^{2}\Delta$ and $W_{r}=G^{I^{\textnormal{B}_{\textnormal{s}}}_{r}}$ are the smearing operators, $P_{s}$ is the projection on time slice $(x_{0})_{s}$ and $F_{\pm}=\exp(\pm i{\bf x}{\bf p})$ is the multiplication with the Fourier phases on the corresponding time slice. Evaluation of the full traces which include the space coordinates becomes possible by representing the projector

[TABLE]

in terms of a random U(1) field $\eta$ with support on on time-slice $(x_{0})_{\mathrm{f}}$. This yields

[TABLE]

and similar for the two-point functions. In this form we see that for each vector $\eta$ two solutions of the Dirac equation are needed in order to compute the fields $\phi_{\mathrm{u}},\phi_{\mathrm{s}}$. The static propagator $S_{\mathrm{h}}$ is inserted by explicit forward propagation. Translation invariance in time is used by averaging over all source time-slices $(x_{0})_{\mathrm{f}}$, each one with a random U(1) field $\eta$. Since the averages Eq. (A.50) and Eq. (A.54) are independent of each other, any number of $\eta$ fields per gauge field is correct; we use a single one per gauge field and time $(x_{0})_{\mathrm{f}}$.

In practice, we compute propagators of periodic quark fields in a gauge field $e^{i\theta_{\mu}a/L}U(x,\mu)$, which includes the constant $U(1)$ background field $\theta^{q}_{\mu}a$ and only use the integer part of the momentum in the Fourier factors $F_{\pm}$. This is equivalent to Eq. (2.18), apart from a phase in the Gaussian smearing, which we set to zero. Choosing a different phase in the smearing along the lines of [55] might be a further optimization.

A.1 Improvement and $1/m$ terms

The $1/m$ terms are simple generalizations of the above. First, $[{\cal C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu}]_{k}\,,k\in$ kin, spin and $[\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}_{rr^{\prime}}]_{k}$, are given by replacing $S_{\mathrm{h}}\to S_{k}$, where the latter are

[TABLE]

and the chromo-magnetic field strength tensor $\hat{F}_{ij}$ is discretized in terms of the clover leaf, see e.g. [22].

Second, the NLO three-point functions $[{\cal C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu}]_{j}$ are given by the substitution $\gamma_{\mu}\to\Gamma_{\mu}^{j}$ in Eq. (A.51) and Eq. (A.56), with $\Gamma_{\mu}^{j}$ listed in Table 5. The O($a$) improvement corrections to the static three-point functions are just linear combinations thereof.

Appendix B Alternative static fit with $N_{\rm K}=1,N_{\textnormal{B}_{\textnormal{s}}}=2$

In this section we discuss a simpler version of the fit which includes only one excited state in the ${\textnormal{B}_{\textnormal{s}}}$ sector. We use only the two highest ${\textnormal{B}_{\textnormal{s}}}$ smearings, therefore we have a total of 12 fit parameters (as opposed to 20 in the fit with $N_{\textnormal{B}_{\textnormal{s}}}=3$).

Similar to Sec. 4 we show the ${\textnormal{B}_{\textnormal{s}}}$ energy states. Here, the growth of the errors at larger $t_{\rm B2/3}^{\rm min}$ is much milder than for the $N_{\textnormal{B}_{\textnormal{s}}}=3$ fit, but on the other hand the stability of the results wrt. changes of $t_{\rm B2/3}^{\rm min}$ is in some cases rather unsatisfying. It can be even more clearly seen for $E_{\rm{\textnormal{B}_{\textnormal{s}}}}^{(0)}$ (cf. Fig. 21) which, in particular for O7, shows a characteristic drift upwards with growing $t_{\rm B3}^{\rm min}$.

On the other hand, considering these drifts in energies (and amplitudes as well), the ground-state form factors show reasonable stability and agreeement with the $N_{\textnormal{B}_{\textnormal{s}}}=3$ fit results, cf. Fig. 22. The case of $\varphi_{0}^{(0)}$ for O7 is by far the worst, but most look more like the stable and consistent $\varphi_{1}^{(0)}$.

When choosing a final form factor value from the fit, one clearly has to keep $t_{\rm B2/3}^{\rm min}$ slightly higher than for the $N_{\textnormal{B}_{\textnormal{s}}}=3$ fit. We use $t_{\rm B2}^{\rm min}\approx 0.52$ fm and the same $t_{\rm B\Delta}^{\rm min}$ as before. The values summarized in Table 6 are consistent within errors with other methods and give similar precision.

Appendix C Note on the convergence of the static summed ratio

In Section 5.2, devoted to static summed ratios, we observe that the plateaux for $\mu=0$ in Fig. 14 seem to start very early. Taking the results for O7 at face value, one could start the plateau fit as early as 0.4 fm. To see whether this is a genuine ground state dominance we investigate truncated sums of the form:

[TABLE]

By putting $t_{\textnormal{B}_{\textnormal{s}}}^{\rm min}$ or $t_{\rm K}^{\rm min}$ larger than 0 we can suppress the excited states in the ${\textnormal{B}_{\textnormal{s}}}$ and ${\rm K}$ sector respectively.121212Note however that this increases the excited state corrections from $\sim\tau\Delta\exp(-\tau\Delta)$ to $\sim\tau^{\prime}\Delta\exp(-\tau^{\prime}\Delta)$ with $\tau^{\prime}=\tau-t_{\textnormal{B}_{\textnormal{s}}}^{\rm min}-t_{\rm K}^{\rm min}$. The results of this procedure are shown in Fig. 23. We see that the excited states in the ${\rm K}$ sector push the result upwards, while the excited states in the ${\textnormal{B}_{\textnormal{s}}}$ sector push the results downwards. Their cancellation results in a flat “fake” plateau which can start at very early times. A similar but much less pronounced effect is at work for $\mu=1$. One should therefore not rely solely on the flatness of the plots but also devise an independent criterion for the beginning of the plateaux, as was done in Section 5.2.

Acknowledgements. We would like to acknowledge useful discussions with Oliver Bär, Stefan Schaefer, Michele Della Morte, Brian Colquhoun.

We gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS share of the supercomputer JUQUEEN at Jülich Supercomputing Centre (JSC). GCS is the alliance of the three national supercomputing centres HLRS (Universität Stuttgart), JSC (Forschungszentrum Jülich), and LRZ (Bayerische Akademie der Wissenschaften), funded by the German Federal Ministry of Education and Research (BMBF) and the German State Ministries for Research of Baden-Württemberg (MWK), Bayern (StMWFK) and Nordrhein-Westfalen (MIWF). We acknowledge PRACE for awarding us access to resource JUQUEEN in Germany at Jülich and to resource SuperMUC in Germany at München. We also thank the LRZ for a CPU time grant on SuperMUC, project pr85ju, and DESY for access to the PAX cluster in Zeuthen.