Lattice QCD form factor for $B_s\to D_s^* l\nu$ at zero recoil with non-perturbative current renormalisation

TL;DR

This paper reports a non-perturbative lattice QCD calculation of the $B_s o D_s^*$ axial form factor at zero recoil, achieving higher precision and reducing uncertainties compared to previous methods, with implications for CKM matrix element determinations.

Contribution

The study introduces a fully non-perturbative renormalization of the axial current in lattice QCD for $B_s o D_s^*$, improving accuracy over prior perturbative approaches.

Findings

Form factor $ ightarrow$ 0.9020(96)(90)

Reduced total uncertainty by over a factor of two

Consistent with earlier lattice results

Abstract

We present details of a lattice QCD calculation of the $B_s\to D_s^*$ axial form factor at zero recoil using the Highly Improved Staggered Quark (HISQ) formalism on the second generation MILC gluon ensembles that include up, down, strange and charm quarks in the sea. Using the HISQ action for all valence quarks means that the lattice axial vector current that couples to the $W$ can be renormalized fully non-perturbatively, giving a result free of the perturbative matching errors that previous lattice QCD calculations have had. We calculate correlation functions at three values of the lattice spacing, and multiple `$b$'-quark masses, for physical $c$ and $s$. The functional dependence on the $b$-quark mass can be determined and compared to Heavy Quark Effective Theory expectations, and a result for the form factor obtained at the physical value of the $b$-quark mass. We find $\mathcal{F}^{B_s\to D_s^*}(1) = h^s_{A_1}(1) = 0.9020(96)_{\text{stat}}(90)_{\text{sys}}$. This is in agreement with earlier lattice QCD results, which use NRQCD $b$ quarks, with a total uncertainty reduced by more than a factor of two. We discuss implications of this result for the $B\to D^*$ axial form factor at zero recoil and for determinations of $V_{cb}$.

Figures

Tables

Table 1. Table 1: Parameters for the ensembles of gluon field configurations that we use  Bazavov  et al. ( 2010b , 2013 ) . a 𝑎 a is the lattice spacing, determined from the Wilson flow parameter, w 0 subscript 𝑤 0 w_{0} . Values for w 0 / a subscript 𝑤 0 𝑎 w_{0}/a are from: set 1,  Chakraborty  et al. ( 2017 ) , sets 2 and 3,  Chakraborty  et al. ( 2015 ) and set 4  McNeile ( 2015 ) . The physical value of w 0 subscript 𝑤 0 w_{0} was determined at 0.1715(9) fm from f π subscript 𝑓 𝜋 f_{\pi}   Dowdall  et al. ( 2013 ) . N x subscript 𝑁 𝑥 N_{x} is the spatial extent and N t subscript 𝑁 𝑡 N_{t} the temporal extent of the lattice in lattice units; n cfg subscript 𝑛 cfg n_{\text{cfg}} is the number of gluon field configurations in the ensemble and n src subscript 𝑛 src n_{\text{src}} the number of different time sources used per configuration. Light, strange and charm quarks are included in the sea, their masses are given in columns 6-8, and the valence quark masses in columns 9-11. The s 𝑠 s and c 𝑐 c valence quarks were tuned in  Chakraborty  et al. ( 2015 ) . We use a number of heavy quark masses to assist the extrapolation to the physical b 𝑏 b mass. Column 12 gives the temporal separations between source and sink, T 𝑇 T , of the 3-point correlation functions computed on each ensemble.

set	handle	$w_0/a$	$N_x^3\times N_t$	$n_{\text{cfg}}\times n_{\text{src}}$	$a{m}_{l0}$	$a{m}_{s0}$	$a{m}_{c0}$	$a{m}_{s0}^{\text{val}}$	$a{m}_{c0}^{\text{val}}$	$a{m}_{h0}^{\text{val}}$	T
1	fine	1.9006(20)	${32}^3\times 96$	$938\times 8$	0.0074	0.037	0.440	0.0376	0.45	0.5, 0.65, 0.8	14,17,20
2	fine-physical	1.9518(7)	${64}^3\times 96$	$284\times 4$	0.0012	0.0363	0.432	0.036	0.433	0.5, 0.8	14,17,20
3	superfine	2.896(6)	${48}^3\times 144$	$250\times 8$	0.0048	0.024	0.286	0.0234	0.274	0.427, 0.525, 0.65, 0.8	22,25,28
4	ultrafine	3.892(12)	${64}^3\times 192$	$249\times 4$	0.00316	0.0158	0.188	0.0165	0.194	0.5, 0.65, 0.8	31,36,41

Table 2. Table 2: Values extracted from correlation function fits for h A 1 s ​ ( 1 ) subscript superscript ℎ 𝑠 subscript 𝐴 1 1 h^{s}_{A_{1}}(1) , along with quantities required in our fits to determine a value at the physical point. Results are given on each gluon field ensemble for each valence heavy quark mass used. Results come from our simultaneous fits to two-point and three-point correlation functions: h A 1 s ​ ( 1 ) subscript superscript ℎ 𝑠 subscript 𝐴 1 1 h^{s}_{A_{1}}(1) values are determined using Eq. ( 19 ) and the ground-state meson masses in columns 4, 5, 6, 8, 9 and 10 from Eq. ( 12 ). f H c subscript 𝑓 subscript 𝐻 𝑐 f_{H_{c}} is the H c subscript 𝐻 𝑐 H_{c} meson decay constant determined from Eq. ( 35 ).

Set	$a{m}_{h0}^{\text{val}}$	$h_{A_1}^s\left(1\right)$	$a{M}_{H_s}$	$a{M}_{D_s^{\ast }}$	$a{M}_{H_c}$	$a{f}_{H_c}$	$a{M}_{{\eta }_h}$	$a{M}_{{\eta }_c}$	$a{M}_{{\eta }_s}$
1	0.5	0.9255(20)	0.95972(12)	0.96616(44)	1.419515(41)	0.186299(70)	1.471675(38)	1.367014(40)	0.313886(75)
	0.65	0.9321(22)	1.12511(16)		1.573302(40)	0.197220(77)	1.775155(34)
	0.8	0.9434(24)	1.28128(21)		1.721226(39)	0.207068(78)	2.064153(30)
2	0.5	0.9231(21)	0.95462(12)	0.93976(42)	1.400034(28)	0.183472(62)	1.470095(25)	1.329291(27)	0.304826(52)
	0.8	0.9402(27)	1.27577(22)		1.702456(23)	0.203407(45)	2.062957(19)
3	0.427	0.9107(46)	0.77453(24)	0.63589(49)	1.067224(46)	0.126564(70)	1.233585(41)	0.896806(48)	0.207073(96)
	0.525	0.9165(49)	0.88487(31)		1.172556(46)	0.130182(72)	1.439515(37)
	0.65	0.9246(65)	1.02008(39)		1.303144(46)	0.133684(75)	1.693895(33)
	0.8	0.9394(66)	1.17487(54)		1.454205(46)	0.137277(79)	1.987540(30)
4	0.5	0.9143(51)	0.80245(24)	0.47164(39)	1.011660(32)	0.098970(52)	1.342639(65)	0.666586(89)	0.15412(17)
	0.65	0.9273(62)	0.96386(33)		1.169761(34)	0.100531(60)	1.650180(56)
	0.8	0.9422(72)	1.11787(43)		1.321647(37)	0.101714(70)	1.945698(48)

Table 3. Table 3: Normalization constants applied to the lattice axial vector current in ( 19 ). Z A subscript 𝑍 𝐴 Z_{A} is found from ( 16 ) and Z disc subscript 𝑍 disc Z_{\text{disc}} from ( 17 ).

Set	$a{m}_{h0}^{\text{val}}$	$Z_A$	$Z_{\text{disc}}$
1	0.5	1.03178(57)	0.99819
	0.65	1.03740(58)	0.99635
	0.8	1.04368(56)	0.99305
2	0.5	1.03184(47)	0.99829
	0.8	1.04390(39)	0.99315
3	0.427	1.0141(12)	0.99931
	0.525	1.0172(12)	0.99859
	0.65	1.0214(12)	0.99697
	0.8	1.0275(12)	0.99367
4	0.5	1.00896(44)	0.99889
	0.65	1.01363(49)	0.99704
	0.8	1.01968(55)	0.99375

Table 4. Table 4: Error budget for h A 1 s ​ ( 1 ) subscript superscript ℎ 𝑠 subscript 𝐴 1 1 h^{s}_{A_{1}}(1) . Errors are given as a percentage of the final answer. The mass mistuning error includes that from valence strange and sea light and strange quarks; we find that taking a ± 10 ​ MeV plus-or-minus 10 MeV \pm 10\,\text{MeV} uncertainty in the physical value of the η b subscript 𝜂 𝑏 \eta_{b} mass has a negligible effect.

Source	% Fractional Error
Statistics $+Z_A$	1.06
$a\to 0$	0.73
$m_h\to m_b$	0.69
mass mistuning	0.20
Total	1.45

Equations

$$\frac{d\Gamma }{dw}$$

$$\times (M_{B_s}^2-M_{D_s^{\ast }}^2)\sqrt{w^2-1}\chi (w)|{\mathcal{F}}^{B_s\to D_s^{\ast }}(w){|}^2.$$

$$⟨D_s^{\ast }(ϵ)|V^{\mu }|B_s⟩$$

$$⟨D_s^{\ast }(ϵ)|A^{\mu }|B_s⟩$$

$$h_{A_2}^s(w)$$

$$⟨D_s^{\ast }(ϵ)|A^{\mu }|B_s⟩$$

$${\mathcal{F}}^{B_s\to D_s^{\ast }}(1)=h_{A_1}^s(1).$$

$$C_M(t)=$$

$${\Phi }_M(t)=\sum _{\mathbf{x}}\bar{q}(\mathbf{x},t)\Gamma q^{\prime }(\mathbf{x},t),$$

$$C_{H_s}(t)=\frac{1}{4}\sum _{\mathbf{x},\mathbf{y}}\text{Tr}\left[g_h(x,y)g_s^{†}(x,y)\right].$$

$$C_{D_s^{\ast }}(t)=\frac{1}{4}\sum _{\mathbf{x},\mathbf{y}}(-1{)}^{x_i+y_i}\text{Tr}\left[g_c(x,y)g_s^{†}(x,y)\right].$$

$$C_{{\hat{H}}_c}(t)=\frac{1}{4}\sum _{\mathbf{x},\mathbf{y}}(-1{)}^{{\bar{x}}_0+{\bar{y}}_0}\text{Tr}\left[g_h(x,y)g_c^{†}(x,y)\right],$$

$$C_{\text{3pt}}(t,T)=$$

$$A^{\mu }(\mathbf{y},t)=\bar{c}(\mathbf{y},t){\gamma }^5{\gamma }^{\mu }h(\mathbf{y},t).$$

$$C_{\text{3pt}}(t,T)=$$

$$\times \text{Tr}\left[g_h(x,y)g_c(y,z)g_s^{†}(z,x)\right],$$

$$C_M(t){|}_{\text{fit}}$$

$$\displaystyle-(-1)^{t}|a^{M,o}_{n}|^{2}f(E_{n}^{M,o},t)\Big{)};$$

$$f(E,t)$$

$$C_{\text{3pt}}$$

$$\displaystyle\Big{(}\,a^{H_{s}}_{j}J^{nn}_{jk}a^{D_{s}^{*}}_{k}f(E^{H_{s}},t)f(E^{D_{s}^{*}}_{n},T-t)$$

$$+a_j^{H_s,o}{J}_{jk}^{on}{a}_k^{D_s^{\ast }}(-1{)}^{t}f(E_n^{H_s,o},t)f(E^{D_s^{\ast }},T-t)$$

$$+a_j^{H_s}{J}_{jk}^{no}{a}_k^{D_s^{\ast },o}(-1{)}^{T-t}f(E^{H_s},t)f(E_n^{D_s^{\ast },o},T-t)$$

$$\displaystyle+a^{H_{s},o}_{j}J^{oo}_{jk}a^{D_{s}^{*},o}_{k}(-1)^{T}f(E^{H_{s},o}_{n},t)f(E^{D_{s}^{*},o},T-t)\,\Big{)}.$$

$$⟨0|{\Phi }_M|M⟩{|}_{\text{lat}}=\sqrt{2M_M}{a}_0^M.$$

$$⟨D_s^{\ast }(\hat{k})|A^k|H_s⟩{|}_{\text{lat}}=2\sqrt{M_{H_s}{M}_{D_s^{\ast }}}{J}_{00}^{nn}.$$

$$(m_{h0}^{\text{val}}+m_{c0}^{\text{val}})⟨0|P|H_c⟩{|}_{\text{lat}}=M_{{\hat{H}}_c}{Z}_A⟨0|A^0|{\hat{H}}_c⟩{|}_{\text{lat}}.$$

$$Z_{\text{disc}}$$

$${\tilde{C}}_q=\cosh a{m}_{q,\text{tree}}\left(1-\frac{1+{ϵ}_{q,\text{Naik}}}{2}{\sinh }^{2}a{m}_{q,\text{tree}}\right).$$

$$\displaystyle am_{q,\text{tree}}=am_{q0}\Big{(}1-{3\over 80}am_{q0}^{4}+{23\over 2240}am_{q0}^{6}$$

$$\displaystyle+{1783\over 537600}am_{q0}^{8}-{76943\over 23654400}am_{q0}^{10}+\mathcal{O}(am_{q0}^{12})\Big{)},$$

$$h_{A_1}^s(1)=\frac{1}{3}\sum _{k=1}^3\frac{Z_A{Z}_{\text{disc}}⟨D_s^{\ast }(\hat{k})|A^k|H_s⟩{|}_{\text{lat}}}{2\sqrt{M_{H_s}{M}_{D_s^{\ast }}}}.$$

$$h_{A_1}^s(1)(a,m_l,m_h)$$

$$h_{A_1}(1)$$

$$+\mathcal{O}\left(\frac{1}{m_c^n{m}_h^m},n+m\ge 3\right),$$

$${\eta }_A=1-\frac{{\alpha }_s}{\pi }\left(\frac{1+m_c/m_h}{1-m_c/m_h}\ln \frac{m_c}{m_h}+\frac{8}{3}\right).$$

$$m_h=M_{H_s}-{\bar{\Lambda }}_{\text{MRS}}-\frac{{\mu }_{\text{MRS}}^2}{M_{H_s}-{\bar{\Lambda }}_{\text{MRS}}}+\mathcal{O}\left(\frac{1}{m_h^2}\right).$$

$${\epsilon }_q\equiv \frac{2}{M_{{\eta }_q}}.$$

Full text

I Introduction

The study of quark flavour-changing interactions is a key component of the search for physics beyond the Standard Model (SM). There are currently a number of related tensions between experimental measurements and SM predictions Wei et al. (2009); Lees et al. (2012a, b, 2013); Aaij et al. (2014a, b); Huschle et al. (2015); Aaij et al. (2016a, 2015a, 2015b, b); Wehle et al. (2017); Sato et al. (2016); Hirose et al. (2018); Aaij et al. (2018a, b, 2017); Sirunyan et al. (2018); Aaboud et al. (2018), along with discrepancies between systematically independent determinations of Cabibbo-Kobayashi-Maskawa (CKM) matrix elements Amhis et al. (2017); Bevan et al. (2014); Alberti et al. (2015). A more precise understanding of these processes is needed to resolve these issues.

The $\bar{B}^{0}\to D^{*+}\ell^{-}\bar{\nu}$ decay (and its charge conjugate, that we simply abbreviate to $B\to D^{*}\ell\nu$ from now on) supplies one of the three methods used for precisely determining the CKM element $|V_{cb}|$ Schroder (1994); Bortoletto et al. (1990); Fulton et al. (1991); Albrecht et al. (1992); Barish et al. (1995); Buskulic et al. (1996, 1994); Abbiendi et al. (2000); Abreu et al. (2001); Adam et al. (2003); Abdallah et al. (2004); Aubert et al. (2008a, b, 2009); Dungel et al. (2010); Abdesselam et al. (2017); Bailey et al. (2014); Abdesselam et al. (2018). Measurements of branching fractions are extrapolated through $q^{2}$ to the zero recoil point to deduce $\mathcal{F}(1)|V_{cb}|$, where $\mathcal{F}(1)$ is the value of the only form factor contributing at zero recoil. Then a determination of $\mathcal{F}(1)$ in the Standard Model (via Lattice QCD Bailey et al. (2014); Harrison et al. (2018)) can be divided out to infer $|V_{cb}|$.

$|V_{cb}|$ is an important quantity and needs to be determined accurately. It constrains one side of the unitarity triangle via the ratio $|V_{ub}|/|V_{cb}|$. It is also a dominant uncertainty in the determination of the $CP$-violation parameter $\epsilon_{K}$ (where there is currently tension between the SM and experiment, see for example Bailey et al. (2018a)).

Previous determinations of $|V_{cb}|$ have shown systematic discrepancies with each other. The two competing values were those derived from exclusive decays ($B\to D^{*}\ell\nu$ and $B\to D\ell\nu$ with $B\rightarrow D^{*}$ giving the more accurare result), and inclusive ($B\to X_{c}\ell\nu$, where $X_{c}$ is any charmed hadronic state). In 2016 the Heavy Flavour Averaging Group (HFLAV) gave a value derived from exclusive $B\rightarrow D^{*}$ decays of $|V_{cb}|_{\text{excl}}=(39.05\pm 0.47_{\text{exp}}\pm 0.58_{\text{th}})\times 10^{-3}$ and from inclusive decays, using the kinetic scheme, of $|V_{cb}|_{\text{incl}}=(42.19\pm 0.78)\times 10^{-3}$ Amhis et al. (2017). It has since been suggested, based on unfolded Belle data Abdesselam et al. (2017), that the tension seen here arose (at least partly) from the use of a very constrained parameterization in the extrapolation of the experimental $B\rightarrow D^{*}$ decay rates to zero recoil Bernlochner et al. (2017); Bigi et al. (2017); Grinstein and Kobach (2017). Recent exclusive determinations of $V_{cb}$ have then used a less constrained parameterisation to give a larger, and less precise, result for $V_{cb}$ that is no longer in tension with the inclusive result. For example, the Particle Data Group quote $|V_{cb}|_{\text{excl}}=(41.9\pm 2.0)\times 10^{-3}$ Tanabashi et al. (2018). However, an even more recent $V_{cb}$ determination from $B\to D^{*}\ell\nu$ data by the BaBar collaboration Lees et al. (2019) used the less constrained parameterisation but still found a tension with the inclusive result. This clearly points to the need for more work to improve the accuracy of the exclusive result. On the theory side a better understanding of the form factors for $B\rightarrow D^{*}$ from lattice QCD is required, both at zero recoil and away from zero recoil.

Another motivation for studying $B\to D^{*}\ell\nu$ is the tension between SM and experimental determinations of the ratio $R_{D^{(*)}}=\mathcal{B}(\bar{B}\to D^{(*)}\tau\bar{\nu}_{\tau})/\mathcal{B}(\bar{B}\to D^{(*)}\ell\bar{\nu}_{\ell})$ ($\ell=e$ or $\mu$). The latest HFLAV report gives the combined statistical significance of the anomalies in $R_{D}$ and $R_{D^{*}}$ to be $3.8\sigma$ Amhis et al. (2017). A preliminary new analysis from Belle Caria (2019), however, gives results closer to the SM and pulls the global average down to $3.1\sigma$. More precise measurements and predictions will either confirm or dismiss a new physics explanation.

The weak decay process $B_{s}\to D^{*}_{s}\ell\nu$ is very similar to $B\to D^{*}\ell\nu$ and could also be used to determine $|V_{cb}|$ and test the SM. It is feasible to study this decay at LHC and from the theoretical side it is a more attractive channel than $B\to D^{*}$. The absence of valence light quarks means that lattice QCD results have smaller statistical errors and are less computationally expensive. Finite-volume effects and the dependence on $u/d$ quark masses (for quarks in the sea) are also smaller. The $D_{s}^{*}$ has no Zweig-allowed strong decay mode, unlike the $D^{*}$, and is in fact a relatively long-lived particle Donald et al. (2014a) that can be considered ‘gold-plated’ in lattice QCD. This makes the $B_{s}\to D^{*}_{s}\ell\nu$ both a useful test bed for lattice techniques (that may be later used to study $B\to D^{*}\ell\nu$ decays) and a key decay process for which to make predictions ahead of experimental results.

Lattice QCD calculations have shown that several weak decay form factors are relatively insensitive to whether the spectator quark is a $u/d$ or $s$ quark Bailey et al. (2012); Koponen et al. (2013); Monahan et al. (2017). A combination of chiral perturbation theory and Heavy Quark Symmetry Jenkins and Savage (1992) backs up this expectation for $B$ decays. We can therefore expect the form factors to be very similar for $B_{s}\to D^{*}_{s}$ and $B\to D^{*}$. A recent lattice calculation Harrison et al. (2018) found an insignificant $\mathcal{O}(1\%)$ difference at zero recoil: $\mathcal{F}^{B\to D^{*}}(1)/\mathcal{F}^{B_{s}\to D_{s}^{*}}(1)=1.013(14)_{\text{stat}}(17)_{\text{sys}}$. Information from the study of $B_{s}\to D_{s}^{*}$ can then be applicable to $B\to D^{*}$.

Lattice QCD calculations of the $B_{(s)}\to D_{(s)}^{*}$ form factors at zero recoil have so far been performed by two collaborations using different methods. The Fermilab Lattice and MILC collaborations calculated $\mathcal{F}^{B\to D^{*}}(1)$ in Bernard et al. (2009); Bailey et al. (2014) using the Fermilab action for both $b$ and $c$ quarks El-Khadra et al. (1997) and asqtad $u/d$ quarks Bazavov et al. (2010a). More recently the HPQCD collaboration computed both $\mathcal{F}^{B\to D^{*}}(1)$ and $\mathcal{F}^{B_{s}\to D_{s}^{*}}(1)$ Harrison et al. (2018) using improved NRQCD $b$ quarks Lepage et al. (1992); Dowdall et al. (2012) and Highly Improved Staggered (HISQ) $c$ and $u/d/s$ quarks Follana et al. (2007). The RBC/UKQCD Flynn et al. (2016) and LANL-SWME Bailey et al. (2018b) collaborations are also working towards these form factors using variants of the Fermilab action for heavy quarks and JLQCD has a calculation in progress using Möbius domain-wall quarks Kaneko et al. (2018).

The formalism to use for the heavy quarks is a major consideration in designing a lattice QCD calculation to determine these form factors. Most of the calculations discussed in the previous paragraph (apart from the JLQCD calculation) use approaches that make use of the nonrelativistic nature of heavy quark bound states to tune the $b$ (and in some cases also $c$) quark masses. This avoids potentially large discretisation effects appearing in the results in the form of a systematic error of size $(am_{b})^{n}$, where $n$ is an integer that depends on the level of improvement in the action. The absence of these discretisation errors means that $b$ quarks can be handled on relatively coarse lattices where $am_{b}>1$. However the price to be paid is that the current operator that couples to the $W$ boson is also implemented within a nonrelativistic framework and must then be renormalised to match the appropriate operator in continuum QCD. This matching can be done using perturbative lattice QCD but has only be done through $\mathcal{O}(\alpha_{s})$ for these actions Harada et al. (2002); Monahan et al. (2013). This leaves a substantial source of uncertainty from missing higher-order terms in the perturbative matching that is not easily reduced. This matching uncertainty contributes $\sim$ 80% of the final error in the HPQCD calculation Harrison et al. (2018) and $\sim$ 30% in the Fermilab/MILC calculation Bailey et al. (2014) because of the differing allowances for missing higher-order terms.

Here we report details and results of a calculation of the $B_{s}\to D^{*}_{s}$ form factor at zero recoil using an approach free of perturbative matching uncertainties. We perform our calculation on the second-generation MILC ensembles Bazavov et al. (2010b, 2013), including effects from 2+1+1 flavours in the sea using the HISQ action. We also use the HISQ action for all valence quarks. We obtain results at a number of differing masses for the $b$ (we refer to this generically as the heavy quark $h$), and perform an extrapolation to $m_{h}=m_{b}$. By using only HISQ quarks, we can obtain the normalizations of all required currents fully non-perturbatively. We refer to this as the heavy-HISQ approach. By using many heavy masses and multiple values of the lattice spacing, including very fine lattices, we can model both the form factor dependence on the heavy mass, and the discretisation effects associated with using large $am_{h}$ values.

The heavy-HISQ approach was developed by HPQCD to compute $B$ meson masses and decay constants McNeile et al. (2012a, b) and the $b$ quark mass McNeile et al. (2010); Chakraborty et al. (2015). It is also now being used by other collaborations for these calculations Bazavov et al. (2018a); Petreczky and Weber (2019). A proof-of-principle application of heavy-HISQ to form factors was given in Lytle et al. (2016); Colquhoun et al. (2016) for $B_{c}$ decays, showing that the full $q^{2}$ range of the decay could be covered. Here we extend the approach to form factors for $B_{s}$ decays but working only at zero recoil, a straightforward extension of earlier work. Using the heavy-HISQ approach also has the added benefit of eludicating the dependence of form factors on heavy quark masses, meaning that we can test expectations from Heavy Quark Effective Theory (HQET).

This article is structured in the following way: Section II defines the form factor and gives details of the lattice calculation, including the nonperturbative normalisation and extrapolation in heavy-quark mass; Section III presents our results and compares to earlier calculations and Section IV gives our conclusions and outlook. In the appendix, we give details of a number of tests we performed on the correlator fits and the continuum, chiral and heavy-quark extrapolations.

II Calculation Details

II.1 Form Factors

The differential decay rate for the $\bar{B}_{s}^{0}\to D_{s}^{*+}l^{-}\bar{\nu}_{l}$ decay is given in the SM by

[TABLE]

where $w=v_{B_{s}}\cdot v_{D^{*}_{s}}$, $v=p/M$ is the 4-velocity of each meson, and $\chi(w)$ is a known function of $w$ with $\chi(1)=1$ (see, for example, appendix G of Harrison et al. (2018)). $\bar{\eta}_{\text{EW}}$ accounts for electroweak corrections from diagrams where photons or $Z$s are exchanged in addition to a $W^{-}$, as well as the Coulomb attraction of the final-state charged particles Sirlin (1982); Ginsberg (1968); Atwood and Marciano (1990). The differential decay rate for the $B_{s}^{0}\to D_{s}^{*-}l^{+}\bar{\nu}_{l}$ is identical.

The form factor $\mathcal{F}^{B_{s}\to D_{s}^{*}}(w)$ is a linear combination of hadronic form factors that parameterize the vector and axial-vector matrix elements between initial and final state hadrons. A common choice of parameterization used in the context of Heavy Quark Effective Theory (HQET) is Richman and Burchat (1995)

[TABLE]

where $V^{\mu}=\bar{c}\gamma^{\mu}b$ is the vector $b\to c$ current and $A^{\mu}=\bar{c}\gamma^{\mu}\gamma^{5}b$ is the axial-vector current. $\epsilon$ is the polarization 4-vector of the $D_{s}^{*}$ final state.

At zero recoil ($w=1$), the vector matrix element vanishes, the axial-vector element simplifies to

[TABLE]

and $\mathcal{F}^{B_{s}\to D_{s}^{*}}(w)$ reduces to

[TABLE]

Our goal is to compute $h^{s}_{A_{1}}(1)$.

All we need to do this is the matrix element $\langle D^{*}_{s}(\epsilon)|A^{\mu}|B_{s}\rangle$ with both the $B_{s}$ and $D_{s}^{*}$ at rest, with the $D_{s}^{*}$ polarization $\epsilon$ in the same direction as the (spatial) axial-vector current.

II.2 Lattice Calculation

The gluon field configurations that we use were generated by the MILC collaboration Bazavov et al. (2010b, 2013). Table 1 gives the relevant parameters for the specific ensembles that we use. The gluon field is generated using a Symanzik-improved gluon action with coefficients calculated through $\mathcal{O}(\alpha_{s}a^{2},n_{f}\alpha_{s}a^{2})$ Hart et al. (2009). The configurations include the effect of 2+1+1 flavours of dynamical quarks in the sea ($u$,$d$,$s$,$c$, with $m_{u}=m_{d}\equiv m_{l}$), using the HISQ action Follana et al. (2007). In three of the four ensembles (fine, superfine and ultrafine), the bare light quark mass is set to $m_{l0}/m_{s0}=0.2$. The fact that the $m_{l0}$ value is unphysically high is expected to have only a small effect on $h^{s}_{A_{1}}(1)$, because there are no valence light quarks. The effect is quantified here by including a fourth ensemble (fine-physical) with (approximately) physical $m_{l0}$.

We use a number of different masses for the valence heavy quark, $h$. This is in order to resolve the dependence of $h_{A_{1}}^{s}(1)$ on the heavy mass, so that an extrapolation to $m_{h}=m_{b}$ can be performed. By varying the heavy mass on each ensemble and by using ensembles at varying small lattice spacing, we can resolve both the discretisation effects that grow with heavy quark mass ($am^{\text{val}}_{h0}\lesssim 1$) and the physical dependence of the continuum form factor on $m_{h}$.

Staggered quarks have no spin degrees of freedom, so that solution of the Dirac equation on each gluon field is numerically fast. The remnant of the doubling problem means that quark bilinears of specific spin-parity have multiple copies, called ‘tastes’ Follana et al. (2007). They differ in the amount of point-splitting between the fields and the space-time dependent phase needed to substitute for the appropriate $\gamma$ matrix. In this calculation we can use only local (non point-split) bilinears, which is an advantage in terms of statistical noise, since no gluon fields are included in the current operator. In the standard staggered spin-taste notation, the operators that we use are: pseudoscalar, $\Gamma_{P}=(\gamma^{5}\otimes\gamma^{5})$; vector, $\Gamma_{V}^{\mu}=(\gamma^{\mu}\otimes\gamma^{\mu})$ and axial-vector, $\Gamma_{A}^{\mu}=(\gamma^{\mu}\gamma^{5}\otimes\gamma^{\mu}\gamma^{5})$.

We compute several ‘two-point’ correlation functions on the ensembles detailed in table 1, combining HISQ propagators from solving the Dirac equation for each random wall time source. These correlation functions take the form

[TABLE]

where $\langle\rangle$ represents a functional integral, $q,q^{\prime}$ are valence quark fields of the flavours the $M$ meson is charged under, $\Gamma$ is the spin-taste structure of $M$ and $1/N_{\text{taste}}$ is the staggered quark normalisation for closed loops. The random-wall source and the sum over $\bf{x}$ at the sink project onto zero spatial momentum. We compute the correlation functions for all $t$ values, i.e. $0\leq t\leq N_{t}$.

The correlation function for the heavy-strange pseudoscalar meson, $H_{s}$, with valence quark content $h\overline{s}$ and spin-taste structure $\Gamma_{P}$ is constructed from HISQ propagators as:

[TABLE]

Here $g_{q}(x,y)$ is a HISQ propagator for flavour $q$, the trace is over color and $1/4$ is the staggered quark normalisation. $x_{0}=0$ and $y_{0}=t$ and the sum is over spatial sites $\bf{x}$, $\bf{y}$. We also compute correlators for a charm-strange vector meson $D_{s}^{*}$, with structure $\Gamma^{i}_{V}$, using

[TABLE]

We average over polarisations, $i=1,2,3$.

We also compute correlation functions for two tastes of pseudoscalar heavy-charm mesons denoted $H_{c}$ and $\hat{H}_{c}$ respectively. $H_{c}$ has spin-taste structure $\Gamma_{P}$, and $\hat{H}_{c}$ has structure $\Gamma^{0}_{A}$. $H_{c}$ correlators are computed using Eq. (7) (with $g_{s}$ replaced with $g_{c}$), while $\hat{H}_{c}$ correlators are given by

[TABLE]

where we use the notation $\bar{z}_{\mu}=\sum_{\nu\neq\mu}z_{\nu}$. These correlators will be used to normalise the axial vector $b\overline{c}$ current as discussed in Section II.4.

A useful physical proxy (that does not run) for the quark mass is that of the pseudoscalar meson made from that flavour of quark. It is therefore also useful, for our heavy quark mass extrapolation, to calculate correlation functions for heavy-heavy pseudoscalars, denoted $\eta_{h}$, with spin-taste structure $\Gamma_{P}$ using Eq. (7). Likewise, to test the impact of any mistuning of the charm and strange quark masses, we also determine $\eta_{c}$ and $\eta_{s}$ correlators analogously. We can tune the $c$ and $b$ masses using the experimental values for the $\eta_{c}$ and $\eta_{b}$ masses, allowing for slight shifts from missing QED effects and the fact that we do not allow these mesons to annihilate to gluons McNeile et al. (2010). The mass of the $\eta_{s}$ meson (which is not a physical state) can be fixed in lattice QCD from the $K$ and $\pi$ meson masses Davies et al. (2010a); Dowdall et al. (2013).

We then generate the ‘three-point’ correlation functions that contain the $H_{s}$ to $D_{s}^{*}$ transition.

[TABLE]

Our $H_{s}$ source is given spin-taste $\Gamma_{P}$, the $D_{s}^{*}$ sink, $\Gamma^{i}_{V}$, and the current insertion $\Gamma^{i}_{A}$. This gives the required cancellation of tastes within the three-point function Donald et al. (2014b). In terms of HISQ propagators

[TABLE]

where we fix $x_{0}=0$, $y_{0}=t$ and $z_{0}=T$. We compute the three-point correlation functions for all $t$ values within $0\leq t\leq T$, and 3 $T$ values that vary between ensembles and are given in Table 1. We average over the 3 directions for $i$ for increased statistical precision.

II.3 Analysis of Correlation Functions

We use simultaneous Bayesian fits Lepage et al. (2002); cor (2018) to extract the axial vector matrix element and meson masses from the two- and three-point correlation functions. This allows us to include the covariance between results at different heavy quark masses on a given ensemble into our subsequent fits in Section II.5.

We fit the two-point correlation functions using the functional form

[TABLE]

where $N_{t}$ is the temporal extent of the lattice, and $E^{M,(o)}_{n}$,$a^{M,(o)}_{n}$ are fit parameters, with the excited-state energy parameters implemented as energy differences to the state below Lepage et al. (2002). The second term accounts for the presence of opposite-parity states that contribute an oscillating term to the correlation function when using staggered quarks Follana et al. (2007). These terms do not appear when $M$ is a pseudoscalar with a quark and antiquark of the same mass, so in the $M=\eta_{h},\eta_{c},$ and $\eta_{s}$ cases the second term is not required. For all correlator fits we set $N_{\text{exp}}=5$; this allows the impact of systematic effects from excited states to be included in the ground-state parameters that we are interested in.

The three-point correlation functions have the fit form

[TABLE]

This includes fit parameters common to the fits of the $H_{s}$ and $D_{s}^{*}$ two-point correlators, along with new fit parameters $J_{jk}$.

We perform a single simultaneous fit containing each correlator computed ($H_{s},D_{s}^{*},\eta_{h},\eta_{c},\eta_{s},H_{c},\hat{H}_{c}$, and three-point) for each ensemble. We set gaussian priors for the parameters $J_{jk}$, and log-normal priors for all other parameters. Using log-normal distributions forbids energy differences $E_{n+1}^{M}-E_{n}^{M}$ and amplitudes $a_{n}^{M}$ (which can be taken to be positive here) from moving too close to zero or changing sign, improving stability of the fit.

Ground state energies $E_{0}^{M}$ are given priors of $(am_{q0}+am_{q^{\prime}0}+a\Lambda_{\text{QCD}})\pm 2a\Lambda_{\text{QCD}}$, where $m_{q0}$ and $m_{q^{\prime}0}$ are the masses of the appropriate quarks, and $\Lambda_{\text{QCD}}$ is the confinement scale, which we set to 0.5GeV. For $q=h$ or $c$, this corresponds to the leading order HQET expression for a heavy meson mass. Ground-state energies of oscillating states, $E_{0}^{M,o}$, are given priors of $(am_{q0}+am_{q^{\prime}0}+2a\Lambda_{\text{QCD}})\pm 2a\Lambda_{\text{QCD}}$. Excited state energy differences, $E_{i+1}^{M}-E_{i}^{M}$, $i>0$ are given prior values $2a\Lambda_{\text{QCD}}\pm a\Lambda_{\text{QCD}}$. Priors for ground state amplitudes $a_{0}^{M}$, are set from plots of effective amplitudes. The resulting priors always have a variance at least 10 times that of the final result for the ground-state. We use log(amplitude) priors of -1.20(67) for non-oscillating excited states and -3.0(2.0) for oscillating excited states. The ground-state non-oscillating to non-oscillating 3-point parameter, $J_{00}^{nn}$ is given a prior of $1\pm 0.6$, and the rest of the 3-point parameters $J_{jk}^{nn}$ are given $0\pm 1$.

$E_{0}^{M}=aM_{M}$ is the mass of the ground-state meson $M$ in lattice units. The masses $M_{H_{s}}$ and $M_{\eta_{h}}$ can both be used as proxys for $m_{h}$ in the extrapolation to $m_{h}=m_{b}$. The annihilation amplitude for an $M$-meson is given (in lattice units) by

[TABLE]

The (as yet unnormalised) matrix element that we need to obtain $h_{A_{1}}^{s}(1)$ is given by

[TABLE]

To ensure that truncating the sum over states at $N_{\text{exp}}$ accounts for the full systematic error from excited states, we cut out some data very close to the sources and sinks, where even higher excited states might have some effect. To do this we only include data with $t\geq t_{\text{cut}}$ and $t\leq N_{t}-t_{\text{cut}}$ in the two-point case and $t\leq T-t_{\text{cut}}$ in the three-point case. We can in principle use a different $t_{\text{cut}}$ for every correlation function included in our fit, but we do not use a big range of $t_{\text{cut}}$ values. They range from 1 to 3 for the three-point functions and up to 8 for the two-point functions.

The determination and minimisation of the $\chi^{2}$ function in our fit procedure requires the inversion of the covariance matrix that captures correlations between the different pieces of ‘data’ (correlation functions) in our fit. The low eigenmodes of the correlation matrix are not well determined with the statistics that we have and so we implement an SVD (singular value decomposition) cut in the inversion of the correlation matrix to avoid underestimating the uncertainty in the parameters of the fit cor (2018). This replaces correlation matrix eigenvalues below $\lambda_{\text{min}}$, equal to svdcut times the largest eigenvalue, with $\lambda_{\text{min}}$. $\lambda_{\text{min}}$ is estimated using the diagnosis tools in the Corrfitter package cor (2018) and corresponds typically to an svdcut of $10^{-3}$ here.

Figure 1 summarises stability tests of our fits, focussing on the key parameter $J_{00}^{nn}$ that is converted to the ground-state to ground-state transition amplitude using Eq. (15).

The fit parameters determined by our fits that we will use to calculate the physical value for $h^{s}_{A_{1}}(1)$ are given in Table 2. Notice that the statistical errors on the results grow with the heavy quark mass. This is a well understood problem in lattice heavy-light meson physics (see, for example Davies et al. (2010b)). Our method here has the advantage of including information from lighter-than-$b$ heavy quarks with improved statistical precision.

II.4 Normalisation of the Axial Current

The partially-conserved axial-vector current for the HISQ action is a complicated linear combination of one-link and three-link lattice currents. In this study we use only local axial vector currents. This simplifies the lattice QCD calculation but creates the need for our resulting current matrix element to be multiplied by a matching factor $Z_{A}$ to produce the appropriate continuum current. We determine $Z_{A}$ via a fully non-perturbative method Donald et al. (2014b).

We use the fact that the staggered local pseudoscalar current of spin-taste $(\gamma^{5}\otimes\gamma^{5})$, multiplied by the sum of its valence quark masses, is absolutely normalized via the PCAC relation. From the two-point $H_{c}$ and $\hat{H}_{c}$ correlator fits we can extract the decay amplitudes: $\langle 0|\bar{c}(\gamma^{5}\otimes\gamma^{5})h|H_{c}\rangle\equiv\langle 0|P|H_{c}\rangle$ and $\langle 0|\bar{c}(\gamma^{0}\gamma^{5}\otimes\gamma^{0}\gamma^{5})h|\hat{H}_{c}\rangle=\langle 0|A^{0}|\hat{H}_{c}\rangle$ as in Eq. (14). Then, the normalization for the local $A^{0}$ current (common to that of the local spatial axial-vector current $A^{k}$ up to discretisation effects), $Z_{A}$, is fixed by demanding that

[TABLE]

The $Z_{A}$ values found on each ensemble and $am^{\text{val}}_{h0}$ are given in Table 3.

There is an ambiguity in what mass to use on the right hand side of Eq. (16). We use the non-goldstone mass $M_{\hat{H}_{c}}$, but one could just as well replace this with $M_{H_{c}}$ since the difference is a discretisation effect. The meson mass difference is very small for heavy mesons Follana et al. (2007) and so we find the effect of changing the taste of meson mass used never exceeds 0.15% of $Z_{A}$ throughout the range of ensembles and heavy masses that we use and has no impact on the continuum result.

We also remove tree-level mass-dependent discretisation effects coming from the wavefunction renormalisation Follana et al. (2007) by multiplying by a factor $Z_{\text{disc}}$. This is derived in Monahan et al. (2013) as:

[TABLE]

See also Bazavov et al. (2018a). $m_{q,\text{tree}}$ is the tree-level pole mass in the HISQ action. It has an expansion in terms of the bare mass Follana et al. (2007)

[TABLE]

$\epsilon_{q,\text{Naik}}$ fixes the Naik parameter Naik (1989) ($N=1+\epsilon$ is the coefficient of the tree-level improvement term for the derivative) in the HISQ action when it is being used for heavy quarks Follana et al. (2007). $\epsilon_{q,\text{Naik}}$ is set to its tree-level value, removing the leading tree-level errors from the dispersion relation. As an expansion in $am_{q,\text{tree}}$ it begins at $\mathcal{O}(am_{q,\text{tree}})^{2}$ Follana et al. (2007). To determine $\epsilon_{q,\text{Naik}}$ we use the closed form expression for it given in Monahan et al. (2013) and this can also be used along with Eq. (18) to evaluate $Z_{\text{disc}}$. The pole condition can be used to show that the expansion of $\tilde{C}_{q}$ begins at $am^{4}_{q0}$ as $1-3am^{4}_{q0}/80+\ldots$. The effect of $Z_{\text{disc}}$ is then very small, never exceeding $0.2\%$. $Z_{\text{disc}}$ values on each ensemble for each $am^{\text{val}}_{h0}$ are given in table 3.

Combining these normalizations with the lattice current from the 3-point fits, we find a value for the form factor at a given heavy mass and lattice spacing:

[TABLE]

II.5 Obtaining a Result at the Physical Point

We now discuss how we fit our results for the zero recoil form factor, $h_{A_{1}}^{s}(1)$, as a function of valence heavy quark mass, sea light quark mass and lattice spacing to obtain a result at the physical point where the heavy quark mass is that of the $b$, the sea quark masses are physical and the lattice spacing is zero.

In summary, we fit our results for $h_{A_{1}}^{s}(1)$ to the following form

[TABLE]

The terms in the first line allow for dependence on the valence heavy quark and charm quark masses (with $\varepsilon_{q}\equiv 1/m_{q}$) using input from HQET, to be discussed below. $\mathcal{N}_{\text{disc}}$ and $\mathcal{N}_{\text{mistuning}}$ account for discretisation and mass mistuning effects, also discussed below. The physical result is then $h^{s}_{A_{1}}(1)(0,m_{l,\text{phys}},m_{b})$.

II.5.1 Dependence on the heavy valence quark mass

Our fit of the $m_{h}$ dependence is guided by HQET, which considers both the $c$ quark and the heavy quark of mass $m_{h}$ to be heavy here. In particular, for the parameter $h^{(s)}_{A_{1}}(1)$, HQET forbids terms of $\mathcal{O}(1/m_{Q})$ where $m_{Q}$ can be $m_{c}$ or $m_{b}$ Luke (1990). The HQET expression for $h_{A_{1}}(1)$ is then given by Falk and Neubert (1993); Mannel (1994):

[TABLE]

where $l_{V}$, $l_{A}$ and $l_{P}$ are $\mathcal{O}(\Lambda^{2}_{\text{QCD}})$ (with possible mild dependence on whether the spectator quark is $s$ or $u/d$). $\eta_{A}$ accounts for ultraviolet matching between HQET and QCD, and has been computed to 2-loops in perturbative QCD Czarnecki (1996). It has mild dependence on $m_{h}$ through logarithms of $m_{c}/m_{h}$; at one-loop $\eta_{A}$ has explicit form Close et al. (1984)

[TABLE]

The coefficient of $\alpha_{s}/\pi$ then varies from -0.66 to -0.29 across the range of $m_{h}$ from $m_{h}=m_{c}$ to $m_{h}=m_{b}$, taking $m_{b}/m_{c}=4.577(8)$ Bazavov et al. (2018b). The two-loop correction is small Czarnecki (1996). $\eta_{A}$ is then close to 1 and differs by a few percent across our range of $m_{h}$.

Our calculation has results at multiple values of $m_{h}$, and could therefore in principle provide information on the coefficients $l_{A}$ and $l_{P}$ of the $m_{h}$-dependent terms in the HQET expansion. The charm quark mass is fixed to its physical value and so we cannot access the value of $l_{V}$ independent of a choice of $\eta_{A}$ at $m_{h}=m_{c}$. The terms in round brackets in Eq. (21), multiplying $\eta_{A}$, are all very small because of the suppression by heavy-quark masses. To constrain them tightly requires very precise data and, as we will see, we are not able to determine $l_{A}$, $l_{P}$ or $l_{V}$ accurately with our results. It therefore does not make sense to attempt to compare them accurately to HQET expectations. To do so would require using an appropriate quark mass definition (since different definitions will move quark mass dependence between the $l_{A}$ term and the others in Eq. (21) ) and the two-loop expression for $\eta_{A}$ with appropriate value for $\alpha_{s}$ (since logarithmic $m_{h}$ dependence of $\eta_{A}$ can be misinterpreted as part of a polynomial in $1/m_{h}$).

Instead we simply take an HQET-inspired form for the $m_{h}$-dependence and set $\eta_{A}$ to 1, resulting in the first line of our fit form, Eq. (20). This is sufficient to test, through the results we obtain for $l_{A}$, $l_{V}$ and $l_{P}$ using this expression, that the HQET expectation for the approximate size of these coefficients is fulfilled. We take priors on $l_{A,V,P}$ in our fit of $0\pm 1\,\text{GeV}^{2}$.

We have several different proxies, derived from heavy meson masses, that we can take for the heavy quark mass that appears in $\varepsilon_{h}$ in Eq. (20). We do not expect our physical result for $h^{s}_{A_{1}}$ to vary significantly depending on which meson mass we use, but the results for $l_{A}$, $l_{V}$ and $l_{P}$ will vary because of different sub-leading terms in the relationship between meson and quark mass. The most obvious substitutions to use for the heavy quark mass are the mass of the pseudoscalar heavy-strange meson, $M_{H_{s}}$, and half the mass of the pseudoscalar heavyonium meson, $M_{\eta_{h}}$. We also tested using the quark mass in the minimal renormalon subtracted (MRS) scheme suggested in Brambilla et al. (2018). This takes

[TABLE]

where $\mu_{\text{MRS}}^{2}=\mu^{2}_{\pi,\text{MRS}}-d_{H^{(*)}}\mu^{2}_{G,\text{MRS}}$ with $d_{H^{(*)}}=1$ for pseudoscalar mesons and $-1/3$ for vectors. For this case we use parameters determined in Bazavov et al. (2018b) for the MRS scheme: $\bar{\Lambda}_{\text{MRS}}=0.552(30)\text{GeV}$, $\mu^{2}_{\pi,\text{MRS}}=0.06(22)\text{GeV}^{2}$ and $\mu^{2}_{G,\text{MRS}}=0.38(1)\text{GeV}^{2}$. We take $m_{h}$ from Eq. (23) using our results for the mass of the pseudoscalar heavy-strange meson and $m_{c}$ from our results for the mass of the $D_{s}^{*}$ meson.

We take our central fit, for simplicity, from the result of using half the pseudoscalar heavyonium mass for $m_{h}$ and half the pseudoscalar charmonium mass for $m_{c}$ i.e. taking

[TABLE]

We test the stability of the fit results under the different choices discussed above in Section III.2.

II.5.2 Mistuning of other quark masses

Our calculation has results for multiple different heavy quark masses on each gluon field configuration. The valence charm and strange quark masses, however, are tuned to their physical values. This is done by fixing the $\eta_{c}$ and $\eta_{s}$ meson masses to their physical values in a pure QCD world allowing, for example, for $\eta_{c}$ annihilation as discussed in Chakraborty et al. (2015). Any possible mistuning of the charm quark mass is accounted for in our fit function by the dependence on the charm quark mass that is included in the first line of Eq. (20). When the fit function is evaluated at the physical point we set $\varepsilon_{c}$ from the physical $\eta_{c}$ mass.

The strange (valence and sea) and light (sea) mass mistunings are accounted for using the tuning in Chakraborty et al. (2015). For the strange quark, we define $\delta_{s}=m_{s}-m_{s}^{\text{tuned}}$, where $m_{s}^{\text{tuned}}$ is given by

[TABLE]

$M_{\eta_{s}}^{\text{physical}}$ is determined in lattice simulations from the masses of the pion and kaon Dowdall et al. (2013). The ratio $\delta_{s}/m_{s}^{\text{tuned}}$ then gives the fractional mistuning. The valence strange quark masses are very well tuned, but the sea strange quark masses less so.

We similarly account for mistuning of the masses of the (sea) light quarks by defining $\delta_{l}=m_{l0}-m_{l}^{\text{tuned}}$. We find $m_{l}^{\text{tuned}}$ from $m_{s}^{\text{tuned}}$, leveraging the fact that the ratio of quark masses is regularization independent, and the ratio was calculated in Bazavov et al. (2018a):

[TABLE]

We set $m_{l}^{\text{tuned}}$ to $m_{s}^{\text{tuned}}$ divided by this ratio.

The full term we include to account for mistuning is then given by

[TABLE]

where $c_{l}$, $c_{s}$ and $c_{s}^{\text{val}}$ are fit parameters with prior distributions $0\pm 1$. We neglect $\delta_{s,l}^{2}$ contributions since these are an order of magnitude smaller and are not resolved in the results of our lattice calculation.

The gluon field configurations that we use have $m_{u}=m_{d}\equiv m_{l}$ in the sea. In the real world this is not true. We test the impact of possible isospin-breaking on our fits by testing for sensitivity to the sea light quark masses. We do this by changing the $m_{l}^{\text{tuned}}$ value up and down by the expected value for $m_{d}-m_{u}$ Tanabashi et al. (2018). We find the effect to be completely negligible in comparison to the other sources of error.

II.5.3 Discretisation Effects

Discretisation effects in our lattice QCD results are accounted for following the methodology of McNeile et al. (2012b). We take

[TABLE]

The leading terms, with $i=0$, allow for discretisation effects that are set by the heavy quark mass and also discretisation effects that are set by the charm quark mass (or indeed any other lighter scale that is independent of heavy quark mass). The $i>0$ terms allow for discretisation effects to vary as the heavy quark mass is varied, with $M_{\eta_{h}}$ being used here as a proxy for the heavy quark mass. $d_{ijk}$ are fit parameters with prior distributions $0\pm 1.0$. All discretisation effects are of even order by construction of the HISQ action Follana et al. (2007).

We tested the impact on the fit of including extra discretisation effects set by the scale $\Lambda_{\text{QCD}}$ but this made no difference (since such effects are much smaller than those already included by the $am^{\text{val}}_{c0}$ terms). We also tested the effects of increasing the number of terms in each sum, but the final result remained unchanged.

II.5.4 Finite-volume Effects

The finite volume effects in our lattice results are expected to be negligible, because we are working with heavy mesons that have no valence light quarks and no Zweig-allowed strong decay modes. Coupling to chiral loops or decay channels with pions that could produce significant finite-volume effects Laiho and Van de Water (2006) is therefore absent and we can safely ignore finite-volume effects here.

In section III.2 we detail the results of several tests of the stability of our final result under changes to the details of the fit.

II.5.5 Topological Charge Effects

It has been observed that the finest MILC ensembles ($a\simeq 0.45$fm and finer) suffer from slow variation in the topological charge Bernard and Toussaint (2018) with Monte Carlo time. The question then arises if physical observables obtained by averaging over the ensemble could be biassed by being measured in only a small range of topological charge sectors. This issue is addressed in Bernard and Toussaint (2018) through a calculation of the ‘topological adjustment’ needed for meson masses and decay constants on the ultrafine lattices used here (with $m_{l}/m_{s}=0.2$). The adjustment found for the $D_{s}$ decay constant is 0.002%. We might expect the impact of a frozen topological charge on $h_{A_{1}}^{s}(1)$ to be of a similar size to this, given that it involves a transition between heavy-strange mesons. Allowing for this systematic uncertainty (or even ten times it) has a negligible effect on our final result.

III Results and Discussion

III.1 Result for $h^{s}_{A_{1}}(1)$

The results of our correlation function fits (discussed in Section II.3) are given in Table 2. We tabulate values for $h^{s}_{A_{1}}(1)$ at each heavy quark mass that we have used on each gluon field ensemble from Table 1. We also tabulate the meson masses needed to allow determination of $h^{s}_{A_{1}}(1)$ at the physical point, using the fit form of Eq. 20.

The fit function of Eq. (20) is readily applied, giving a $\chi^{2}/[\text{dof}]$ of 0.21 for 12 degrees of freedom. Figure 2 shows our results for $h^{s}_{A_{1}}(1)$ along with the fit function at zero lattice spacing and physical $u/d$, $s$ and $c$ quark masses as the grey band. Evaluating the fit function at the physical $b$ mass, as determined by $M_{\eta_{b}}$, gives our final result

[TABLE]

Adding the statistical and systematic errors in quadrature, we find a total fractional error of $1.5\%$. The error budget for this result is given in table 4. Note that we allow for an additional $\pm 10\,{\text{MeV}}$ uncertainty in the physical value of the $\eta_{b}$ mass beyond the experimental uncertainty, since our lattice QCD results do not include the effect of $\eta_{b}$ annihilation and QED McNeile et al. (2012b). This has no effect, however, since the heavy quark mass dependence is so mild.

Our total uncertainty is dominated by the statistical errors in our lattice results. The systematic error is dominated by that from the continuum extrapolation.

We include in Figure 2 the value from the only other lattice determination of $h^{s}_{A_{1}}(1)$ Harrison et al. (2018). This calculation also used MILC $n_{f}=2+1+1$ gluon field ensembles, but with the bulk of the ensembles used having coarser lattice spacing. This was made possible by the use of the NRQCD action for the $b$ quark Dowdall et al. (2012). The HISQ action was used for all the other quarks. The result of this calculation was: $h_{A_{1}}^{s}(1)=0.883(12)_{\text{stat}}(28)_{\text{sys}}$. Our result is in agreement with this, but with substantially smaller errors. The NRQCD uncertainty of 3.4% is dominated by the systematic error from the $\mathcal{O}(\alpha_{s})$ matching factor used to normalise the NRQCD-HISQ current and this error is absent from our calculation.

In addition to a value for $h^{s}_{A_{1}}(1)$ our calculation is able to give information on the physical dependence on the heavy quark mass of $\mathcal{F}^{H_{s}\rightarrow D^{*}_{s}}(1)$. We see from Figure 2 that this dependence is very mild to the point of being absent. We can determine the ratio of $\mathcal{F}^{H_{s}\rightarrow D^{*}_{s}}(1)$ for $m_{h}=m_{b}$ to $m_{h}=m_{c}$ (albeit that this latter point corresponds to an unphysical $D_{s}\rightarrow D_{s}^{*}$ decay) and find the value 0.998(23). Each of the terms (including $\eta_{A}$) in the HQET expectation of Eq. 21 can give effects of order a few percent to this ratio. The fact that we find no heavy quark mass dependence at the level of 2% shows that these effects must tend to cancel out.

The fit of our lattice results to Eq. (20) gives fit parameters $l_{V,A,P}$ which, as discussed in Section  II.5.1, provide a test of HQET. We find

[TABLE]

from our baseline fit. These results are compatible with values of $\mathcal{O}(\Lambda^{2}_{\text{QCD}})$ as expected by HQET. As discussed in Section II.5.1 these fit parameters change depending on the proxy that we use for the quark mass as well as our treatment of $\eta_{A}$. However, as we show in the tests performed in the next section (see Figure 5) this has little impact on our value for $h^{s}_{A_{1}}(1)$.

III.2 Further tests of our fit

Because we tune our $b$ and $c$ valence quark masses using the pseudoscalar heavyonium meson mass, we can independently test our results by comparing both our heavy-strange and $D_{s}^{*}$ meson masses against experiment. These results are shown in Figures 3 and 4. In each case we subtract half the corresponding pseudoscalar heavyonium mass to reduce lattice spacing uncertainties in the comparison to experiment Davies et al. (2010b).

Figure 3 shows that our $D_{s}^{*}$ meson mass agrees with experiment on all our ensembles at the level of our 5 MeV uncertainties. Systematic effects from missing QED and $\eta_{c}$ annihilation are expected to be of size a few MeV Davies et al. (2010b).

Figure 4 shows our results for the heavy-strange pseudoscalar meson mass as a function of the pseudoscalar heavyonium mass. We show the difference $\Delta_{h}=M_{H_{s}}-M_{\eta_{h}}/2$ to remove the leading $m_{h}$ dependence and also to reduce uncertainties from the value of the lattice spacing. We fit $\Delta_{h}$ to a simple function of $\varepsilon_{h}$ (Eq. (24)) :

[TABLE]

The leading, linear, term in $\varepsilon_{h}$ allows for the fact that the heavyonium ($\eta_{h}$) binding energy grows linearly with $m_{h}$ in a $1/r$ potential. We take priors on the $c_{i}$ of: $c_{-1}:0.05(5)$; $c_{0}:0.5(5)$; $c_{1}:0(1)$. $\mathcal{N}_{\text{disc}}$ takes the same form as in Eq. (28) with $a\Lambda_{\text{QCD}}$ (where $\Lambda_{\text{QCD}}$ is taken as 0.5 GeV) replacing $am_{c0}^{\text{val}}$, which is not relevant here. $\mathcal{N}_{\text{mistuning}}$ takes the same form as in Eq. (27).

Our result for the difference $M_{H_{s}}-M_{\eta_{h}}/2$ in the continuum at $m_{h}=m_{c}$ is 0.4755(37) GeV and at $m_{h}=m_{b}$ is 0.6588(61) GeV. These agree well with the earlier HPQCD results on $n_{f}=2+1$ gluon field configurations of 0.4753(22) GeV Davies et al. (2010b) and 0.658(11) GeV McNeile et al. (2012a). They also agree well with the experimental values of 0.4764(3) GeV and 0.6674(12) GeV Tanabashi et al. (2018), allowing for the $\sim$ 3–5 MeV effect from missing QED and $\eta_{b}$ and $\eta_{c}$ annihilation processes in the lattice QCD results.

We also performed a number of tests of our continuum/heavy-quark mass dependence fit to our results for $h^{s}_{A_{1}}(1)$. These are tabulated graphically in Figure 5.

One of the tests, denoted ‘Ratio with $f_{H_{c}}$’ in Figure 5, is described in more detail in Appendix A. It involves fitting the ratio of $h^{s}_{A_{1}}(1)$ to the $H_{c}$ decay constant, as a function of heavy quark mass and, after determining the continuum result at $m_{h}=m_{b}$, multiplying by the value for the $B_{c}$ decay constant determined from lattice QCD to obtain $h^{s}_{A_{1}}(1)$. The reason for doing this is because this ratio has smaller discretisation effects than $h^{s}_{A_{1}}(1)$ alone, as is clear from Figure 8 in Appendix A. It has stronger dependence on $m_{h}$, however, coming from the $H_{c}$ decay constant, along with sizeable uncertainties introduced from the uncertainty in the lattice spacing. Another disadvantage is that the physical result for $H_{c}$ decay constant must also be obtained. We find that this method gives results in agreement with our standard fit but with significantly larger uncertainties. It provides a good test, however, because it has very different $m_{h}$ dependence.

III.3 Implications for $B\to D^{*}$

As discussed in Section I, $h^{s}_{A_{1}}(1)$ is expected to be close in value to the equivalent $B\to D^{*}$ form factor, since they only differ in the mass of the light spectator quark and in effects arising from the strong decay of the $D^{*}$ to $D\pi$. In Harrison et al. (2018) the ratio of the two form factors was found to be: $h_{A_{1}}(1)/h^{s}_{A_{1}}(1)=1.013(14)_{\text{stat}}(17)_{\text{sys}}$. Note that systematic effects from the perturbative matching of the NRQCD-HISQ current largely cancel in this ratio.

Multiplying this by our result for $h^{s}_{A_{1}}(1)$, we can determine $h_{A_{1}}(1)$ as

[TABLE]

adding all the uncertainties in quadrature.

In Figures 6 and 7, we compare current lattice results for $h_{A_{1}}(1)$ and $h_{A_{1}}^{s}(1)$. Figure 6 compares final results for $h^{s}_{A_{1}}(1)$ from the HPQCD calculation using NRQCD $b$ quarks and HISQ lighter quarks Harrison et al. (2018) with our full HISQ result given here (Eq. (29)). It also compares final results for $h_{A_{1}}(1)$ from using the Fermilab approch Bailey et al. (2014) for $b$ and $c$ quarks and asqtad light quarks, NRQCD $b$ quarks and HISQ lighter quarks Harrison et al. (2018) and our result from Eq. 32 using the strange to light ratio from Harrison et al. (2018). Good agreement between all results is seen, well within the uncertainties quoted.

In Figure 7, we show more detail of the comparison by plotting the lattice results from the previous Fermilab/MILC Bailey et al. (2014) and NRQCD $b$ Harrison et al. (2018) calculations as a function of the valence spectator light quark mass (given by the square of the pion mass). Note that, for the results for $h_{A_{1}}(1)$ to the left of the plot, the valence light and sea masses are the same. For the $h^{s}_{A_{1}}(1)$ points from Harrison et al. (2018) to the right of the plot, the sea light (along with $s$ and $c$) quark masses take their physical values. Although agreement for $h_{A_{1}}(1)$ is seen at physical light quark mass in the continuum limit from all approaches, the NRQCD-HISQ results show systematic light quark mass dependence away from this point that is not visible in the Fermilab/MILC results. The two sets of results move apart as the spectator quark mass increases, and it is therefore not clear how well they would agree for spectator $s$ quarks.

Our results, shown in Figure 7 with black stars, agree with the NRQCD-HISQ results for $h^{s}_{A_{1}}(1)$. The smaller uncertainties from using a fully nonperturbative current normalisation here show that the perturbative matching uncertainty allowed for in Harrison et al. (2018) was conservative. Using the $s/l$ ratio from this calculation, where the perturbative matching uncertainty cancels, allows us to obtain an $h_{A_{1}}(1)$ result that agrees well with both earlier values. Our uncertainty on $h_{A_{1}}(1)$ is similar to that from Harrison et al. (2018) once we have combined the uncertainty from the ratio with that from our value for $h^{s}_{A_{1}}(1)$. However we have removed the perturbative matching uncertainty that dominates the NRQCD-HISQ error.

IV Conclusions

We have calculated the form factor at zero recoil, $\mathcal{F}^{B_{s}\to D^{*}_{s}}(1)$ or $h^{s}_{A_{1}}(1)$, using the relativistic HISQ formalism in full lattice QCD. This allows us to normalise the $b\rightarrow c$ current fully nonperturbatively for the first time and to determine how the form factor depends on the heavy quark mass (at physical charm quark mass). Our results show that dependence on the heavy quark mass is very mild (see Figure 2).

Our result

[TABLE]

agrees with an earlier lattice QCD result Harrison et al. (2018), but with half the uncertainty because of the nonperturbative normalisation of the current. Using the strange to light quark ratio from the earlier paper we are able to obtain a result for $\mathcal{F}^{B\rightarrow D^{*}}(1)$

[TABLE]

which is also free of perturbative matching uncertainties.

$h^{s}_{A_{1}}(1)$ will be a useful value to compare to experimental results in future to determine $V_{cb}$. It has some advantages from a lattice QCD perspective over $h_{A_{1}}$ as discussed in Section I. However, $h_{A_{1}}(1)$ can be combined with existing experimental results to obtain a value for the CKM element $V_{cb}$. The method of combination has been questioned recently when it was realised that the HQET constraints on the extrapolation of the exclusive experimental data to the zero recoil point were having a significant effect. Loosening these constraints gives a higher, but less precise, value for the combination $|\overline{\eta}_{\text{EW}}V_{cb}|h_{A_{1}}(1)$ (see, for example, the $V_{ub}/V_{cb}$ mini-review in Tanabashi et al. (2018)). Combining this experimental value with lattice QCD results for $h_{A_{1}}(1)$ then gives a result for $V_{cb}$ from the $B\rightarrow D^{*}\ell\nu$ exclusive decay that agrees with, but is less accurate than, that from inclusive $b\rightarrow c$ decays. We do not convert our $h_{A_{1}}$ result into a value for $V_{cb}$ here since it is clear from Figure 6 that we will agree with existing results (such as that in Harrison et al. (2018)) and, on its own, our new result does not have sufficient accuracy to reduce uncertainties in $V_{cb}$.

In future lattice QCD form factor calculations for both $B_{s}\rightarrow D_{s}^{*}$ and $B\rightarrow D^{*}$ need to work away from zero recoil to improve overlap with experimental results without the need for extrapolation111Preliminary results using the Fermilab formalism for $b$ and $c$ quarks and asqtad light quarks have already appeared Avilés-Casco et al. (2018), as have results using Möbius domain-wall quarks with a range of $m_{h}$ values up to 2.44$m_{c}$ Kaneko et al. (2018).. Our results here demonstrate the efficacy of HPQCD’s ‘heavy HISQ’ approach for form factors at zero recoil. Away from zero recoil we expect it to be even more useful because it is possible to map out the full $q^{2}$ range of the decay Lytle et al. (2016), where nonrelativistic approaches must stay close to zero recoil because of systematic errors that grow with the magnitude of the daughter meson momentum. Heavy HISQ calculations are underway for the form factors for $B_{(s)}\rightarrow D^{*}_{(s)}$, $B_{c}\rightarrow J/\psi$ decay, and the related $B_{s}\rightarrow D_{s}$ decay over the full $q^{2}$ range, using the techniques developed for $c\rightarrow s$ decays to normalise the currents nonperturbatively Koponen et al. (2013); Donald et al. (2014b). Initial results Colquhoun et al. (2016); McLean et al. (2018) look very promising.

Acknowledgements

We are grateful to the MILC collaboration for the use of their configurations and their code. Computing was done on the Cambridge Service for Data Driven Discovery (CSD3) supercomputer, part of which is operated by the University of Cambridge Research Computing Service on behalf of the UK Science and Technology Facilities Council (STFC) DiRAC HPC Facility. The DiRAC component of CSD3 was funded by BEIS via STFC capital grants and is operated by STFC operations grants. We are grateful to the CSD3 support staff for assistance. Funding for this work came from STFC. We would also like to thank C. Bouchard, B. Colquhoun, D. Hatton, J. Harrison, P. Lepage and M. Wingate for useful discussions.

Appendix A Ratio method for determining $h^{s}_{A_{1}}(1)$

It turns out that the significant discretisation effects visible in our results for $h^{s}_{A_{1}}(1)$ (Figure 2) are largely cancelled when we divide them by lattice QCD results for the decay constant of the heavy-charm pseudoscalar meson, $f_{H_{c}}$. This was also observed in Colquhoun et al. (2016) for vector form factors involving a $b\overline{c}$ current. $f_{H_{c}}$ is determined from the matrix element between the vacuum and the $H_{c}$ meson of the temporal axial vector $b\overline{c}$ current, whereas $h^{s}_{A_{1}}(1)$ is the matrix element between the $H_{s}$ and $D^{*}_{s}$ mesons of the spatial axial vector $b\overline{c}$ current. They behave very differently as a function of heavy quark mass but in practice have similar discretisation errors (compare Figures 2 and 9). We can make use of this in fitting the heavy quark mass dependence of their ratio with reduced discretisation effects. We also then need to fit the $H_{c}$ decay constant on its own in order to determine a physical value for the $B_{c}$ that we can use to determine $h^{s}_{A_{1}}(1)$ at the physical point.

$f_{H_{c}}$ is found using the PCAC relation for HISQ quarks

[TABLE]

where $\langle 0|P|H_{c}\rangle$ is determined in the fit to the $H_{c}$ two-point correlation functions via (14). We use a pseudoscalar operator, $P$, with spin-taste $\gamma_{5}\otimes\gamma_{5}$ so $f_{H_{c}}$ is absolutely normalised. Results for $f_{H_{c}}$ for each ensemble are given in Table 2 and plotted in Figure 9.

On each ensemble, at each heavy quark mass, we form the ratio $h^{s}_{A_{1}}(1)/(f_{H_{c}}\sqrt{M_{H_{c}}})$, plotted in Figure 8. Although discretisation effects largely cancel, the ratio varies strongly with changing heavy quark mass. This makes fitting this ratio as a function of heavy quark mass and lattice spacing very different to that of $h^{s}_{A_{1}}(1)$, with different systematic effects.

We use a fit function of the same form for both $h^{s}_{A_{1}}(1)/(f_{H_{c}}\sqrt{M_{H_{c}}})$ and $f_{H_{c}}$. Denoting the quantity being fit by $F$, we write (following McNeile et al. (2012b)):

[TABLE]

$\alpha_{s}(M)$ is the QCD coupling constant evaluated at scale $M$ and the ratio of $\alpha_{s}$ factors resums leading logarithms in HQET in the decay constant Neubert (1994). We take $\alpha_{s}$ in the $\overline{\text{MS}}$ scheme from lattice QCD Chakraborty et al. (2015). The power $p$ is then -6/25 (for $n_{f}=4$) for the $f_{H_{c}}$ fit and +6/25 for the fit to the ratio $h^{s}_{A_{1}}(1)/(f_{H_{c}}\sqrt{M_{H_{c}}})$. The leading power of $M_{\eta_{h}}$, $n$, is -1 for the fit to $f_{H_{c}}$ based on HQET expectations, but 0 for the fit to the ratio because we have used $f_{H_{c}}\sqrt{M_{H_{c}}}$ in the denominator to remove half-integer powers of $\varepsilon_{h}$ from the fit. The remainder of the fit function allows for inverse powers of $m_{h}$ and discretisation effects. $\mathcal{N}_{\text{mistuning}}$ is the same as that defined earlier for our $h^{s}_{A_{1}}$ fit and is given in Eq. (27). The final term allows for $c$ quark mistuning with prior on $c_{c}$ of $0\pm 1$. We take a prior on the overall constant $A$ of $0\pm 4(\text{GeV}^{3/2})$ in the $f_{H_{c}}$ fit and $0\pm 2(\text{GeV}^{-3/2})$ in the ratio fit. Priors on the $d_{ijk}$ are taken as $0\pm 2$ except for $d_{000}$ which is defined to have value 1.0.

The fit to the ratio is shown in Figure 8 and the fit to $f_{H_{c}}$ in Figure 9. For the ratio fit $\chi^{2}/\text{dof}$ is 0.27 for 12 degrees of freedom and for the $f_{H_{c}}$ fit, 0.53 for 16. Our final result for $f_{H_{c}}$ at $m_{h}=m_{b}$ agrees with a previous HPQCD heavy HISQ determination on gluon field configurations including $n_{f}=2+1$ flavours of sea quarks McNeile et al. (2012b) (shown as the red triangle in Figure 9). Our final result for the ratio $h^{s}_{A_{1}}(1)/(f_{H_{c}}\sqrt{M_{H_{c}}})$ at $m_{h}=m_{b}$ can then be multiplied by our value for $f_{B_{c}}$ and the square root of the mass of the $B_{c}$ meson from the Particle Data Tables Tanabashi et al. (2018), to give $h^{s}_{A_{1}}(1)$. This value is shown as the bottom point in Figure 5. Figure 8 compares the result from the ratio fit given by the grey band to the value (shown by inverted red triangle) obtained by taking our baseline fit result for $h^{s}_{A_{1}}(1)$ from Eq. (33) and calculating from it the value of the ratio $h^{s}_{A_{1}}(1)/(f_{H_{c}}\sqrt{M_{H_{c}}})$ using our value for $f_{B_{c}}$ and the experimental $M_{B_{c}}$. The agreement is good, showing the consistency of the two different approaches.