# Extraction of bare Form Factors for
$\mathrm{B}_{\mathrm{s}}\to\mathrm{K}\ell\nu$ Decays in non-perturbative HQET

Felix Bahr  Debasish Banerjee  Fabio Bernardoni  Mateusz Koren  Hubert Simma
Rainer Sommer

###### Abstract

We discuss the extraction of the ground state $\langle{\rm K}({\bf
p})|V_{\mu}(0)|{\textnormal{B}_{\textnormal{s}}}({\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.

###### Contents

  1. 1 Introduction
  2. 2 Setup for the non-perturbative evaluation
    1. 2.1 Form factors from Euclidean correlation functions
    2. 2.2 Lattice Setup
  3. 3 Static correlation functions
    1. 3.1 Kaon correlation functions
    2. 3.2 ${\textnormal{B}_{\textnormal{s}}}$ correlation functions
      1. 3.2.1 Energies from the GEVP
      2. 3.2.2 Amplitudes
    3. 3.3 Three-point static correlation functions
  4. 4 Matrix elements at static order from a combined fit
    1. 4.1 Fit ranges
    2. 4.2 Fit quality and results
  5. 5 Matrix elements at static order from ratios
    1. 5.1 Ordinary ratios
    2. 5.2 Summed ratios
    3. 5.3 GEVP method
  6. 6 Matrix elements at $1/m$ order
    1. 6.1 Numerical results for the kin and spin insertions at $1/m$ order
    2. 6.2 Current insertions at $1/m$ order
  7. 7 Discussion
  8. A Computation of the correlation functions
    1. A.1 Improvement and $1/m$ terms
  9. B Alternative static fit with $N_{\rm K}=1,N_{\textnormal{B}_{\textnormal{s}}}=2$
  10. 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

Figure 1: The correlation function used to extract the
${\textnormal{B}_{\textnormal{s}}}\to{\rm K}$ decay matrix element.

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

$\displaystyle f_{+}(q^{2})$ $\displaystyle=$
$\displaystyle\displaystyle\frac{1}{\sqrt{2m_{\textnormal{B}_{\textnormal{s}}}}}\,\left[h_{\parallel}+(m_{\textnormal{B}_{\textnormal{s}}}-E_{\rm
K})\,h_{\perp}\right]\,,$ (2.1) $\displaystyle f_{0}(q^{2})$ $\displaystyle=$
$\displaystyle\displaystyle\frac{\sqrt{2m_{\textnormal{B}_{\textnormal{s}}}}}{m_{\textnormal{B}_{\textnormal{s}}}^{2}-m_{\rm
K}^{2}}\,\Bigl{[}(m_{\textnormal{B}_{\textnormal{s}}}-E_{\rm
K})\,h_{\parallel}+\left(E_{\rm K}^{2}-m_{\rm
K}^{2}\right)\,h_{\perp}\Bigr{]}\,.$ (2.2)

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

$\displaystyle(2m_{{\textnormal{B}_{\textnormal{s}}}})^{-1/2}\langle{\rm
K}({\bf p}_{\rm
K})|\hat{V}^{0}_{\mathrm{M}}(0)|{\textnormal{B}_{\textnormal{s}}}({\bf
0})\rangle$ $\displaystyle=$ $\displaystyle h_{\parallel}(E_{\rm K})\,,$ (2.3)
$\displaystyle(2m_{{\textnormal{B}_{\textnormal{s}}}})^{-1/2}\langle{\rm
K}({\bf p}_{\rm
K})|\hat{V}^{k}_{\mathrm{M}}(0)|{\textnormal{B}_{\textnormal{s}}}({\bf
0})\rangle$ $\displaystyle=$ $\displaystyle p^{k}_{\rm K}h_{\perp}(E_{\rm
K})\,,$ (2.4)

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,

$\displaystyle\begin{split}{\cal C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K}}_{\mu}({t_{\rm K}},{t_{{\textnormal{B}_{\textnormal{s}}}}};{\bf p}_{\rm
K})=&\frac{a^{9}}{L^{3}}\sum_{{\bf x}_{\mathrm{f}},{\bf x}_{\mathrm{v}},{\bf
x}_{\mathrm{i}}}e^{-i{\bf p}_{\rm K}\,({\bf x}_{\mathrm{f}}-{\bf
x}_{\mathrm{v}})}\times\langle\overline{{\cal O}}_{\rm
us}(x_{\mathrm{f}})V_{\mu}(x_{\mathrm{v}}){\cal O}_{\rm
bs}(x_{\mathrm{i}})\rangle,\\\ &{t_{\rm
K}}=(x_{\mathrm{f}}-x_{\mathrm{v}})_{0}>0\,,\quad{t_{{\textnormal{B}_{\textnormal{s}}}}}=(x_{\mathrm{v}}-x_{\mathrm{i}})_{0}>0\,,\end{split}$
(2.5)

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

$V_{\mu}(x)=\overline{\psi}_{\rm u}(x)\gamma_{\mu}\psi_{\rm b}(x)\,.$ (2.6)

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

$\displaystyle\mathcal{C}^{\rm K}(t;{\bf p}_{\rm K})$
$\displaystyle=\frac{a^{6}}{L^{3}}\sum_{{\bf x}_{\mathrm{f}},{\bf
x}_{\mathrm{i}}}e^{-i{\bf p}_{\rm K}\,({\bf x}_{\mathrm{f}}-{\bf
x}_{\mathrm{i}})}\langle\overline{{\cal O}}_{\rm us}(x_{\mathrm{f}}){\cal
O}_{\rm us}(x_{\mathrm{i}})\rangle,\quad
t=(x_{\mathrm{f}}-x_{\mathrm{i}})_{0}$ (2.7)
$\displaystyle\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}(t)$
$\displaystyle=\frac{a^{6}}{L^{3}}\sum_{{\bf x}_{\mathrm{f}},{\bf
x}_{\mathrm{i}}}\langle\overline{{\cal O}}_{\rm bs}(x_{\mathrm{f}}){\cal
O}_{\rm bs}(x_{\mathrm{i}})\rangle,\quad
t=(x_{\mathrm{f}}-x_{\mathrm{i}})_{0}\,.$ (2.8)

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

$\displaystyle\langle{\rm K}({\bf
p})|\hat{V}_{\mathrm{M}}^{0}(0)|{\textnormal{B}_{\textnormal{s}}}({\bf
0})\rangle$ $\displaystyle=$ $\displaystyle\lim_{t\to\infty}{{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{0}(t,t;{\bf
p})\over[\mathcal{C}^{\rm K}(2t;{\bf
p})\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}(2t)]^{1/2}}\,,$
$\displaystyle\langle{\rm K}({\bf
p})|\hat{V}_{\mathrm{M}}^{k}(0)|{\textnormal{B}_{\textnormal{s}}}({\bf
0})\rangle$ $\displaystyle=$ $\displaystyle i\,\lim_{t\to\infty}{{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{k}(t,t;{\bf
p})\over[\mathcal{C}^{\rm K}(2t;{\bf
p})\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}(2t)]^{1/2}}\,.$

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]

