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

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

|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 K (where there is currently tension between the SM and experiment, see for example [42]).
Previous determinations of |V cb | have shown systematic discrepancies with each other. The two competing values were those derived from exclusive decays (B → D * ν * e.mclean.1@research.gla.ac.uk † christine.davies@glasgow.ac.uk ‡ URL: http://www.physics.gla.ac.uk/HPQCD and B → D ν with B → D * giving the more accurare result), and inclusive (B → X c ν, where X c is any charmed hadronic state). In 2016 the Heavy Flavour Averaging Group (HFLAV) gave a value derived from exclusive B → D * decays of |V cb | excl = (39.05 ± 0.47 exp ± 0.58 th ) × 10 −3 and from inclusive decays, using the kinetic scheme, of |V cb | incl = (42.19 ± 0.78) × 10 −3 [20]. It has since been suggested, based on unfolded Belle data [38], that the tension seen here arose (at least partly) from the use of a very constrained parameterization in the extrapolation of the experimental B → D * decay rates to zero recoil [43][44][45]. 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 | excl = (41.9 ± 2.0) × 10 −3 [46]. Although this means that inclusive and exclusive determinations now agree, it 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 → D * from lattice QCD is required, both at zero recoil and away from zero recoil. Another motivation for studying B → D * ν is the tension between SM and experimental determinations of the ratio R D ( * ) = B(B → D ( * ) τν τ )/B(B → D ( * ) ν ) ( = e or µ). The combined statistical significance of the anomalies in R D and R D * currently stands at 3.9σ [20]. More precise measurements and predictions will either confirm or dismiss a new physics explanation.
The weak decay process B s → D * s ν is very similar to B → D * ν 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 → D * . The absence of valence light quarks means that lattice QCD results have smaller statistical arXiv:1904.02046v1 [hep-lat] 3 Apr 2019 errors and are less computationally expensive. Finitevolume 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 [47] that can be considered 'gold-plated' in lattice QCD. This makes the B s → D * s ν both a useful test bed for lattice techniques (that may be later used to study B → D * ν 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 [48][49][50]. A combination of chiral perturbation theory and Heavy Quark Symmetry [51] backs up this expectation for B decays. We can therefore expect the form factors to be very similar for B s → D * s and B → D * . A recent lattice calculation [41] found an insignificant O(1%) difference at zero recoil: F Bs→D * s (1)/F B→D * (1) = 1.013 (14) stat (17) sys . Information from the study of B s → D * s can then be applicable to B → D * .
Lattice QCD calculations of the B (s) → 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 F B→D * (1) in [39,52] using the Fermilab action for both b and c quarks [53] and asqtad u/d quarks [54]. More recently the HPQCD collaboration computed both F B→D * (1) and F Bs→D * s (1) [41] using improved NRQCD b quarks [55,56] and Highly Improved Staggered (HISQ) c and u/d/s quarks [57].
The RBC/UKQCD [58] and LANL-SWME [59] 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 [60].
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 O(α s ) for these actions [61,62]. 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 ∼ 80% of the final error in the HPQCD calcu-lation [41] and ∼ 30% in the Fermilab/MILC calculation [39] because of the differing allowances for missing higher-order terms.
Here we report details and results of a calculation of the B s → 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 [63,64], 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 [65,66] and the b quark mass [67,68]. It is also now being used by other collaborations for these calculations [69,70]. A proof-of-principle application of heavy-HISQ to form factors was given in [71,72] 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.

A. Form Factors
The differential decay rate for theB 0 s → D * + s l −ν l decay is given in the SM by χ(1) = 1 (see, for example, appendix G of [41]).η EW accounts for electroweak corrections from diagrams where photons or Zs are exchanged in addition to a W − , as well as the Coulomb attraction of the final-state charged particles [73][74][75]. The differential decay rate for the B 0 s → D * − s l +ν l is identical. The form factor F Bs→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 [76] is the axial-vector current. 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 and F Bs→D * s (w) reduces to Our goal is to compute h s A1 (1). All we need to do this is the matrix element D * s ( )|A µ |B s with both the B s and D * s at rest, with the D * s polarization in the same direction as the (spatial) axial-vector current.

B. Lattice Calculation
The gluon field configurations that we use were generated by the MILC collaboration [63,64]. Table I 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 O(α s a 2 , n f α s a 2 ) [77]. 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 ≡ m l ), using the HISQ action [57]. 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 A1 (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 s A1 (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 val h0 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' [57]. They differ in the amount of point-splitting between the fields and the space-time dependent phase needed to substitute for the appropriate γ 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, Γ P = (γ 5 ⊗ γ 5 ); vector, Γ µ V = (γ µ ⊗ γ µ ) and axial-vector, Γ µ A = (γ µ γ 5 ⊗ γ µ γ 5 ). We compute several 'two-point' correlation functions on the ensembles detailed in table I, combining HISQ propagators from solving the Dirac equation for each random wall time source. These correlation functions take the form where represents a functional integral, q, q are valence quark fields of the flavours the M meson is charged under, Γ is the spin-taste structure of M and 1/N taste is the staggered quark normalisation for closed loops. The random-wall source and the sum over x at the sink project onto zero spatial momentum. We compute the correlation functions for all t values, i.e. 0 ≤ t ≤ N t .
The correlation function for the heavy-strange pseudoscalar meson, H s , with valence quark content hs and spin-taste structure Γ P is constructed from HISQ propagators as: 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 x, y. We also compute correlators for a charm-strange vector meson D * s , with structure Γ i V , using 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Ĥ c respectively. H c has spin-taste structure Γ P , andĤ c has structure Γ 0 A . H c correlators are computed using Eq. (7)   [63,64]. a is the lattice spacing, determined from the Wilson flow parameter, w0. Values for w0/a are from: set 1, [78], sets 2 and 3, [68] and set 4 [79]. The physical value of w0 was determined at 0.1715(9) fm from fπ [80]. Nx is the spatial extent and Nt the temporal extent of the lattice in lattice units; n cfg is the number of gluon field configurations in the ensemble and nsrc 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 and c valence quarks were tuned in [68]. We use a number of heavy quark masses to assist the extrapolation to the physical b mass. Column 12 gives the temporal separations between source and sink, T , of the 3-point correlation functions computed on each ensemble.
(with g s replaced with g c ), whileĤ c correlators are given by where we use the notationz µ = ν =µ z ν . These correlators will be used to normalise the axial vector bc current as discussed in Section II D.
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 η h , with spin-taste structure Γ P using Eq. (7). Likewise, to test the impact of any mistuning of the charm and strange quark masses, we also determine η c and η s correlators analogously. We can tune the c and b masses using the experimental values for the η c and η b masses, allowing for slight shifts from missing QED effects and the fact that we do not allow these mesons to annihilate to gluons [67]. The mass of the η s meson (which is not a physical state) can be fixed in lattice QCD from the K and π meson masses [80,81].
We then generate the 'three-point' correlation functions that contain the H s to D * s transition.
Our H s source is given spin-taste Γ P , the D * s sink, Γ i V , and the current insertion Γ i A . This gives the required cancellation of tastes within the three-point function [82]. In terms of HISQ propagators 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 ≤ t ≤ T , and 3 T values that vary between ensembles and are given in Table I. We average over the 3 directions for i for increased statistical precision.

C. Analysis of Correlation Functions
We use simultaneous Bayesian fits [83,84] 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 E.
We fit the two-point correlation functions using the functional form where N t is the temporal extent of the lattice, and E are fit parameters, with the excited-state energy parameters implemented as energy differences to the state below [83]. The second term accounts for the presence of opposite-parity states that contribute an oscillating term to the correlation function when using staggered quarks [57]. These terms do not appear when M is a pseudoscalar with a quark and antiquark of the same mass, so in the M = η h , η c , and η s cases the second term is not required. For all correlator fits we set N 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 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 , η h , η s , H c ,Ĥ c , and threepoint) 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 M n+1 − E M n and amplitudes a M n (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 M 0 are given priors of (am q0 + am q 0 + aΛ QCD ) ± 2aΛ QCD , where m q0 and m q 0 are the masses of the appropriate quarks, and Λ 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 M,o 0 , are given priors of (am q0 + am q 0 + 2aΛ QCD ) ± 2aΛ QCD . Excited state energy differences, E M i+1 − E M i , i > 0 are given prior values 2aΛ QCD ± aΛ QCD . Priors for ground state amplitudes a M 0 , 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 3point parameter, J nn 00 is given a prior of 1 ± 0.6, and the rest of the 3-point parameters J nn jk are given 0 ± 1.
The (as yet unnormalised) matrix element that we need to obtain h s A1 (1) is given by To ensure that truncating the sum over states at N 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 ≥ t cut and t ≤ N t − t cut in the two-point case and t ≤ T − t cut in the three-point case. We can in principle use a different t cut for every correlation function included in our fit, but we do not use a big range of t 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 χ 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 [84]. This replaces correlation matrix eigenvalues below λ min , equal to svdcut times the largest eigenvalue, with λ min . λ min is estimated using the diagnosis tools in the Corrfitter package [84] and corresponds typically to an svdcut of 10 −3 here. Figure 1 summarises stability tests of our fits, focussing on the key parameter J nn 00 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 A1 (1) are given in Table II. Notice that the statistical errors on the results grow with the heavy quark mass. This is a well under-stood problem in lattice heavy-light meson physics (see, for example [85]). Our method here has the advantage of including information from lighter-than-b heavy quarks with improved statistical precision.

D. Normalisation of the Axial Current
The partially-conserved axial-vector current for the HISQ action is a complicated linear combination of onelink 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 [82].
We use the fact that the staggered local pseudoscalar current of spin-taste (γ 5 ⊗ γ 5 ), multiplied by the sum of its valence quark masses, is absolutely normalized via the PCAC relation. From the two-point H c andĤ c correlator fits we can extract the decay amplitudes: 0|c(γ 5 ⊗ γ 5 )h|H c ≡ 0|P |H c and 0|c(γ 0 γ 5 ⊗ γ 0 γ 5 )h|Ĥ c = 0|A 0 |Ĥ c 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 The Z A values found on each ensemble and am val h0 are given in Table III. There is an ambiguity in what mass to use on the right hand side of Eq. (16). We use the non-goldstone mass MĤ c , but one could just as well replace this with M Hc since the difference is a discretisation effect. The meson mass difference is very small for heavy mesons [57] 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 [57] by multiplying by a factor Z disc . This is derived in [62] as: C q = cosh am q,tree 1 − 1 + q,Naik 2 sinh 2 am q,tree .
See also [69]. m q,tree is the tree-level pole mass in the HISQ action. It has an expansion in terms of the bare mass [57] am q,tree = am q0 1 − am 10 q0 + O(am 12 q0 ) , q,Naik fixes the Naik parameter [86] (N = 1+ is the coefficient of the tree-level improvement term for the derivative) in the HISQ action when it is being used for heavy quarks [57]. q,Naik is set to its tree-level value, removing the leading tree-level errors from the dispersion relation. As an expansion in am q,tree it begins at O(am q,tree ) 2 [57].
To determine q,Naik we use the closed form expression for it given in [62] and this can also be used along with Eq. (18) to evaluate Z disc . The pole condition can be used to show that the expansion ofC q begins at am 4 q0 as 1 − 3am 4 q0 /80 + . . .. The effect of Z disc is then very small, never exceeding 0.2%. Z disc values on each ensemble for each am val h0 are given in table III. 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:

E. Obtaining a Result at the Physical Point
We now discuss how we fit our results for the zero recoil form factor, h s A1 (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 s A1 (1) to the following form The terms in the first line allow for dependence on the valence heavy quark and charm quark masses (with ε q ≡ 1/m q ) using input from HQET, to be discussed below. N disc and N mistuning account for discretisation and mass mistuning effects, also discussed below. The physical result is then h s A1 (1)(0, m l,phys , m b ).

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  [87]. The HQET expression for h A1 (1) is then given by [88,89]:   where l V , l A and l P are O(Λ 2 QCD ) (with possible mild dependence on whether the spectator quark is s or u/d). η A accounts for ultraviolet matching between HQET and QCD, and has been computed to 2-loops in perturbative QCD [90]. It has mild dependence on m h through logarithms of m c /m h ; at one-loop η A has explicit form [91] The coefficient of α s /π 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) [92]. The two-loop correction is small [90]. η 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 η A at m h = m c . The terms in round brackets in Eq. (21), multiplying η 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 η A with appropriate value for α s (since logarithmic m h dependence of η 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 η 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 ± 1 GeV 2 .
We have several different proxies, derived from heavy meson masses, that we can take for the heavy quark mass that appears in ε h in Eq. (20). We do not expect our physical result for h s A1 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 Hs , and half the mass of the pseudoscalar heavyonium meson, M η h . We also tested using the quark mass in the minimal renormalon subtracted (MRS) scheme suggested in [93]. This takes where µ 2 S = µ 2 π,S − d H ( * ) µ 2 G,S with d H ( * ) = 1 for pseudoscalar mesons and −1/3 for vectors. For this case we use parameters determined in [92] for the MRS scheme: Λ MRS = 0.552(30)GeV, µ 2 π,MRS = 0.06(22)GeV 2 and µ 2 G,MRS = 0.38(1)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 We test the stability of the fit results under the different choices discussed above in Section III B.

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 η c and η s meson masses to their physical values in a pure QCD world allowing, for example, for η c annihilation as discussed in [68]. 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 ε c from the physical η c mass.
The strange (valence and sea) and light (sea) mass mistunings are accounted for using the tuning in [68].
M physical ηs is determined in lattice simulations from the masses of the pion and kaon [80]. The ratio δ s /m tuned s 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 δ l = m l0 − m tuned l . We find m tuned l from m tuned s , leveraging the fact that the ratio of quark masses is regularization independent, and the ratio was calculated in [69]: m s m l phys = 27.18 (10).
We set m tuned l to m tuned s divided by this ratio. The full term we include to account for mistuning is then given by where c l , c s and c val s are fit parameters with prior distributions 0 ± 1. We neglect δ 2 s,l 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 ≡ 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 tuned l value up and down by the expected value for m d − m u [46]. We find the effect to be completely negligible in comparison to the other sources of error.

Discretisation Effects
Discretisation effects in our lattice QCD results are accounted for following the methodology of [66]. We take 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 η h being used here as a proxy for the heavy quark mass. d ijk are fit parameters with prior distributions 0 ± 1.0. All discretisation effects are of even order by construction of the HISQ action [57]. We tested the impact on the fit of including extra discretisation effects set by the scale Λ QCD but this made no difference (since such effects are much smaller than those already included by the am val c0 terms). We also tested the effects of increasing the number of terms in each sum, but the final result remained unchanged.

Finite-volume Effects
The finite volume effects in our lattice results are expected to be negligible, because we are working with   (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 uncertainty in the physical value of the η b mass has a negligible effect.
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 [94] is therefore absent and we can safely ignore finite-volume effects here.
In section III B we detail the results of several tests of the stability of our final result under changes to the details of the fit.

III. RESULTS AND DISCUSSION
The results of our correlation function fits (discussed in Section II C) are given in Table II. We tabulate values for h s A1 (1) at each heavy quark mass that we have used on each gluon field ensemble from Table I  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 IV. Note that we allow for an additional ±10 MeV uncertainty in the physical value of the η b mass beyond the experimental uncertainty, since our lattice QCD results do not include the effect of η b annihilation and QED [66]. 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 A1 (1) [41]. This calculation also  Table I. Solid lines simply join the points on a given ensemble for added clarity. The red inverted triangle gives the determination of the same quantity from a previous study using the NRQCD action for the b quark [41].
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 [56]. The HISQ action was used for all the other quarks. The result of this calculation was: h s A1 (1) = 0.883(12) stat (28) 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 O(α 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 A1 (1) our calculation is able to give information on the physical dependence on the heavy quark mass of F Hs→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 F Hs→D * s (1) for m h = m b to m h = m c (albeit that this latter point corresponds to an unphysical D s → D * s decay) and find the value 0.998 (23). Each of the terms (including η 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 E 1, provide a test of HQET. We find l s V = 0.71(28)GeV 2 , l s A = −0.34 (32) l s P = −0.53(34)GeV 2 , from our baseline fit. These results are compatible with values of O(Λ 2 QCD ) as expected by HQET. As discussed in Section II E 1 these fit parameters change depending on the proxy that we use for the quark mass as well as our treatment of η 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 A1 (1).

B. 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 [85]. 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 η c annihilation are expected to be of size a few MeV [85].  shows our results for the heavy-strange pseudoscalar meson mass as a function of the pseudoscalar heavyonium mass. We show the difference ∆ h = M Hs − M η h /2 to remove the leading m h dependence and also to reduce uncertainties from the value of the lattice spacing. We fit ∆ h to a simple function of ε h (Eq. (24)) : × (1 + N disc + N mistuning ). The leading, linear, term in ε h allows for the fact that the heavyonium (η 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). N disc takes the same form as in Eq. (28) with aΛ QCD (where Λ QCD is taken as 0.5 GeV) replacing am val c0 , which is not relevant here. N mistuning takes the same form as in Eq. (27).
Our result for the difference M Hs − M η 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 [85] and 0.658(11) GeV [65]. They also agree well with the experimental values of 0.4764(3) GeV and 0.6674(12) GeV [46], allowing for the ∼ 3-5 MeV effect from missing QED and η b and η 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 A1 (1). These are tabulated graphically in Figure 5.
One of the tests, denoted 'Ratio with f Hc ' in Figure 5, is described in more detail in Appendix A. It involves fitting the ratio of h s A1 (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 A1 (1). The reason for doing this is because this ratio has smaller discretisation effects than h s A1 (1) alone, as is clear from Figure 8 (20) of the form p/M 3 η h where p is a fit parameter with the same prior as l s V,A,P . In this case the Bayes factor falls by a factor of 7, suggesting that the results do not contain a cubic dependence on the heavy mass. The next two points show the results of including specific implementations of ηA described in Section II E (rather than the value 1). In the upper variant parameter ρ is given prior 0 ± 1. 'A + 1/m b mc + 1/m 2 b ' is the result of replacing 1 + lV /m 2 c in the fit with a fit parameter A with prior distribution 1 ± 1. The fact that this does not affect the fit shows that mistuning of the charm quark mass is a negligible effect here. The points with labels beginning 'ε h =' show the result of replacing the heavy mass proxy Mη h /2 with MH s and the MRS quark mass (Eq. (23)) respectively. The bottom point labelled 'Ratio with fH c ' is the result of an alternative extrapolation described in Appendix A.
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.  [41] and the value marked '(Fermilab, Fermilab/MILC)' is from [39].
As discussed in Section I, h s A1 (1) is expected to be close in value to the equivalent B → 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π. In [41] the ratio of the two form factors was found to be: h A1 (1)/h s A1 (1) = 1.013 (14) stat (17) 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 A1 (1), we can determine h A1 (1) as adding all the uncertainties in quadrature. In Figures 6 and 7, we compare current lattice results for h A1 (1) and h s A1 (1). Figure 6 compares final results for h s A1 (1) from the HPQCD calculation using NRQCD b quarks and HISQ lighter quarks [41] with our full HISQ result given here (Eq. (29)). It also compares final results for h A1 (1) from using the Fermilab approch [39] for b and c quarks and asqtad light quarks, NRQCD b quarks and HISQ lighter quarks [41] and our result from Eq. 32 using the strange to light ratio from [41]. 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 [39] and NRQCD b [41] 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 A1 (1) to the left of the plot, the valence light and sea masses are the same. For the h s A1 (1) points from [41] [41] and [39] and are plotted as a function of valence (=sea) light quark mass, given by the square of Mπ. On the right are points for h s A 1 (1) from [41] plotted at the appropriate valence mass for the s quark, but obtained at physical sea light quark masses. The final result for hA 1 (1) from [39], with its full error bar, is given by the inverted blue triangle. The inverted red triangles give the final results for hA 1 (1) and h s A 1 (1) from [41]. Our results here are given by the black stars.
to the right of the plot, the sea light (along with s and c) quark masses take their physical values. Although agreement for h A1 (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 A1 (1). The smaller uncertainties from using a fully nonperturbative current normalisation here show that the perturbative matching uncertainty allowed for in [41] was conservative. Using the s/l ratio from this calculation, where the perturbative matching uncertainty cancels, allows us to obtain an h A1 (1) result that agrees well with both earlier values. Our uncertainty on h A1 (1) is similar to that from [41] once we have combined the uncertainty from the ratio with that from our value for h s A1 (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, F Bs→D * s (1) or h s A1 (1), using the relativistic HISQ for-malism in full lattice QCD. This allows us to normalise the b → 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 agrees with an earlier lattice QCD result [41], 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 F B→D * (1) which is also free of perturbative matching uncertainties. h s A1 (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 A1 as discussed in Section I. However, h A1 (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 |η EW V cb |h A1 (1) (see, for example, the V ub /V cb minireview in [46]). Combining this experimental value with lattice QCD results for h A1 (1) then gives a result for V cb from the B → D * ν exclusive decay that agrees with, but is less accurate than, that from inclusive b → c decays. We do not convert our h A1 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 [41]) 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 → D * s and B → D * need to work away from zero recoil to improve overlap with experimental results without the need for extrapolation 1 . 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 [71], 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) → D * (s) , B c → J/ψ decay, and the related B s → D s decay over the full q 2 range, using the techniques developed for c → s decays to normalise the currents nonperturbatively [49,82]. Initial results [72,96] Table I. The grey band shows the result of the fit described in the text, evaluated at a = 0 and physical l, s and c quark masses to give the physical heavy quark mass dependence of the ratio. At m h = m b we obtain the result given by the black star. For comparison with our previous fit for h s A 1 (1) the inverted red triangle shows our result from Eq. (33) converted to a ratio using the value for fB c from Figure 9 and MB c from experiment [46].
It turns out that the significant discretisation effects visible in our results for h s A1 (1) (Figure 2) are largely The heavy-charm pseudoscalar meson decay constant, fH c , plotted against Mη h (a proxy for the heavy quark mass). Gluon field ensembles listed in the legend follow the order of sets in Table I. The grey band shows the result of the fit described in the text, evaluated at a = 0 and physical l, s and c quark masses to give the physical heavy quark mass dependence of the decay constant. At m h = m b we obtain the result given by the black star. The red triangle shows the result from a previous heavy HISQ determination of fB c on n f = 2 + 1 gluon field ensembles [66].
cancelled when we divide them by lattice QCD results for the decay constant of the heavy-charm pseudoscalar meson, f Hc . This was also observed in [72] for vector form factors involving a bc current. f Hc is determined from the matrix element between the vacuum and the H c meson of the temporal axial vector bc current, whereas h s A1 (1) is the matrix element between the H s and D * s mesons of the spatial axial vector bc 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 A1 (1) at the physical point.
f Hc is found using the PCAC relation for HISQ quarks where 0|P |H c is determined in the fit to the H c twopoint correlation functions via (14). We use a pseudoscalar operator, P , with spin-taste γ 5 ⊗ γ 5 so f Hc is absolutely normalised. Results for f Hc for each ensemble are given in Table II and plotted in Figure 9. On each ensemble, at each heavy quark mass, we form the ratio h s A1 (1)/(f Hc M Hc ), 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 A1 (1), with different systematic effects.
We use a fit function of the same form for both h s A1 (1)/(f Hc M Hc ) and f Hc . Denoting the quantity being fit by F , we write (following [66]):  [97]. We take α s in the MS scheme from lattice QCD [68]. The power p is then -6/25 (for n f = 4) for the f Hc fit and +6/25 for the fit to the ratio h s A1 (1)/(f Hc M Hc ). The leading power of M η h , n, is -1 for the fit to f Hc based on HQET expectations, but 0 for the fit to the ratio because we have used f Hc M Hc in the denominator to remove halfinteger powers of ε h from the fit. The remainder of the fit function allows for inverse powers of m h and discretisa-tion effects. N mistuning is the same as that defined earlier for our h s A1 fit and is given in Eq. (27). The final term allows for c quark mistuning with prior on c c of 0±1. We take a prior on the overall constant A of 0 ± 4(GeV 3/2 ) in the f Hc fit and 0 ± 2(GeV −3/2 ) in the ratio fit. Priors on the d ijk are taken as 0 ± 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 Hc in Figure 9. For the ratio fit χ 2 /dof is 0.27 for 12 degrees of freedom and for the f Hc fit, 0.53 for 16. Our final result for f Hc 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 [66] (shown as the red triangle in Figure 9). Our final result for the ratio h s A1 (1)/(f Hc M Hc ) at m h = m b can then be multiplied by our value for f Bc and the square root of the mass of the B c meson from the Particle Data Tables [46], to give h s A1 (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 A1 (1) from Eq. (33) and calculating from it the value of the ratio h s A1 (1)/(f Hc M Hc ) using our value for f Bc and the experimental M Bc . The agreement is good, showing the consistency of the two different approaches.