$\displaystyle\log({\cal C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K}}_{\mu})=\log Z^{\mathrm{HQET}}_{V_{\mu}}+\log([{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K}}_{\mu}]_{\mathrm{stat}})+\sum_{k}\omega_{k}\rho_{\mu}^{k}+\delta_{\cal
O}+\mathcal{O}\left(1/m_{\mathrm{b}}^{2}\right)\,,$ (2.10)

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

$\displaystyle\rho_{\mu}^{k}={[{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu}]_{k}\over[{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu}]_{\mathrm{stat}}}\,.$
(2.11)

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

ensemble | $\beta$ | $L/a$ |  $a$ [fm] |  $m_{\pi}$ [MeV] | $m_{\pi}L$
---|---|---|---|---|---
A5 | $5.2$ | $32$ | $0.0749(8)$ | $330$ | $4.0$
O7 | $5.5$ | $64$ | $0.0483(4)$ | $270$ | $4.2$
Table 1: Overview of the subset of $N_{\rm f}=2$ CLS ensembles used. They have $T=2\times L$. Lattice spacings, $a$, are taken from Ref. [19] and Ref. [20]. The pion mass is denoted by $m_{\pi}$, and we have $m_{\pi}L\geq 4$ such that finite-volume effects are sufficiently (exponentially) suppressed at the level of accuracy we are aiming for.  ensemble | $\kappa_{s}$ | $N_{\rm cfg}$ | $I^{\rm K}$ | $(I^{\textnormal{B}_{\textnormal{s}}}_{1},I^{\textnormal{B}_{\textnormal{s}}}_{2},I^{\textnormal{B}_{\textnormal{s}}}_{3})$ | ${\bf n}$ | $\vartheta$
---|---|---|---|---|---|---
A5 | $0.13535(2)$ | $1000$ | 35 | $(30,60,155)$ | $(1,0,0)^{T}$ | $+0.2141$
O7 | $0.13627(1)$ | $490$ | 75 | $(68,135,338)$ | $(1,1,0)^{T}$ | $-0.3593$
Table 2: Details of the measurements. The hopping parameter of the valence
strange quark, $\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_{\rm cfg}$. The number of smearing iterations applied to the
light quark(s) in the ${\rm K}$ and ${\textnormal{B}_{\textnormal{s}}}$ meson
is given by $I^{\rm K}$ and $I^{\textnormal{B}_{\textnormal{s}}}_{r}$,
respectively. The integer vector ${\bf n}$ and the scalar twist angle
$\vartheta$ determine the momentum of the ${\rm K}$ meson by setting
$\theta^{\mathrm{s}}={\bf n}\vartheta$ in (2.19), i.e. ${\bf p}_{\rm
K}=(2\pi-\vartheta){\bf n}/L$.

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

$\displaystyle V_{0}^{\textnormal{stat}}$
$\displaystyle=\overline{\psi}_{\textnormal{u}}\gamma_{0}\psi_{\textnormal{h}}+ac_{\textnormal{V}_{0}}(g_{0})\overline{\psi}_{\textnormal{l}}\sum_{l}\overleftarrow{\nabla}^{\textnormal{S}}_{l}\gamma_{l}\psi_{\textnormal{h}},$
(2.12a) $\displaystyle V_{k}^{\textnormal{stat}}$
$\displaystyle=\overline{\psi}_{\textnormal{u}}\gamma_{k}\psi_{\textnormal{h}}-ac_{\textnormal{V}_{k}}(g_{0})\overline{\psi}_{\textnormal{l}}\sum_{l}\overleftarrow{\nabla}^{\textnormal{S}}_{l}\gamma_{l}\gamma_{k}\psi_{\textnormal{h}}\,.$
(2.12b)

with [27]

$c_{\textnormal{V}_{k}}=0.0518g_{0}^{2}\,,\quad
c_{\textnormal{V}_{0}}=0.0380g_{0}^{2}\,$ (2.13)

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$)

$\displaystyle{\cal O}_{\mathrm{bs},r}(x)$
$\displaystyle\equiv\overline{\psi}_{\mathrm{h}}(x)\gamma_{5}\psi^{\mathrm{sm}}_{\mathrm{s},r}(x)\,,$
(2.14a) $\displaystyle\overline{{\cal O}}_{\mathrm{bs},r}(x)$
$\displaystyle\equiv\overline{\psi}^{\mathrm{sm}}_{\mathrm{s},r}(x)\gamma_{0}\gamma_{5}\gamma_{0}\psi_{\mathrm{h}}(x)\,.$
(2.14b)

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

$\psi^{\mathrm{sm}}_{\mathrm{s},r}=(1-\kappa_{\mathrm{G}}\,a^{2}\Delta)^{I^{\textnormal{B}_{\textnormal{s}}}_{r}}\psi_{s}\,.$
(2.15)

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,

$\displaystyle{\cal O}_{\mathrm{su}}(x_{\mathrm{f}})$
$\displaystyle\equiv\overline{\psi}^{\mathrm{sm}}_{\mathrm{s}}(x_{\mathrm{f}})\gamma_{5}\psi^{\mathrm{sm}}_{\mathrm{u}}(x_{\mathrm{f}})\,,$
(2.16a) $\displaystyle\overline{{\cal O}}_{\mathrm{su}}(x_{\mathrm{f}})$
$\displaystyle\equiv\overline{\psi}^{\mathrm{sm}}_{\mathrm{u}}(x_{\mathrm{f}})\gamma_{0}\gamma_{5}\gamma_{0}\psi^{\mathrm{sm}}_{\mathrm{s}}(x_{\mathrm{f}})\,,$
(2.16b)

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

$q^{2}=(m^{\rm phys}_{\textnormal{B}_{\textnormal{s}}})^{2}+(m^{\rm phys}_{\rm
K})^{2}-2m^{\rm phys}_{\textnormal{B}_{\textnormal{s}}}\sqrt{(m^{\rm
phys}_{\rm K})^{2}+{\bf p}_{\rm K}^{2}}$ (2.17)

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

$\displaystyle\psi_{\mathrm{s}}(x_{0},{\bf x}+\hat{{\bf k}}L)$
$\displaystyle=\,e^{i\theta_{k}^{\mathrm{s}}}\,\psi_{\mathrm{s}}(x_{0},{\bf
x})\,,$ (2.18a) $\displaystyle\psi_{\mathrm{b}}(x_{0},{\bf x}+\hat{{\bf k}}L)$
$\displaystyle=\,e^{i\theta_{k}^{\mathrm{b}}}\,\psi_{\mathrm{b}}(x_{0},{\bf
x})\,,$ (2.18b)

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

$(p_{\rm K})_{k}=(2\pi n_{k}-\theta^{\mathrm{s}}_{k})/L\,,$ (2.19)

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):

$\displaystyle\mathcal{C}^{\rm K}({t_{\rm K}})=$
$\displaystyle\sum_{m=0}^{\infty}(\kappa^{(m)})^{2}e^{-E_{\rm K}^{(m)}{t_{\rm
K}}}\cong\sum_{m=0}^{N_{\rm K}-1}(\kappa^{(m)})^{2}e^{-E_{\rm K}^{(m)}{t_{\rm
K}}},$ (3.20)
$\displaystyle\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}_{rs}({t_{{\textnormal{B}_{\textnormal{s}}}}})=$
$\displaystyle\sum_{n=0}^{\infty}\beta_{r}^{(n)}\beta_{s}^{(n)}e^{-E_{{\textnormal{B}_{\textnormal{s}}}}^{(n)}{t_{{\textnormal{B}_{\textnormal{s}}}}}}\cong\sum_{n=0}^{N_{\textnormal{B}_{\textnormal{s}}}-1}\beta_{r}^{(n)}\beta_{s}^{(n)}e^{-E_{{\textnormal{B}_{\textnormal{s}}}}^{(n)}{t_{{\textnormal{B}_{\textnormal{s}}}}}},$
(3.21) $\displaystyle\begin{split}{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu,r}({t_{\rm
K}},{t_{{\textnormal{B}_{\textnormal{s}}}}})=&\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\kappa^{(m)}\varphi^{(m,n)}_{\mu}\beta_{r}^{(n)}e^{-E_{\rm
K}^{(m)}t_{\rm
K}}e^{-E_{\textnormal{B}_{\textnormal{s}}}^{(n)}t_{\textnormal{B}_{\textnormal{s}}}}\cong\\\
&\sum_{n=0}^{N_{\rm
K}-1}\sum_{m=0}^{N_{\textnormal{B}_{\textnormal{s}}}-1}\kappa^{(m)}\varphi^{(m,n)}_{\mu}\beta_{r}^{(n)}e^{-E_{\rm
K}^{(m)}t_{\rm
K}}e^{-E_{\textnormal{B}_{\textnormal{s}}}^{(n)}t_{\textnormal{B}_{\textnormal{s}}}},\end{split}$
(3.22)

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.

$\displaystyle\langle{\rm K}({\bf p})|\hat{V}_{\rm
M}^{0}(0)|{\textnormal{B}_{\textnormal{s}}}({\bf 0})\rangle=$
$\displaystyle~{}\varphi^{(0,0)}_{0}.$ (3.23)

In the following we describe two methods to determine them:

  1. 1.

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

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

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.

Figure 2: Effective masses, $E_{\rm K,eff}(t_{\rm K})=-\partial_{t}\log{\cal
C}^{{\rm K}}(t_{\rm K})$ for ensemble A5 and two values of momentum, together
with the fitted energies. In this graph as well as all following ones we use
MeV and fm units just for illustration. We thus convert from lattice units to
physical ones without taking into account the error in the lattice spacings.
The lattice parameters for the Kaon with momentum $|{\bf p}|=535$ MeV are
given in Table 2. For the Kaon at $|{\bf p}|=757$ MeV, we have used ${\bf
n}=(1,0,0)$ and $\vartheta=-0.9263$. Figure 3: Kaon energies at three values
of momentum for ensemble A5. The lattice parameters are indicated in the
caption of Fig 2. The dashed curve shows the continuum dispersion relation
evaluated with the rest mass from Ref. [19] (it is not a fit).

###  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

$\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}(t)v^{(n)}(t,t_{0})=\lambda^{(n)}(t,t_{0})\,\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}(t_{0})\,v^{(n)}(t,t_{0})\,,\quad
n=1,\ldots,N_{\textnormal{B}_{\textnormal{s}}}\,,$ (3.25)

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

$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{)}.$
(3.26)

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.

Figure 4: Dependence of the ${\rm K}$ and ${\textnormal{B}_{\textnormal{s}}}$
ground-state energies on $t$ for ensemble A5 (left) and O7 (right).
$E_{{\textnormal{B}_{\textnormal{s}}},{\rm eff}}^{(0)}$ is shifted vertically
for presentation. Dashed (and full) horizontal lines show the plateau fit
values (and ranges).

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

${\cal
C}_{rr}^{\textnormal{B}_{\textnormal{s}}}(t)\cong\sum_{n=0}^{N_{\textnormal{B}_{\textnormal{s}}}-1}(\beta_{r}^{(n)})^{2}e^{-E_{\rm{\textnormal{B}_{\textnormal{s}}}}^{(n)}t}\,$
(3.27)

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

Figure 5: Contributions to the three-point function corresponding to the
${\textnormal{B}_{\textnormal{s}}}\to{\rm K}$ decay (left), and a wrap-around
Kaon (right). Figures taken from Ref. [32].

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

$\displaystyle{\cal C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K}}_{\mu,r}(t_{\rm K},t_{\textnormal{B}_{\textnormal{s}}})\cong$
$\displaystyle\kappa^{(0)}\varphi^{(0,0)}_{\mu}\beta_{r}^{(0)}e^{-E_{\rm{\textnormal{B}_{\textnormal{s}}}}^{(0)}t_{\textnormal{B}_{\textnormal{s}}}}e^{-E_{\rm
K}^{(0)}t_{\rm K}}+$ (3.28) $\displaystyle\kappa^{(0)}\xi_{\mu,r}e^{-E_{{\rm
B}^{*}}t_{\textnormal{B}_{\textnormal{s}}}}e^{-E_{\rm
K}^{(0)}(T-t_{\textnormal{B}_{\textnormal{s}}}-t_{\rm K})},$

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

${\cal C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu,r}(t_{\rm
K})\cong B_{\mu,r}e^{-Et_{\rm K}}+C_{\mu,r}e^{+Et_{\rm K}},$ (3.29)

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:

$\frac{C_{\mu,r}e^{+Et_{\rm K}}}{B_{\mu,r}e^{-Et_{\rm K}}+C_{\mu,r}e^{+Et_{\rm
K}}}<c_{\rm wr}\,\frac{\delta{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu,r}(t_{\rm
K},t_{\textnormal{B}_{\textnormal{s}}})}{{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu,r}(t_{\rm
K},t_{\textnormal{B}_{\textnormal{s}}})}\,,$ (3.30)

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

Figure 6: The effective mass $am^{\rm 3pt}_{\rm eff}({t_{\rm
K}})=\log(C^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{1,3}({t_{\rm
K}},t_{\textnormal{B}_{\textnormal{s}}})/(C^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K}}_{1,3}({t_{\rm K}}+a,t_{\textnormal{B}_{\textnormal{s}}}))$ with
$t_{\textnormal{B}_{\textnormal{s}}}=8a\approx 0.6\,{\rm fm}$ for ensemble A5,
together with the fitted function. The start of the finite-$T$ fit and $t_{\rm
K3}^{\rm max,wr}$ are also shown.

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.

Figure 7: Time ranges $t_{\rm K3}^{\rm max,wr}$, $t_{\rm K3}^{\rm max,noise}$,
and $t_{\rm K3}^{\rm min}$ for lattices A5 and O7. The different lines with
same color show the values for the three smearings. The vertical dashed lines
are the two choices of $t_{\rm B3}^{\rm min}$ for the combined fit of Sec. 4
and the dots correspond to the data points included in the fit. To guide the
eye, we also show the line ${t_{\rm
K}}={t_{{\textnormal{B}_{\textnormal{s}}}}}$ which is relevant for the ratios
in Sec. 5.

##  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

$t_{\rm K3}^{\rm max}(t_{\textnormal{B}_{\textnormal{s}}})=\min\big{\\{}t_{\rm
K3}^{\rm max,wr},\,t_{\rm K3}^{\rm max,noise}\big{\\}}\,,$ (4.31)

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

$\frac{\delta{\cal C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K}}_{\mu,r}(t_{\rm K},t_{\textnormal{B}_{\textnormal{s}}})}{{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu,r}(t_{\rm
K},t_{\textnormal{B}_{\textnormal{s}}})}<c_{\rm noise}\,,$ (4.32)

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

$t_{\rm K3}^{\rm
min}(t_{\textnormal{B}_{\textnormal{s}}})=\max\big{\\{}0.8\;{\rm fm},\,t_{\rm
K3*}^{\rm min}(t_{\textnormal{B}_{\textnormal{s}}})\big{\\}}$ (4.33)

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.

Figure 8: Stability plots for all three ${\textnormal{B}_{\textnormal{s}}}$
energies on ensemble A5 (left) and O7 (right). The selected $t_{\rm B2/3}^{\rm
min}$ values are highlighted with a black color, while the dashed lines are
values from the GEVP, cf. Sec. 3.2.1. The behaviour at large $t_{\rm B2}^{\rm
min}$ and $t_{\rm B3}^{\rm min}$ can be traced to the fact that the fit can no
longer resolve three separate states. The selected fit values are discussed in
Sec. 4.1. A zoom into the behaviour of
$E_{\rm{\textnormal{B}_{\textnormal{s}}}}^{(0)}$ is shown in Fig. 9.

Figure 9: Stability plots for $E_{\rm{\textnormal{B}_{\textnormal{s}}}}^{(0)}$
for ensemble A5 (left) and O7 (right).

Figure 10: Stability plots of the ground-state form factors for ensemble A5
(left) and O7 (right). The results are rather stable with respect to changes
of $t_{\rm B2/3}^{\rm min}$, although errors grow quickly.

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.

Figure 11: Excited-state form factors at $\mu=1$ for ensemble O7. id | $\mu$ | Fit | Ratio | Summed ratio
---|---|---|---|---
|  | $t_{\rm B2}^{\rm min}\approx 0.32\,{\rm fm}$ | $t_{\rm B2}^{\rm min}\approx 0.38\,{\rm fm}$ | $\mathcal{R}^{\,\textrm{II}}$ | $\mathcal{M}^{\textrm{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 3: Results for the ground-state form factors $\varphi_{\mu}^{(0,0)}$.
The second error in the result for $\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
$\mathcal{S}^{\textrm{I}}_{\mu,3}(\tau)$, see the shaded band in Fig. (14).

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}}}}}$):

$\mathcal{R}^{\,\textrm{I}}_{\mu,r}({t_{\rm
K}},{t_{{\textnormal{B}_{\textnormal{s}}}}})=\frac{{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu,r}({t_{\rm
K}},{t_{{\textnormal{B}_{\textnormal{s}}}}})}{\big{[}{\cal C}^{\rm
K}(\tau){\cal
C}^{\textnormal{B}_{\textnormal{s}}}_{rr}(\tau)\big{]}^{1/2}}\exp\big{\\{}(\tilde{E}^{(0)}_{\textnormal{B}_{\textnormal{s}}}(\tau)-\tilde{E}^{(0)}_{\rm
K}(\tau))\tfrac{{t_{{\textnormal{B}_{\textnormal{s}}}}}-{t_{\rm
K}}}{2}\big{\\}}\,.$ (5.34)

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

$\mathcal{R}^{\,\textrm{II}}_{\mu,r}({t_{\rm
K}},{t_{{\textnormal{B}_{\textnormal{s}}}}})=\frac{{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu,r}({t_{\rm
K}},{t_{{\textnormal{B}_{\textnormal{s}}}}})}{\big{[}{\cal C}^{\rm K}({t_{\rm
K}}){\cal
C}^{\textnormal{B}_{\textnormal{s}}}_{rr}({t_{{\textnormal{B}_{\textnormal{s}}}}})\big{]}^{1/2}}\exp\big{\\{}\tilde{E}^{(0)}_{\textnormal{B}_{\textnormal{s}}}(\tau)\tfrac{{t_{{\textnormal{B}_{\textnormal{s}}}}}}{2}+\tilde{E}^{(0)}_{\rm
K}(\tau)\tfrac{{t_{\rm K}}}{2}\big{\\}}$ (5.35)

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

$\mathcal{R}^{\,\textrm{III}}_{\mu,i}({t_{\rm
K}},{t_{{\textnormal{B}_{\textnormal{s}}}}})=\frac{{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu,r}({t_{\rm
K}},{t_{{\textnormal{B}_{\textnormal{s}}}}})}{\mathcal{N}^{\rm K}{\cal C}^{\rm
K}({t_{\rm K}})\,\mathcal{N}^{\textnormal{B}_{\textnormal{s}}}_{r}{\cal
C}^{\textnormal{B}_{\textnormal{s}}}_{rr}({t_{{\textnormal{B}_{\textnormal{s}}}}})},$
(5.36)

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.

Figure 12: Overview of the ratios $\mathcal{R}^{\,\textrm{I}}$ with ${t_{\rm
K}}={t_{{\textnormal{B}_{\textnormal{s}}}}}$ for ensembles A5 (left) and O7
(right). Only the highest smearing, $r=3$, is presented. For comparison we
also show the results from the combined fit in Sec. 4: the dashed horizontal
lines are the fitted values of $\varphi^{(0,0)}_{\mu}$ and the curves
($1\sigma$ bands) are the respective ratios of correlation functions
reconstructed from the fitted parameters.

Figure 13: Overview of the ratios $\mathcal{R}^{\,\textrm{II}}$ with ${t_{\rm
K}}={t_{{\textnormal{B}_{\textnormal{s}}}}}$ for ensemble A5 (left) and O7
(right), with the fit curves analogous as in Fig. 12. Selected plateaux are
also shown. For $\mu=1$ on ensemble A5 the wrapper contribution limits the
available time range to $t\leq t_{\rm K3}^{\rm max,wr}$, cf. Fig. 7.

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

$\mathcal{M}^{\textrm{I}}_{\mu,r}(\tau)=\partial_{\tau}\,\mathcal{S}^{\textrm{I}}_{\mu,r}(\tau)=\partial_{\tau}\,a\sum_{{t_{{\textnormal{B}_{\textnormal{s}}}}}=0}^{\tau}\mathcal{R}^{\,\textrm{I}}_{\mu,r}(\tau-{t_{{\textnormal{B}_{\textnormal{s}}}}},{t_{{\textnormal{B}_{\textnormal{s}}}}})\,.$
(5.37)

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.

Figure 14: Overview of the $\mathcal{M}^{\textrm{I}}$ for ensemble A5 (left)
and O7 (right), together with the fit curves analogous as in Figs. 12 and 13.
The shaded bands represents $c_{2}$ of fits
$\mathcal{S}^{\textrm{I}}_{\mu,r}(\tau)=c_{1}+c_{2}\tau$ to the sum in the
indicated range of $\tau$, chosen such that wrapper effects are negligible.

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

$\displaystyle\log(\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}(t))$
$\displaystyle=-m_{\textnormal{bare}}t+\log({{\cal
C}}^{{\textnormal{B}_{\textnormal{s}}},{\rm
stat}}(t))+\sum_{k\in\\{\mathrm{kin},\mathrm{spin}\\}}\omega_{k}\left(\frac{{{\cal
C}}^{{\textnormal{B}_{\textnormal{s}}},k}(t)}{{{\cal
C}}^{{\textnormal{B}_{\textnormal{s}}},{\rm stat}}(t)}\right).$ (6.38)

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]

$\displaystyle-\partial_{t}\log{{\cal
C}}^{{\textnormal{B}_{\textnormal{s}}},{\rm stat}}_{rr}(t)=E^{{\rm
stat}}+\mathcal{O}\left(e^{-\Delta E^{{\rm stat}}t}\right),$ (6.39)
$\displaystyle-\partial_{t}\frac{{{\cal
C}}^{{\textnormal{B}_{\textnormal{s}}},\textnormal{k}}_{rr}(t)}{{{\cal
C}}^{{\textnormal{B}_{\textnormal{s}}},{\rm
stat}}_{rr}(t)}=E^{\textnormal{k}}+\mathcal{O}\left(te^{-\Delta E^{{\rm
stat}}t}\right)\,,\quad k\in\\{\mathrm{kin},\mathrm{spin}\\}.$ (6.40)

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:

$\displaystyle\frac{{{\cal
C}}^{{\textnormal{B}_{\textnormal{s}}},\textnormal{k}}_{rr}(t)}{{{\cal
C}}^{{\textnormal{B}_{\textnormal{s}}},{\rm stat}}_{rr}(t)}={\rm
A}^{{\textnormal{B}_{\textnormal{s}}},\textnormal{k}}_{r}-E^{\textnormal{k}}t+\mathcal{O}\left(te^{-\Delta
E^{{\rm stat}}t}\right)\,.$ (6.41)

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}$):

$\displaystyle\frac{{{\cal C}}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K},k}_{\mu,r}(t,t)}{{{\cal C}}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K},{\rm stat}}_{\mu,r}(t,t)}={\rm A}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K},k}_{\mu,r}-E^{k}t+\mathcal{O}\left(te^{-\Delta E^{{\rm stat}}t}\right).$
(6.42)

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

$\displaystyle\mathcal{R}^{\,\textrm{I}}_{\mu,r}(t)$
$\displaystyle=\mathcal{R}^{\,\textrm{I},{{\rm
stat}}}_{\mu,r}(t)\left[1+\sum_{k}\omega_{k}\,\rho^{\textnormal{I},{k}}_{\mu,r}(t)+\sum_{j}\omega_{\mu
j}\,\rho^{\textnormal{I},{j}}_{\mu,r}(t)\right]$

with

$\displaystyle\rho^{\textnormal{I},{k}}_{\mu,r}(t,t)$
$\displaystyle=\frac{{{\cal C}}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K},k}_{\mu,r}(t,t)}{{{\cal C}}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K},{\rm stat}}_{\mu,r}(t,t)}-\frac{1}{2}\frac{{{\cal
C}}^{{\textnormal{B}_{\textnormal{s}}},k}_{r}(2t)}{{{\cal
C}}^{{\textnormal{B}_{\textnormal{s}}},{\rm stat}}_{r}(2t)}\,,\quad
k\in\\{\mathrm{kin},\mathrm{spin}\\}\,,$ (6.43)
$\displaystyle\rho^{\textnormal{I},{j}}_{\mu,r}(t,t)$
$\displaystyle=\frac{{{\cal C}}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K},j}_{\mu,r}(t,t)}{{{\cal C}}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K},j}_{\mu,r}(t,t)}\,.$ (6.44)

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,

$\displaystyle\rho^{k}_{\mu}$
$\displaystyle=\lim_{t\to\infty}\rho^{\textnormal{I},{k}}_{\mu,r}(t)={\rm
A}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K},k}_{\mu,r}-\frac{1}{2}{\rm
A}^{{\textnormal{B}_{\textnormal{s}}},k}_{r}\,.$ (6.45)
$\displaystyle\rho^{j}_{\mu}$
$\displaystyle=\lim_{t\to\infty}\rho^{\textnormal{I},{j}}_{\mu,r}(t)\,,$
(6.46)

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

Figure 15: $a^{2}E_{\rm kin}$ for ensembles A5 (left) and O7 (right). The
highest smearing and the GEVP result are shown. The band shows the GEVP
plateaux average.

Figure 16: $a^{2}E_{\rm spin}$ for ensembles A5 (left) and O7 (right). The
highest smearing and the GEVP result are shown. The band shows the GEVP
plateaux average.

###  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

$E_{\rm
eff}^{k}(t,t_{0})=-{\partial}_{t}\log\big{(}\lambda^{k}(t,t_{0})\big{)}\,.$
(6.47)

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.

$\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 4: Matrix elements at $\mathcal{O}\left(1/m\right)$: kin and spin
contributions.

Figure 17: Fitted ${\rm A}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K},\mathrm{kin}}_{\mu,3}$ and ${\rm
A}^{{\textnormal{B}_{\textnormal{s}}},\mathrm{kin}}_{3}$ for ensembles A5
(left) and O7 (right). See text for explanations.

Figure 18: Absolute values of fitted ${\rm
A}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K},\mathrm{spin}}_{\mu,3}$ and
${\rm A}^{{\textnormal{B}_{\textnormal{s}}},\mathrm{spin}}_{3}$ for ensembles
A5 (left) and O7 (right). See text for explanations.

Figure 19: Overview of the current insertions for ensembles A5 (left) and O7
(right). Fit bands are plateaux averages starting at 0.86 fm. For A5 the
$V_{1,4}$ contribution is exactly equal to $V_{1,3}$, while on the O7 they are
not exactly identical but still indistinguishable on the plot, therefore
$V_{1,4}$ was not plotted.

###  6.2 Current insertions at $1/m$ order

The $1/m$ vector current contributions are

$\displaystyle
V_{0}^{\mathrm{HQET}}(x)=Z_{V_{0}}^{\mathrm{HQET}}\big{(}V_{0}^{\mathrm{stat}}(x)+{\textstyle\sum}_{j=1}^{2}\omega_{0,j}V_{0,j}(x)\big{)},$
(6.48) $\displaystyle
V_{i}^{\mathrm{HQET}}(x)=Z_{V_{i}}^{\mathrm{HQET}}\big{(}V_{i}^{\mathrm{stat}}(x)+{\textstyle\sum}_{j=1}^{4}\omega_{i,j}V_{i,j}(x)\big{)}\,,\quad
V_{\mu,j}(x)=\overline{\psi}_{u}\Gamma_{\mu,j}\psi_{\rm h}\,,$ (6.49)

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.

$\mu$ | $j$ | $\Gamma_{\mu}^{j}$ | $\omega^{\textnormal{tree}}_{\mu,j}\cdot m_{\rm h}$ |  $a\rho_{\mu}^{j}$ A5 |  $a\rho_{\mu}^{j}$ O7
---|---|---|---|---|---
0 | 1 | $\sum_{l}\gamma_{l}\tfrac{1}{2}(\nabla_{l}^{\textnormal{S}}-\overleftarrow{\nabla}_{l}^{\textnormal{S}})$ | 1/2 | -0.0182(12) | -0.0241(07)
0 | 2 | $\sum_{l}\gamma_{l}\tfrac{1}{2}(\nabla_{l}^{\textnormal{S}}+\overleftarrow{\nabla}_{l}^{\textnormal{S}})$ | 1/2 | -0.0434$(\phantom{0}4)$ | -0.0289(02)
$i$ | 1 | $\sum_{l}\tfrac{1}{2}(\nabla_{l}^{\textnormal{S}}-\overleftarrow{\nabla}_{l}^{\textnormal{S}})\gamma_{l}\gamma_{i}$ | 1/2 | 0.3508(27) | 0.2852(15)
$i$ | 2 | $\tfrac{1}{2}(\nabla_{i}^{\textnormal{S}}-\overleftarrow{\nabla}_{i}^{\textnormal{S}})$ | -1 | 0.0759(17) | 0.0709(11)
$i$ | 3 | $\sum_{l}\tfrac{1}{2}(\nabla_{l}^{\textnormal{S}}+\overleftarrow{\nabla}_{l}^{\textnormal{S}})\gamma_{l}\gamma_{i}$ | 1/2 | 0.1737(16) | 0.1132(08)
$i$ | 4 | $\tfrac{1}{2}(\nabla_{i}^{\textnormal{S}}+\overleftarrow{\nabla}_{i}^{\textnormal{S}})$ | -1 | 0.1737(16) | 0.1133(08)
Table 5: Overview of the $1/m$ vector current insertions, their tree-level
matching coefficients, and the results extracted from the highest light-quark
smearing, $r=3$. We use symmetric covariant derivatives
$\nabla^{\textnormal{S}}_{i}$.

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

${\cal C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm K}}_{\mu,r}({t_{\rm
K}},{t_{{\textnormal{B}_{\textnormal{s}}}}};{\bf p})=\langle\widehat{\cal
C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K}}_{\mu,r}\rangle_{U}\,,\quad\mathcal{C}^{\rm K}(t;{\bf
p})=\langle\widehat{\cal C}^{\rm
K}\rangle_{U}\,,\quad\mathcal{C}^{\textnormal{B}_{\textnormal{s}}}_{rr^{\prime}}(t)=\langle\widehat{\cal
C}^{\textnormal{B}_{\textnormal{s}}}_{rr^{\prime}}\rangle_{U}\ $ (A.50)

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

$\displaystyle\widehat{\cal C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K}}_{\mu,r}$
$\displaystyle=\frac{a^{9}}{L^{3}}\mathrm{Tr}\left[\gamma_{5}P_{\mathrm{f}}F_{-}W_{\rm
K}S_{\mathrm{u}}\gamma_{\mu}P_{\mathrm{v}}F_{+}S_{\mathrm{h}}P_{\mathrm{i}}\gamma_{5}W_{r}S_{\mathrm{s}}W_{\rm
K}\right]$ (A.51) $\displaystyle\widehat{\cal C}^{\rm K}$
$\displaystyle=\frac{a^{6}}{L^{3}}\mathrm{Tr}\left[\gamma_{5}P_{\mathrm{f}}F_{-}W_{\rm
K}S_{\mathrm{u}}W_{\rm K}\gamma_{5}P_{\mathrm{i}}F_{+}W_{\rm
K}S_{\mathrm{s}}W_{\rm K}\right]$ (A.52) $\displaystyle\widehat{\cal
C}^{\textnormal{B}_{\textnormal{s}}}_{rr^{\prime}}$
$\displaystyle=\frac{a^{6}}{L^{3}}\mathrm{Tr}\left[\gamma_{5}P_{\mathrm{f}}S_{\mathrm{h}}P_{\mathrm{i}}\gamma_{5}W_{r^{\prime}}S_{\mathrm{s}}W_{r}\right]$
(A.53)

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

$P_{\mathrm{f}}=\langle\eta\eta^{\dagger}\rangle_{\eta}$ (A.54)

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

$\displaystyle\widehat{\cal C}^{{\textnormal{B}_{\textnormal{s}}}\to{\rm
K}}_{\mu,r}$
$\displaystyle=\frac{a^{9}}{L^{3}}\left\langle\eta^{\dagger}F_{-}\gamma_{5}W_{\rm
K}S_{\mathrm{u}}P_{\mathrm{v}}\gamma_{\mu}F_{+}S_{\mathrm{h}}P_{\mathrm{i}}\gamma_{5}W_{r}S_{\mathrm{s}}W_{\rm
K}\,\eta\right\rangle_{\eta}$ (A.55)
$\displaystyle=\left\langle\phi_{\mathrm{u}}^{\dagger}\gamma_{5}P_{\mathrm{v}}\gamma_{\mu}F_{+}S_{\mathrm{h}}P_{\mathrm{i}}\gamma_{5}W_{r}\phi_{\mathrm{s}}\right\rangle_{\eta}$
(A.56) $\displaystyle\phi_{\mathrm{s}}=S_{\mathrm{s}}W_{\rm
K}\,\eta\,,\quad\phi_{\mathrm{u}}=S_{\mathrm{u}}W_{\rm
K}F_{-}^{\dagger}\,\eta\,$ (A.57)

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

$\displaystyle S_{\mathrm{kin}}$
$\displaystyle=S_{\mathrm{h}}\nabla_{i}^{*}\nabla_{i}S_{\mathrm{h}}$ (A.58)
$\displaystyle S_{\mathrm{spin}}$
$\displaystyle=\frac{i}{4}S_{\mathrm{h}}[\gamma_{i},\gamma_{j}]\hat{F}_{ij}S_{\mathrm{h}}$
(A.59)

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

Figure 20: Stability plots for the ${\textnormal{B}_{\textnormal{s}}}$
energies on ensemble A5 (left) and O7 (right). The selected $t_{\rm B2/3}^{\rm
min}$ values are highlighted with a black color, while the dashed lines are
values from the GEVP, cf. Sec. 3.2.1. A zoom into the behaviour of
$E_{\rm{\textnormal{B}_{\textnormal{s}}}}^{(0)}$ is shown in Fig. 21.

Figure 21: Stability plots for
$E_{\rm{\textnormal{B}_{\textnormal{s}}}}^{(0)}$ for ensemble A5 (left) and O7
(right).

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.

id | $\mu$ | Fit value
---|---|---
A5 | 0 | 1.099(12)
A5 | 1 | 0.580(9)
O7 | 0 | 1.100(13)
O7 | 1 | 0.606(9)
Table 6: Results for the ground-state form factors $\varphi_{\mu}^{(0)}$ using
the fit with $N_{\textnormal{B}_{\textnormal{s}}}=2$.

Figure 22: Stability plots of the ground-state form factors for ensemble A5
(left) and O7 (right). Figure 23: The summed ratio for ensemble O7, $\mu=0$
with three versions: with symmetric sum limits and with either only
$t_{\textnormal{B}_{\textnormal{s}}}^{\rm min}$ or $t_{\rm K}^{\rm min}$
larger than 0\. Note that due to different (worse) convergence properties,
these plots cannot be directly compared with the ones in Fig. 14.

##  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:

$\tilde{\mathcal{M}}^{\textrm{I}}_{\mu,r}(\tau)=\partial_{\tau}\,a\sum_{{t_{{\textnormal{B}_{\textnormal{s}}}}}=t_{\textnormal{B}_{\textnormal{s}}}^{\rm
min}}^{\tau-t_{\rm K}^{\rm
min}}\mathcal{R}^{\,\textrm{I}}_{\mu,r}(\tau-{t_{{\textnormal{B}_{\textnormal{s}}}}},{t_{{\textnormal{B}_{\textnormal{s}}}}}).$
(C.60)

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.

## References

  * [1] Fermilab Lattice, MILC collaboration, J. A. Bailey et al., _$|V_{ub}|$ from $B\to\pi\ell\nu$ decays and (2+1)-flavor lattice QCD_, _Phys. Rev._ D92 (2015) 014024, [1503.07839].
  * [2] J. M. Flynn, T. Izubuchi, T. Kawanai, C. Lehner, A. Soni, R. S. Van de Water et al., _$B\to\pi\ell\nu$ and $B_{s}\to K\ell\nu$ form factors and $|V_{ub}|$ from 2+1-flavor lattice QCD with domain-wall light quarks and relativistic heavy quarks_, _Phys. Rev._ D91 (2015) 074510, [1501.05373].
  * [3] C. M. Bouchard, G. P. Lepage, C. Monahan, H. Na and J. Shigemitsu, _$B_{s}\to K\ell\nu$ form factors from lattice QCD_, _Phys. Rev._ D90 (2014) 054506, [1406.2279].
  * [4] B. Colquhoun, R. J. Dowdall, J. Koponen, C. T. H. Davies and G. P. Lepage, _$B\rightarrow\pi\ell\nu$ at zero recoil from lattice QCD with physical $u/d$ quarks_, _Phys. Rev._ D93 (2016) 034502, [1510.07446].
  * [5] B. Colquhoun, S. Hashimoto and T. Kaneko, _$B\rightarrow\pi\ell\nu$ with Möbius Domain Wall Fermions_, _EPJ Web Conf._ 175 (2018) 13004, [1710.07094].
  * [6] Z. Gelzer et al., _Semileptonic $B$-meson decays to light pseudoscalar mesons on the HISQ ensembles_, _EPJ Web Conf._ 175 (2018) 13024, [1710.09442].
  * [7] C. J. Monahan, C. M. Bouchard, G. P. Lepage, H. Na and J. Shigemitsu, _Form factor ratios for $B_{s}\rightarrow K\,\ell\,\nu$ and $B_{s}\rightarrow D_{s}\,\ell\,\nu$ semileptonic decays and $|V_{ub}/V_{cb}|$_, _Phys. Rev._ D98 (2018) 114509, [1808.09285].
  * [8] Fermilab Lattice, MILC collaboration, Y. Liu et al., _$B_{s}\to K\ell\nu$ Form Factors with 2+1 Flavors_, _EPJ Web Conf._ 175 (2018) 13008, [1711.08085].
  * [9] Fermilab Lattice, MILC collaboration, A. Bazavov et al., _$B_{s}\to K\ell\nu$ decay from lattice QCD_, 1901.02561.
  * [10] F. Bahr, D. Banerjee, F. Bernardoni, A. Joseph, M. Koren, H. Simma et al., _Continuum limit of the leading-order HQET form factor in $B_{s}\to K\ell\nu$ decays_, _Phys. Lett._ B757 (2016) 473–479.
  * [11] J. Heitger and R. Sommer, _Nonperturbative heavy quark effective theory_ , _JHEP_ 0402 (2004) 022.
  * [12] B. Blossier, M. Della Morte, N. Garron and R. Sommer, _HQET at order $1/m$: I. Non-perturbative parameters in the quenched approximation_, _JHEP_ 1006 (2010) 002.
  * [13] B. Blossier, M. Della Morte, P. Fritzsch, N. Garron, J. Heitger, R. Sommer et al., _Parameters of Heavy Quark Effective Theory from $N_{\rm f}=2$ lattice QCD_, _JHEP_ 1209 (2012) 132.
  * [14] M. Della Morte, S. Dooling, J. Heitger, D. Hesse and H. Simma, _Matching of heavy-light flavor currents between HQET at order 1/ $m$ and QCD: I. Strategy and tree-level study_, _JHEP_ 1405 (2014) 060.
  * [15] D. Hesse and R. Sommer, _A one-loop study of matching conditions for static-light flavor currents_ , _JHEP_ 1302 (2013) 115.
  * [16] C. Best, M. Gockeler, R. Horsley, E.-M. Ilgenfritz, H. Perlt, P. E. L. Rakow et al., _Pion and rho structure functions from lattice QCD_ , _Phys. Rev._ D56 (1997) 2743–2754, [hep-lat/9703014].
  * [17] R. Sommer, _Introduction to Non-perturbative Heavy Quark Effective Theory_ , in _Modern perspectives in lattice QCD: Quantum field theory and high performance computing. Proceedings, International School, 93rd Session, Les Houches, France, August 3-28, 2009_ , pp. 517–590, 2010.
  * [18] B. Blossier et al., _HQET at order 1/m: III. Decay constants in the quenched approximation_ , _JHEP_ 1012 (2010) 039.
  * [19] P. Fritzsch, F. Knechtli, B. Leder, M. Marinkovic, S. Schaefer, R. Sommer et al., _The strange quark mass and Lambda parameter of two flavor QCD_ , _Nucl. Phys. B_ 865 (2012) 397.
  * [20] G. P. Engel, L. Giusti, S. Lottini and R. Sommer, _Chiral Symmetry Breaking in QCD with Two Light Flavors_ , _Phys. Rev. Lett._ 114 (2015) 112001.
  * [21] P. Fritzsch, F. Knechtli, B. Leder, M. Marinkovic, S. Schaefer et al., _The strange quark mass and the $\Lambda$ parameter of two flavor QCD_, _Nucl.Phys._ B865 (2012) 397–429, [1205.5380].
  * [22] ALPHA collaboration, K. Jansen and R. Sommer, _O(a) improvement of lattice QCD with two flavors of Wilson quarks_ , _Nucl. Phys._ B530 (1998) 185–203, [hep-lat/9803017].
  * [23] K. G. Wilson, _Confinement of quarks_ , _Phys. Rev._ D10 (1974) 2445–2459.
  * [24] ALPHA collaboration, S. Lottini, _Approaching the chiral point in two-flavour lattice simulations_ , _Acta Phys. Polon. Supp._ 7 (2014) 565, [1406.2939].
  * [25] ALPHA collaboration, M. Della Morte et al., _Lattice HQET with exponentially improved statistical precision_ , _Phys.Lett._ B581 (2004) 93–98, [hep-lat/0307021].
  * [26] M. Della Morte, A. Shindler and R. Sommer, _On lattice actions for static quarks_ , _JHEP_ 0508 (2005) 051, [hep-lat/0506008].
  * [27] A. Grimbach, D. Guazzini, F. Knechtli and F. Palombi, _O(a) improvement of the HYP static axial and vector currents at one-loop order of perturbation theory_ , _JHEP_ 03 (2008) 039, [0802.0862].
  * [28] C. Alexandrou, F. Jegerlehner, S. Gusken, K. Schilling and R. Sommer, _B meson properties from lattice QCD_ , _Phys. Lett._ B256 (1991) 60–67.
  * [29] S. Güsken, U. Löw, K. H. Mütter, R. Sommer, A. Patel and K. Schilling, _Nonsinglet Axial Vector Couplings of the Baryon Octet in Lattice QCD_ , _Phys. Lett._ B227 (1989) 266–269.
  * [30] APE collaboration, M. Albanese et al., _Glueball Masses and String Tension in Lattice QCD_ , _Phys. Lett._ B192 (1987) 163–169.
  * [31] M. Lüscher and P. Weisz, _On-Shell Improved Lattice Gauge Theories_ , _Commun. Math. Phys._ 97 (1985) 59.
  * [32] F. T. Bahr, _Form factors for semileptonic ${\textnormal{B}_{\textnormal{s}}}\to{\rm K}\ell\nu$ decays in lattice QCD_.  PhD thesis, Humboldt-Universität zu Berlin, Mathematisch-Naturwissenschaftliche Fakultät, 2015.  http://dx.doi.org/10.18452/17384.
  * [33] M. Lüscher and U. Wolff, _How to Calculate the Elastic Scattering Matrix in Two-dimensional Quantum Field Theories by Numerical Simulation_ , _Nucl. Phys. B_ 339 (1990) 222.
  * [34] B. Blossier, M. Della Morte, G. von Hippel, T. Mendes and R. Sommer, _On the generalized eigenvalue method for energies and matrix elements in lattice field theory_ , _JHEP_ 0904 (2009) 094.
  * [35] F. Bernardoni, B. Blossier, J. Bulava, M. Della Morte, P. Fritzsch et al., _The b-quark mass from non-perturbative $N_{f}=2$ Heavy Quark Effective Theory at $O(1/m_{h})$_, _Phys. Lett. B_ 730 (2014) 171.
  * [36] ALPHA collaboration, U. Wolff, _Monte Carlo errors with less errors_ , _Comput. Phys. Commun._ 156 (2004) 143–153, [hep-lat/0306017].
  * [37] ALPHA collaboration, S. Schaefer, R. Sommer and F. Virotta, _Critical slowing down and error analysis in lattice QCD simulations_ , _Nucl. Phys._ B845 (2011) 93–119, [1009.5228].
  * [38] C. Michael, _Fitting correlated data_ , _Phys. Rev._ D49 (1994) 2616–2619.
  * [39] C. Michael and A. McKerrell, _Fitting correlated hadron mass spectrum data_ , _Phys. Rev._ D51 (1995) 3745–3750.
  * [40] B. Bunk _, private communications_ (1985) .
  * [41] M. Bruno and R. Sommer _, in preparation_ .
  * [42] L. Maiani, G. Martinelli, M. L. Paciello and B. Taglienti, _Scalar Densities and Baryon Mass Differences in Lattice QCD With Wilson Fermions_ , _Nucl. Phys._ B293 (1987) 420.
  * [43] J. Bulava, M. Donnellan and R. Sommer, _On the computation of hadron-to-hadron transition matrix elements in lattice QCD_ , _JHEP_ 01 (2012) 140, [1108.3774].
  * [44] S. Capitani, B. Knippschild, M. Della Morte and H. Wittig, _Systematic errors in extracting nucleon properties from lattice QCD_ , _PoS_ LATTICE2010 (2010) 147, [1011.1358].
  * [45] F. Bernardoni, B. Blossier, J. Bulava, M. Della Morte, P. Fritzsch et al., _The b-quark mass from non-perturbative $N_{f}=2$ Heavy Quark Effective Theory at $O(1/m_{h})$_, _Phys.Lett._ B730 (2014) 171–177.
  * [46] B. Blossier, M. Della Morte, G. von Hippel, T. Mendes and R. Sommer, _On the generalized eigenvalue method for energies and matrix elements in lattice field theory_ , _JHEP_ 0904 (2009) 094, [0902.1265].
  * [47] M. Cè, L. Giusti and S. Schaefer, _A local factorization of the fermion determinant in lattice QCD_ , _Phys. Rev._ D95 (2017) 034503, [1609.02419].
  * [48] M. Cè, L. Giusti and S. Schaefer, _Local multiboson factorization of the quark determinant_ , _EPJ Web Conf._ 175 (2018) 11005, [1711.01592].
  * [49] S. Hashimoto, _Hints and challenges in heavy flavor physics_ , in _36th International Symposium on Lattice Field Theory (Lattice 2018) East Lansing, MI, United States, July 22-28, 2018_ , 2019.  1902.09119.
  * [50] Bär, Oliver and M. Golterman, _Chiral perturbation theory for gradient flow observables_ , _Phys. Rev._ D89 (2014) 034505, [1312.4999].
  * [51] Bär, Oliver, _Multi-hadron-state contamination in nucleon observables from chiral perturbation theory_ , _EPJ Web Conf._ 175 (2018) 01007, [1708.00380].
  * [52] J. Green, _Systematics in nucleon matrix element calculations_ , in _36th International Symposium on Lattice Field Theory (Lattice 2018) East Lansing, MI, United States, July 22-28, 2018_ , 2018.  1812.10574.
  * [53] O. Bär, _Nucleon-pion-state contribution in lattice calculations of the nucleon charges $g_{A},g_{T}$ and $g_{S}$_, _Phys. Rev._ D94 (2016) 054505, [1606.09385].
  * [54] Bär, Oliver, _$N\pi$ -state contamination in lattice calculations of the nucleon axial form factors_, 1812.09191.
  * [55] G. S. Bali, B. Lang, B. U. Musch and A. Schäfer, _Novel quark smearing for hadrons with high momenta in lattice QCD_ , _Phys. Rev._ D93 (2016) 094515, [1602.05525].