Lattice QCD calculation of the ${{B}_{(s)}\to D_{(s)}^{*}\ell{\nu}}$ form factors at zero recoil and implications for ${|V_{cb}|}$

We present results of a lattice QCD calculation of $B\to D^*$ and $B_s\to D_s^*$ axial vector matrix elements with both states at rest. These zero recoil matrix elements provide the normalization necessary to infer a value for the CKM matrix element $|V_{cb}|$ from experimental measurements of $\bar{B}^0\to D^{*+}\ell^-\bar{\nu}$ and $\bar{B}^0_s\to D_s^{*+}\ell^-\bar{\nu}$ decay. Results are derived from correlation functions computed with highly improved staggered quarks (HISQ) for light, strange, and charm quark propagators, and nonrelativistic QCD for the bottom quark propagator. The calculation of correlation functions employs MILC Collaboration ensembles over a range of three lattice spacings. These gauge field configurations include sea quark effects of charm, strange, and equal-mass up and down quarks. We use ensembles with physically light up and down quarks, as well as heavier values. Our main results are $\mathcal{F}^{B\to D^*}(1)= 0.895\pm 0.010_{\mathrm{stat}}\pm{{0.024}_{\mathrm{sys}}}$ and $\mathcal{F}^{B_s\to D_s^*}(1)= 0.883\pm 0.010_{\mathrm{stat}}\pm{0.028_{\mathrm{sys}}}$. We discuss the consequences for $|V_{cb}|$ in light of recent investigations into the extrapolation of experimental data to zero recoil.


I. INTRODUCTION
Precise measurements of quark flavour-changing interactions offer one way to uncover physics beyond the Standard Model. As successful as the Standard Model appears to be so far, there will continue to be progress reducing experimental and theoretical uncertainties, as well as making new measurements. Existing tensions in the global fits to the Cabibbo-Kobyashi-Maskawa (CKM) parameters may become outright inconsistencies, or new measurements of rare decays may differ significantly from Standard Model predictions.
Measurements of the inclusive b → c decays B → X c ν, combined with an operator product expansion offer a complementary method. The latest estimate is |V cb | = (42.21±0.78)×10 −3 [19,20]. The discrepancy between the inclusive and exclusive result described above is at the 3σ level.
Recently it has been suggested that the inclusive/exclusive difference could be due to modela jgihh2@cam.ac.uk b http://www.physics.gla.ac.uk/HPQCD dependence implicit in extrapolating experimental data for B → D * ν to the zero recoil point. The CLN parametrization [21] has been used in recent analyses since it takes advantage of heavy quark symmetries to improve unitarity constraints in the form factor shape function. This had several advantages for some time, but with increased precision in the experimental data, it is possible that uncertainties arising from these constraints are no longer negligible. In fact, recent work [22][23][24][25][26][27] has shown that replacing the CLN parametrization by the BGL parametrization [28] yields a determination of |V cb | which is as much as 10% higher, in much better agreement with the |V cb | from inclusive decays.
One can also use the exclusive decay B → D ν to estimate |V cb |. Historically this has not given as precise a determination due to having to contend with background from B → D * ν. Recent progress has come from new measurements and joint fits to experimental and lattice [29,30] data over a range of q 2 using socalled z-parametrizations [31,32]. The latest result using B → D ν results is |V cb | = (40.85 ± 0.98) × 10 −3 [33], in acceptable agreement with either the B → D * ν or B → X c ν determinations.
It is worth noting that a determination of |V cb | is important beyond semileptonic b → c decays. Due to insufficient direct knowledge of the top-strange coupling, Standard Model predictions which depend on V ts rely on CKM unitarity, and therefore on V cb . For example the K 0 -K 0 mixing parameter K depends sensitively on V cb ; taking sin 2β as an input, then K ∝ |V cb | 4 at leading order [34].
In this article we present the details and results of a lattice calculation of the zero-recoil form factor needed to extract |V cb | from experimental measurements of the B → D * ν and B s → D * s ν decay rates. This work differs from the Fermilab/MILC calculation [18] in the following respects: (1) the gauge field configurations are the next generation MILC ensembles [35][36][37] which include effects of 2 + 1 + 1 flavours of sea quarks using the highly improved staggered quark (HISQ) action [38]; (2) the fully relativistic HISQ action is used for valence light, strange, and charm quarks; (3) the nonrelativistic QCD (NRQCD) action [39] is used for the bottom quark. Therefore, this work represents a statistically independent, complementary calculation to Ref. [18], with different formulations in many respects. The two main advantages of using the HISQ action is that discretization errors are reduced and that the MILC HISQ ensembles include configurations with physically light u/d quark effects. We reported preliminary results in recent proceedings [40].
Other groups are applying different methods to calculate B (s) → D ( * ) (s) form factors. Two-flavour twisted-mass configurations have been used to estimate the B s → D s form factors near zero recoil [41]; however the uncertainties with this formulation are quite large. Work has also recently begun using the domain wall action for light, strange, and charm quarks [42]. Having results for the form factors from several groups, each using different approaches, would be very helpful and could lead to a further reduction in uncertainties by allowing global fits to uncorrelated numerical and experimental data.
This paper is structured as follows. Sec. II briefly introduces the hadronic matrix elements of interest and sets some notation. In Sec. III we list the inputs to our computation and summarize the correlation functions calculated. The matching between lattice and continuum currents is discussed in Sec. IV. Sec. V is the most important section for the expert reader; there we discuss the fits to the correlation functions, the treatment of discretization and quark mass errors, and estimates of other systematic uncertainties. We summarize the result of the lattice calculation in Sec. VI. In Sec. VII we investigate the implications of the new lattice result in the context of renewed scrutiny of the extrapolation of experimental data to zero recoil; there we propose a simplified series expansion as the one least likely to introduce hidden theoretical uncertainties into a form factor parametrization. We offer brief conclusions in Sec. VIII. Several appendices are provided which contain further definitions and details in hopes of making the manuscript as self-contained as possible. These are noted at appropriate places in the body of the paper.
A reader more interested in the results and consequences than the details of our calculation can safely focus a first reading on Sections II, VI, VII, and VIII, possibly referring to Appendix G.

II. FORM FACTORS
This section simply summarizes standard notation relating the differential decay rate, the relevant hadronic matrix elements, and the corresponding form factors. Throughout the section we refer toB 0 → D * + −ν de-cay, but the expressions forB 0 s → D * + s −ν are the same, mutatis mutandis.
The differential decay rate, integrated over angular variables, is given in the Standard Model by where w = v · v is the scalar product of the B and D * 4-velocities, and χ(w) is a known kinematic function normalized so that χ(1) = 1 (see App. G). The coefficient η EW accounts for electroweak corrections due to box diagrams in which a photon or Z boson is exchanged in addition to a W boson as well as the Coulomb attraction of the final-state charged particles [43][44][45]. The form factor F(w) is a linear combination of hadronic form factors parametrizing the matrix elements of the V −A weak current, i.e.
The only contribution to F(w) at zero recoil, w = 1, is from the matrix element of the axial vector current; this reduces to D * (p, )|cγ j γ 5 b|B(p) = (M B + M D * )A 1 (q 2 max ) * j (3) for j = 1, 2, 3. It is sometimes conventional to work with form factors defined within heavy quark effective theory (HQET). Of relevance to this work, we write At zero recoil, where w = 1 and q 2 = q 2 max , For brevity throughout the paper we will usually use the h A1 notation. When we wish to specifically refer to the B s → D * s form factor, we write h s A1 , so These are the quantities we calculate here.

III. LATTICE PARAMETERS AND METHODOLOGY
Here we give specific details about the lattice calculation. Once again many of the expressions will refer to I: Details of the gauge configurations used in this work. We refer to sets 1, 2 and 3 as 'very coarse', sets 4, 5 and 6 as 'coarse' and sets 7 and 8 as 'fine'. The lattice spacings were determined from the Υ(2S − 1S) splitting in [46]. Sets 3, 6 and 8 use light quarks with their physical masses. u0 is the tadpole improvement factor, here we use the Landau gauge mean link. The final column specifies the total number of configurations multiplied by the number of different start times used for sources on each. In order to improve statistical precision we use random wall sources.
Set a(fm) L/a × T /a am l ams amc u0 n cfg × nt  [47] and the b mass was taken from [46].
(1 + Naik ) is the coefficient of the charm Naik term and ci are the perturbatively improved coefficients appearing in the NRQCD action correct through O(αsv 4 ) [46]. The last column gives the T values used in three point functions. These have changed from those presented in [40] on the very coarse ensembles as it was found that T = 10, 11, 12, 13 resulted in excessive noise on Set 3, which resulted in poor fit stability and the relatively low value of F(1) on this ensemble. B → D * matrix elements, but they apply for any spectator quark mass. The gluon field configurations that we use were generated by the MILC collaboration and include 2+1+1 flavours of dynamical HISQ quarks in the sea and include 3 different lattice spacings [35][36][37]. The u and d quarks have equal mass, m u = m d ≡ m l , and in our calculations we use the values m l /m s = 0.2, 0.1 and the physical value 1/27.4 [48]. The s and c quarks in the sea are also well-tuned [47] and included using the HISQ action. The gauge action is the Symanzik improved gluon action with coefficients correct to O(α s a 2 , n f α s a 2 ) [49]. Table I gives numerical values for the lattice spacings, quark masses, and other parameters describing the ensembles we used.
In calculating correlation functions, we slightly tune the valence s and c masses closer to their physical values. The d, s, and c quark propagators were computed using the MILC code [50]. The b quark is simulated using perturbatively improved non-relativistic QCD [46,51], which takes advantage of the non-relativistic nature of the b quark dynamics in B mesons and produces very good control over discretization uncertainties. Details of the gauge, NRQCD, and HISQ actions used are given in Appendices A, B, and C, respectively. In Table II we record the parameters used in calculating quark propagators.
In order to extract the form factor from lattice calculations we must compute the set of Euclidean correlation functions where each interpolating operator O i is projected onto zero spatial momentum by summing over spatial lattice points and the current J κ is one of several lattice currents (see Sec. IV). The indices i and j label different smearing functions. We use three different smearing operators on each of the B and D * interpolating operators. In implementing O µ (t) D * i we use an unsmeared operator and two gauge covariant Gaussian smearings, implemented by applying 1 − r 2 D * ∇ 2 n n to the field. Here the derivative is stride-2 in order not to mix the staggeredtaste meson multiplets. r D * is the radius (in lattice units) chosen to give good overlap with the ground state, and n is chosen to give a good approximation to a Gaussian while maintaining numerical stability. For the B we use a local operator as well as two Gaussian smearings, implemented as 1 N e −(x−y) 2 /r 2 B , where again r B is a radius in lattice units and N is an overall normalization. Since the B smearings are not gauge invariant, the gauge fields are fixed to Coulomb gauge. We refer to the local operator as l and the Gaussian smearings as g2 and g4 corresponding to radii of 2a and 4a respectively. We use the same choices of radii for both B and D * smearings. The smearing parameters are given in Table III. The interpolating operators themselves are where ∆(x, y) is the appropriate smearing function discussed above. In distinction to the continuum quark fields b, c, s, and u of Sec. II, here we denote the NRQCD b field by Ψ b and the staggered fields, written as 4component Dirac spinors (see App. C), by ψ with the appropriate flavour subscript. We checked both the point-split and local D * interpolating operators on the very coarse, physical point ensemble (Set 3) and found no significant difference in statistical noise or central value of either the D * mass or the matrix element. We primarily used the point-split current as it was simpler to implement in our framework. The results quoted below for the B → D * fits use the point-split vector current, except for Set 3 where results are given for the local vector current. The results below for B s → D * s form factors were obtained using the local vector current.
In order to improve statistics we multiply our smeared sources with random walls to produce, on average, the effect of multiple sources. Taking the all-to-all 2-point function as an example we have Exact computation requires an inversion for each value of y being summed over. Instead we generate a random vector ξ satisfying lim N →∞ N l ξ al (x)ξ bl (y) * = δ(x, y)δ ab .
N here is the number of random vector wall sources. The average over configurations further suppresses violations of this relation; in practice a single random wall per colour, setting N = N c = 3, is sufficient. Inserting the above relation into the 2-point function where we have used γ 5 hermiticity. The naive propagators G are built from staggered quarks and the full form of the correlation function contractions in terms of NRQCD and staggered propagators is given in Appendix D.
These correlation functions can be expressed in terms of amplitudes and decaying exponentials by inserting a complete basis of states. Projecting onto zero momentum and setting q = (M B − M D * , 0, 0, 0) this gives where Note that we have included contributions from oppositeparity states, which depend on imaginary time like (−1) t and arise from using staggered quarks [38], by introducing the sum over a and b. When either a or b is nonzero the corresponding term in the sum is multiplied by a sign factor which oscillates between 1 and −1 in time. We are only interested in the terms with a = b = 0 here; however in order to extract these, the oscillating terms must be fit away. For our choice of operators the A, B and V parameters are real [52]. We discuss our fits to these correlation functions in Sec. V A.

IV. ONE-LOOP MATCHING
We require a lattice current with the same matrix elements as the continuum current to a given order. The matching of lattice and continuum currents is done in [53] through O(α s , α s /am b , Λ QCD /m b ),where Λ QCD is a typical QCD scale of a few hundred MeV, following the method used in [54]. Using power counting in powers of Λ QCD /m b a set of lattice currents is selected. At the order to which we work in this paper only the following currents contribute It is convenient for us to also compute the matrix elements of operators entering at O(α s Λ QCD /m b ) This allows for a configuration-by-configuration check of the code: namely that at zero recoil, the three-point correlation functions satisfy the relation C 3ptJ (1) +C 3ptJ (2) − 2C 3ptJ (3) = 0. This identity is derived using integration by parts and the fact that The full matching is a double expansion in Λ QCD /m b and in α s . The matched current is given by where Z is a multiplicative factor from the tree-level massive-HISQ wave function renormalization for the HISQ c quark. The one-loop coefficients η and τ respectively account for the renormalization of J  Table IV [53].
Matrix elements of currents of order α n s Λ QCD /m b vanish to all orders in α s according to Luke's theorem [55]. We will denote by V the matrix elements of the currents J latt divided by meson mass factors, as in (13) with a = b = 0 and n = m = 0. Luke's theorem implies the combination which represents the the physical, sub-leading matrix element, should be very small, only different from zero due to systematic uncertainties.

V. ANALYSIS OF NUMERICAL DATA
In this section we discuss the two main aspects of numerical analysis. First we present fits to the correlation  (16), calculated at lattice quark masses appropriate to each of our gauge-field ensembles [56]. We also give values on each ensemble for the strong coupling constant in the V scheme at a scale of 2/a (from results in [53]). functions, allowing us to determine h A1 (1) on each of the 8 ensembles. Second, we discuss how we infer a physical value for h A1 (1) with an error estimate for uncertainties associated with current matching, discretization, and dependence on quark masses.

A. Fits to correlation functions
We fit the three correlation functions defined in (12) simultaneously using the corrfitter package developed by Lepage [57,58]. This minimises (18) with respect to p, where ∆C(t, p) = C(t) − C TH (t, p) and p i is the ith parameter in the theory, p i prior is its prior value with error σ p i prior . The correlation matrix σ t,t includes all correlations between data points. Fitting correlators from all smeared sources and sinks simultaneously requires the use of as SVD cut on the eigenvalues when determining the inverse of σ 2 . We also exclude points close in time to the source and sink to suppress excited state contributions and speed up the fit.
We look at the effectiveness of the various smearings by fitting each smearing diagonal, i.e. equal radii, set of two and three point correlator functions independently and comparing the result to the full fit. Figure 1 shows an example of this; plots for the full data set appear in Fig. 9 in Appendix E. In these plots we only include the results of fits with χ 2 /dof < 1.2. We give the ground state and oscillating state two point fit parameters for our full simultaneous fits in Table V. The C B2pt (t) fit amplitudes, the energies and B n ai parameters of (12), are in good agreement with those in Ref. [59]. Table VI gives results for matrix elements corresponding to the currents J (1) latt and J (2) latt . One can see that Luke's theorem holds, in that V (1) sub is very small. Results are also V: Ground state and oscillating state local amplitudes and masses from our fits. Note that on Set 3 and for all the D * s data we use the local vector operator, otherwise we use the point-split operator; therefore, the amplitudes A are not comparable between different operators. Also note that the tabulated B "masses" are the NRQCD "simulation energies" aE sim , representing the nonperturbative contribution to the B meson binding energy. The B parameters are in good agreement with those in [59].  (13)    Note the approximate cancellation between the mixing down term αsτ V (0) and V (1) to give a small V sub as we would expect from Luke's theorem. Note V (2) is numerically smaller than its parametric estimate αsΛQCD/m b ≈ 0.03.
On each ensemble, we obtain a value for the zero-recoil form factors h (s) A1 (1). As in the continuum expressions (3) and (5) we have and similarly for h s A1 (1)| latt . We write V J here to make clear that we fit combinations of three point correlators that correspond to the insertion of the current given by (16). Results for h A1 (1) on each ensemble are presented in Table VII. We computed h s A1 (1) on the physical-point lattices only since chiral perturbation theory predicts this quantity to be much less sensitive to the sea quark mass than the spectator quark mass. (In fact we will see that the spectator quark mass dependence is also small).

B. Chiral-continuum extrapolation
By carrying out the calculation using 8 ensembles, spanning 3 values of lattice spacing and 3 values of the light quark mass, we can quantify many of the systematic uncertainties by performing a least-squares fit to a function which accounts for unphysical parameters or truncation errors. Below we describe how the fits address each of these sources of uncertainty then present results of the fits.
There are two types of systematic error for which we must account. The first type are truncation errors about  (44) which the numerical data contain no information. In this class are the higher-order (in Λ QCD /m b ) current corrections truncated in the perturbative matching described in Sec. IV. The numerical data contain no information about Λ 2 QCD /m 2 b or α s Λ QCD /m b corrections, so we add to each data point nuisance terms and e 4 , e 5 , e 6 , e 7 , e 8 , and e 9 are Gaussian distributed variables, with mean and standard deviation µ(σ), with e 4 = 0(0.5), e 7 = 0(0.3) and e 5,6,8,9 = 0(1), 100% correlated between each data point. The e 5 , e 6 , e 8 and e 9 terms reflect the fact that the coefficients of the truncated Λ 2 QCD /m 2 b and α s Λ QCD /m b terms will be slowly varying functions of am b . Our choice of e 7 is motivated by the magnitude of V (2) and the expectation that Luke's theorem will hold at this order.
The second type of systematic uncertainties arise from truncation, discretization, or tuning errors about which we can draw inferences from our Monte Carlo calculation. Consider the unknown α 2 s corrections to the current normalization. In contrast to the truncation of the Λ QCD /m b expansion, the numerical data is, at least in principle, sensitive to O(α 2 s ) corrections through the running of the coupling on the different lattice spacings. In addition the results have dependence on the lattice spacing and the light quark mass that can be mapped out using theoretical expectations. For the light quark mass dependence this is based on chiral perturbation theory. Therefore we fit the data points to the functional form The first term accounts for the deviation of the physical h A1 (1) from the static quark limit value of 1. The fit parameter B is given a prior of 0(1). We take as priors γ 1 = 0(0.5), γ 5,6 = 0(1). Discretization and quark mass tuning errors are included in δ B a , to be described further below.
The second and third terms in (21) give the leading dependence on the light quark mass, parametrized by M 2 π divided by the chiral scale Λ χ , which we set to be 1 GeV. The coefficient of the chiral logs depends on the D * Dπ coupling g, which we take as 0.53 (8) following [18], and on the pion decay constant in the physical pion mass limit f = 130 MeV. The D * − D mass splitting, ∆ mc , appearing in the chiral logs is taken as 142 MeV. The uncertainties from f and ∆ mc are negligible compared to the error on g and are not included. Further details about the staggered chiral perturbation theory [60] input to (21) are given in Appendix F. We will return to discuss δ g a shortly.
The fourth term in (21) is present in the fit since the current matching has truncation errors of O(α 2 s ). The truncated term would have some mild dependence on am b , which is reflected in the ansatz for this term.
The δ B a and δ g a in (21) parametrize how discretization and quark mass tuning errors could enter the fit form. These originate from the gauge action, the NRQCD action and the HISQ action. In all three actions discretization errors appear as even powers of a, hence we include multiplicative factors Each factor b i contains a distinct sea quark tuning error dependence where the κ i are given a Gaussian prior 0(0.5). Note that we do not include a κ 0 term for the O(a 0 ) piece as such a term would not represent a mistuning error or discretization effect. The product on the right-hand side allows for effects of small mistunings in the sea quark masses and the valence charm and bottom quark masses. For the sea u/d and s quarks we include a multiplicative factor where m sea = 2m l + m s and δx sea = m sea − m phys sea . The physical masses are taken from [61] and are computed using the η s mass. We take m phys l /m phys s = 27.4 [48]. We also include the multiplicative factor where δm c = m c − m phys c , with physical mass taken from [47], and the factor is determined from the spin-averaged kinetic mass of the Υ and η b [46]. c i , d i , and f i are given prior values of 0(0.5). We neglect the effects of the very small mistuning of the light quark masses from their physical value which we expect to be small.
Finite volume corrections to the staggered chiral perturbation theory are given in [60]. Evaluating these expressions on our lattices, we have found that finite volume effects are at least an order of magnitude smaller than the leading O(α 2 s ) error on the unphysical lattices. On sets 3, 6 and 8 the finite volume effects are larger, around half a percent in size. This is significant at the order to which we work. To account for these effects we subtract the finite volume correction to h A1 (1) from our data for these ensembles. We further discuss finite volume effects in Appendix F.
The calculation on each ensemble of the form factor for B s → D * s decay is equivalent to the B → D * calculation, with the light quark propagator replaced with a strange quark propagator. The analysis is substantially more straightforward, both because the data is less noisy and because no chiral extrapolation is required. Before fitting the lattice data, we include a term to account for the absence of O(Λ 2 QCD /m 2 b ) and O(α s Λ QCD /m b ) effects, as in (20), using the same Gaussian variables e 4 , e 5 , e 6 , e 7 , e 8 , and e 9 .
For the continuum-chiral fit to the h s A1 (1) we take the functional form to be the following, where δ sB a B s has the same form and priors as the term included for the where γ 1 , γ 5 and γ 6 are the same as in (21) because these terms represent the same higher order matching corrections. We run the B s → D * s fit simultaneously with the B → D * fit.
The NRQCD and HISQ systematics are the same as before, and we expect negligible isospin breaking and finite volume effects. In Figure 2 we show the M 2 π dependence of our B → D * data and the extrapolated continuum chiral form.
We present results for the h A1 (1) and h s A1 (1) fit param- Table VIII. Plots showing the a 2 dependence of our B → D * and B s → D * s data are shown in Figures 3 and 4 respectively, together with the result of our fit. The O(a 4 ) and O(a 6 ) parameters default to their prior values, while the O(a 2 ) parameters are consistent with zero. We tried various modifications to our fit, the results of which we present in Appendix F. Table IX presents a summary and combination of the uncertainties in our results for h A1 (1) and h s A1 (1).

C. Isospin breaking effects
The effects of electromagnetic interactions and m u = m d on h A1 (1) are negligible compared to the dominant uncertainties quoted in Table IX. We find only a variation of 0.25% in the chiral-continuum fits to h A1 (1) whether the π 0 or π + mass is used as the input value for the physical limit. Electroweak and Coulomb effects in the decay rate (1) are presently accounted for at leading order by a single multiplicative factorη EW to be discussed below in Sec. VII. As lattice QCD uncertainties are reduced in the future, it will be desirable to more directly calculate the effects of electromagnetism in a lattice QCD+QED calculation, where m u = m d can also be implemented.

VI. RESULTS AND DISCUSSION
We have calculated the zero recoil form factor for B → D * ν decay using the most physically realistic gluon field configurations currently available along with quark discretizations that are highly improved. Our final result for the form factor, including all sources of uncertainty, is It is clear from this treatment that the dominant source of uncertainty is the O(α 2 s ) uncertainty coming from the perturbative matching calculation. In principle this could be reduced by a two-loop matching calculation; however, such calculations in lattice NRQCD have not been done before. It is worth noting that for our calculation this uncertainty is somewhat constrained by the fit, as is reflected in Table IX. It has also been suggested [62] that it could be estimated using heavy-HISQ b quarks on 'ultrafine' lattices with a = 0.045 fm and m b a < 1. There we can use the nonperturbative PCAC relation and the absolute normalization of the pseudoscalar current to normalise J (0) , using (m b + m c )P = Z∂ µÂ µ to find the matching coefficient Z and then comparing matrix elements of this normalized current to the result using perturbation theory.
Within errors, our result agrees with the result from the Fermilab Lattice and MILC Collaborations [18], h A1 (1) = 0.906(4) (12). The higher precision achieved in this work is due to the use of the same lattice discretization for the b and c quarks. This enabled them to avoid the larger current-matching uncertainties present in our NRQCD-b, HISQ-c work. Nevertheless, the value of providing a completely independent lattice QCD result using different formalisms is self-evident.
After combining the statistical and systematic errors in quadrature, a weighted average of the two lattice results yields h A1 (1) = 0.904 (12). We use this value in our discussion in Sec. VII.
Our result for the B s → D * s zero-recoil form factor is This is the first lattice QCD calculations of this quantity. We see no significant difference between the result for B → D * and B s → D * s showing that spectator quark mass effects are very small. Correlated systematic uncertainties cancel in the ratio, which we find to be We find there to be no significant U -spin (d ↔ s) breaking effect at the few percent level.

VII. IMPLICATIONS FOR |V cb |
Until recently, one would simply combine a world average of lattice data for h A1 (1) with the latest HFLAV result for theB 0 → D * + − ν differential branching fraction extrapolated to zero recoil:η EW F(1)|V cb | = 35.61(11)(44) × 10 −3 [17]. Doing so with the weighted average of the Fermilab/MILC result and ours yields |V cb | HFLAV = (38.9 ± 0.7) × 10 −3 , where we have used the estimated charge-averaged value ofη EW = 1.015(5) [18]. The uncertainty in |V cb | HFLAV is due in equal parts to lattice and experimental error. Recent work analyzing unfolded Belle data [16] has called into question the accuracy of what has become the standard method of extrapolating experimental data to zero recoil [22][23][24][25][26][27]. In order to understand our new result for h A1 (1), as well as to prepare for future lattice calculations and experimental measurements, we carry out a similar analysis here. We generally agree with conclusions already in the literature, but we present a few of our own suggestions for how one could proceed in the future.
The method used by experiments to date is due to Caprini, Lellouch, and Neubert (CLN) [21]. Their paramatrization of the form factors entering the differential decay rate and angular observables is an expansion about . A full accounting of the breakdown of systematic errors is made difficult by the fact that smaller priors not well constrained by the data are mixed in a correlated way by the fitter; these are reflected in the total systematic uncertainty. Note that the uncertainty from missing α 2 s terms in the matching for hA 1 (1) and h s A 1 (1) is constrained somewhat by the fit; a naive estimate would give 3.5% on the fine lattices.
Uncertainty In the case of the h A1 (w) form factor it was found that the kinematic variable z gives a more convergent series. Given a specific choice of t 0 , z depends on the t = q 2 as with t ± = (M B ± M D * ) 2 . Usually one takes t 0 = t − , and this is the choice assumed throughout this paper. 1 1 One can express z(t, t − ) as a function of w as The CLN form factors are given as follows with the coefficients computed to be [21] r h2r = 53 , r h2 = −15 , r h3r = −231 , r h3 = 91 , r 11 = −0.12 , r 12 = 0.05 , r 21 = 0.11 , These numbers are the result of a calculation in HQET, using QCD sum rules and neglecting contributions of α s Λ QCD /m c and (Λ QCD /m c ) 2 , as well as smaller effects. Until recently effects of neglecting these terms have not been included in fitting the experimental data.
Ref. [21] claims an accuracy of 2%; however this is based on comparing an expansions in z against some full expressions. While this tests the convergence of the expansions, it does not test the accuracy of numerical factors computed in truncated HQET. In fact the data do not require any higher order terms in z or w−1. We found no effect when including a z 4 term or (w − 1) 3 terms in (33) with Gaussian priors allowing the coefficient r h4 to be up to O(10 3 ) and r 13 , r 23 to be up to O(1).
Nevertheless none of this accounts for higher order terms in the HQET. We can get some idea of how the fit is affected by allowing the r coefficients (34) to be fit parameters with Gaussian priors, with means equal to the CLN values but with widths which we vary. Table X shows the results of fitting to the CLN parametrization. We present six variations, which we describe below. In order to infer |V cb | from the lattice h A1 (1) and the fit to data, the main output is the combination In the first fit, we treat the r-coefficients (34) as pure numbers; this has been the standard treatment until recently. The value of I we obtain agrees with the unfolded fit result of Belle [16], I = 34.9(1.5). X: Fits to the unfolded Belle data using the CLN parametrization. The first fit does not account for any uncertainties in the r coefficients (34). The next three include the r coefficients as Gaussian priors with widths of 10%, 20% or 100% uncertainties, respectively. The final two fits assign 10% or 20% uncertainty to the coefficients in hA 1 (w) and allow the coefficients of R1(w) and R2(w) to be O(1). It would be better to include estimates of HQET truncation errors in the fits. We implement this by treating the r-coefficients as fit parameters, adding Gaussian priors with central values as in (34) and with widths equal to our uncertainty. Unfortunately it is not clear how accurately these are known at this order in HQET. We note both α s Λ QCD /m c and (Λ QCD /m c ) 2 are roughly 0.1, so one approach is to suggest truncated terms could vary each of the r's by 10%. However, some linear algebra has been done after truncating HQET expressions to arrive at the form factors (33). This could enhance (or suppress) the truncation error in some terms, and the opposite in others. The fit does not change much if the uncertainties are 20%, but 100% uncertainties in (34) do affect the fit result. Most notably, the value of I increases by 5%, i.e. one standard deviation.
The smallness of the coefficients in the expansions of R 1 and R 2 is likely due to cancellations in the expansions when ratios are taken. Therefore, assuming a relative error on the r ij (i, j = 1, 2) is probably not correct. We present two fits where these coefficients are given Gaussian priors equal to 0 ± 1, while the coefficients in the expansion of h A1 (w) are given 10% or 20% uncertainties. The resulting values for I lie in between the tightly constrained fits and the 100% uncertainty fit.
The fact that the tightly constrained CLN fits describe the data well, with good χ 2 for example, is a success for HQET. It shows that the important physics has been captured within the accuracy of the theory. However, now that we are in the high precision era of flavour physics, we ought to be wary about the accuracy of the assumptions which go into fitting the data. The observation that I increases under a relaxation of assumptions about the r-coefficients agrees with other authors' findings [22][23][24][25][26][27].
An alternative parametrization for the hadronic form  [65,66] factors is the one proposed by Boyd, Grinstein, and Lebed (BGL) [28]. In their conventions the three form factors entering (assuming the lepton mass can be neglected) are f (q 2 ), F 1 (q 2 ), and g(q 2 ). Two of the form factors are kinematically constrained at q 2 = 0: . Each of these is expanded in a Taylor series about z = 0 after factoring out a function intended to account for nearby resonances. Abbreviating t = q 2 , form factors are parametrized by Throughout this paper we take t 0 = t − . With appropriately chosen Q F , the magnitudes of the coefficients a (F ) n are bounded by unitarity constraints.
Even stronger bounds can be imposed if one is able to include all the B ( * ) → D ( * ) matrix elements, with (pseudo)scalar and (axial)vector initial and final states [25], but this is outside the scope of our analysis here. The two functions in (37) are the outer functions φ F (z), which can be found in the literature e.g. in Refs. [23,24,28], and the Blaschke factor The product is over a set of applicable resonances, the vector B * c states for g and the axial vector B c states for f and F 1 . The resonances included in the Blaschke factor should be those with the appropriate quantum numbers and below scattering threshold. There are 4 B c vector and 4 axial-vector states conjectured to be below BD *     threshold. Table XI lists calculations of the vector and axial vector B c resonances. The model estimate for the mass of the heaviest vector state is very close to threshold, so has been left out of several analyses, including here. The magnitude of the Blaschke factors can be very sensitive to n, so leaving states out reduces the strength of unitarity constraints. This is illustrated in Fig. 5 for the Q F (q 2 ) for F = f , F 1 , and g. Table XII shows the results of BGL fits to the unfolded Belle data [16], varying the number of states included in the Blaschke factor and the number of terms kept in the z-expansion. The fits enforce the q 2 = 0 constraint on F 1 (0) and f (0) at the 1% level. Priors for the coefficients a (F ) k are Gaussians with mean 0 and standard deviation 1. Only the k = 0 and 1 coefficients are tabulated; the others are not constrained by the data and remain statistically consistent with 0. As discussed above, the magnitude of these coefficients depends on the number of states in the Blaschke factor. Nevertheless, the results for I are insensitive to this. On the other hand, I does increase by about 0.001, or approximately 0.7σ, when switching from K = 2 to higher order polynomials in z. (Results remain the same for K > 4.) The fits presented in Table XII do not enforce the unitarity bounds (38), but these bounds are not close to being saturated unless only two resonances are included in the Blaschke factors. Performing the fits with the bounds enforced did not significantly affect results. Considering the n B > 2, K = 4 fits, increasing the standard deviation of the Gaussian priors for the series coefficients a 1 , S f F , S g , increased by factors of 1.5-2. For most of the a (F ) k with k ≥ 2, the posterior distribution is the same as the prior, the exception being that a (F1) 2 < ∼ 0.5 seems preferred by the fit, even with wide prior widths.
To the extent that unitarity constraints do not affect the BGL fits, then a simpler approach would be to represent Q F (t) by a simple pole, as in the simplified series expansion (BCL) [32]. That is, one can parametrize the form factors by (36) with with M P being the mass of the lightest resonance with the appropriate quantum numbers. The normalization N F can be chosen so that the series coefficients are of the same order of magnitude as in a particular BGL expansion. With Fig. 5 as a guide, we take N f = 300, N F1 = 7000, and N g = 5. Once again we fit with priors for a (F ) k equal to 0 ± 1. The results for I show the same behaviour for the BCL fits as for the BGL fits.
The virtue of the BCL fit is in its simplicity. The BGL fit requires theory input for the outer functions φ F : perturbatively calculated derivatives of two-point functions at q 2 = 0 and n I , the number of spectator quarks adjusted to account for SU (3) F breaking. (In the BGL fits here we take the values given in Table 2 of Ref. [23].) The Blaschke factor requires as input model estimates for the excited B c resonances to include in the Blaschke factor. If unitarity bounds become tight enough to have an effect on the fits to data, then the effects of theoretical assumptions needs to be carefully included in the error analysis. On the other hand, the BCL fits only take as  (36), either using the BGL (37) or BCL (40) parametrization. Unitarity constraints are not enforced in the fit, but the sums Sg and S f F (38) are given for reference (see text). The number of 1 + /1 − resonances included in the Blaschke factor is n + B /n − B . Terms up to O(z K−1 ) are included in the fits. Coefficients of higher order terms are consistent with zero. additional input the mass of a single resonance, available to very good precision from lattice QCD. In the future, fits to the BCL simplified z-expansion could provide a clean, benchmark fit. Fig. 6 summarizes the consequences to I of different fitting choices selected from Tables X and XII. The top two points show results from CLN fits including no uncertainties on the coefficients (34), or 10% errors on the r h coefficients and allowing the coefficients in the expansions of R 1,2 (w) to be 0 ± 1. The bottom two points are respectively BGL and BCL fits with K = 4, and n + B = 4, n − B = 3 for the BGL fit. In Fig. 7 we compare the fit results, integrated over the experimental bins, of the tightly constrained CLN fit and the BGL and BCL fits (with K = 4) to the Belle data [16]. The agreement is generally good, with the notable exception of the dΓ/dw in the smallest w bin, where the CLN result is in greater tension with the data than the BGL and BCL results.
For the time being, with only one experimental data set available to carry out these investigations, determinations of |V cb | from B → D * ν are less certain than has been thought. The BGL and BCL fits to Belle data indicate I = 0.038(2). Ref. [18] cites a private communication with C. Schwanda givingη EW = η EW η Coulomb = 1.0182 (16) as the product of the electroweak factor η EW = 1.0066(16) and a term accounting for electromagnetic interactions between the charged D * and lepton in the final state. Combining this with the weighted average for h A1 (1) from Fermilab/MILC [18] and this work, we arrive at where the error is dominated by the experimental and related fitting uncertainty. This determination agrees well with both those from inclusive and exclusive B → D ν decays as shown in Fig. 8.
One may ultimately obtain a more precise determination of |V cb | by including all relevant information, from HQET, by imposing stronger unitarity bounds [25], and including light cone sum rule calculations of form factors at large recoil [68]. Comparison of the different approaches would be helpful to highlight the impact of including different ingredients.  [16]. The binned fit results are slightly offset from the bin midpoints for clarity. See Appendix G and Ref. [16] for definitions.  [19,20] and B → D ν [33].

VIII. CONCLUSIONS
We present new unquenched lattice QCD determinations of the zero-recoil form factors h A1 (1) and h s A1 (1), sometimes denoted F B→D * (1) and F Bs→D * s (1), respectively. We have used 8 ensembles spanning 3 lattice spacings and 3 values of light-to-strange quark mass ratios, including the physical point. Our results are This result for h A1 (1) provides a valuable, independent check of the Fermilab/MILC result [18]. We have used completely independent sets of gauge field configurations and different formulations for the charm and bottom quarks. The two results are in good agreement.
While the determination of |V cb | using these results is complicated by the need to investigate assumptions used in extrapolating experimental data to zero recoil, series expansion fits to the unfolded Belle data yield This is consistent with recent determinations using exclusive B → D ν and inclusive decays (Fig 8).
A reanalysis of BaBar data for the differential decay rate would complement the unfolded Belle data used here. We can also look forward to new data from Belle II, after which the the precision of |V cb | from B → D * ν is likely to be much improved. Lattice QCD data away from zero recoil will also help reduce the uncertainties. Preliminary results from the Fermilab/MILC collaboration were presented at the Lattice 2017 conference [69].
Our result for the B s → D * s form factor is the first complete calculation of h s A1 (1). In the future, measurements of the exclusive decays with a strange spectator, B s → D ( * ) s ν, could also provide a constraint on |V cb |. LHCb has reconstructed B 0 s → D * − s µ + ν µ decays [70]. Eventually, with properly normalized branching fractions, these will also provide a method of constraining |V cb |.
Spectator quark mass effects are bounded by our calculation of the ratio h s A1 (1)/h A1 (1) and its consistency with unity. We find deviations from d ↔ s symmetry in the zero recoil B (s) → D * (s) form factors to be no more than 2-3%.

ACKNOWLEDGMENTS
We thank R. R. Horgan for useful discussions and P. Gambino for correspondence regarding fits to experimen-tal data. We are grateful to the MILC collaboration for making publicly available their gauge configurations and their code MILC-7.7.11 [50]. This work was funded in part by STFC grants ST/L000385/1 and ST/L000466/1. The gauge action used to generate the configurations is the Symanzik and tadpole improved action of [49], which contains additional rectangle and parallelogram loops to cancel radiative O(a 2 ) errors: Where −µ indicates a Hermitian conjugated gauge link. a 1 and a 2 are calculated in terms of a 0 using lattice perturbation theory. The perturbative coefficients are specified in Ref. [49].
For the u/d and c valence quarks in our calculation we use the same HISQ action as for the sea quarks [38]. The advantage of using HISQ is that am q discretization errors are under sufficient control that it can be used both for light and for c quarks [38,73,74]. The HISQ action is also numerically inexpensive as a result of the staggering which means we are able to attain better statistics. The valence u/d masses are the same as those in the sea. The masses are given in Table I. Below we summarize a few relevant facts.
The naive Dirac action has a discrete, space time dependent symmetry and following [38] The conventions for γξ are specified in the appendices.
In momentum space this then gives the relation for the naive quark propagator: One can diagonalise the naive action in spin indices using a position dependent transformation of the fields. There are several choices for such a transformation, here we use: with Ω(x) = γ x this yields the action with propagator We then need only do the inversion for a single component of χ and the full naive propagator can be reconstructed trivially by inserting Ω matrices: In order to remove discretization errors and taste exchange violations the operator ∆ µ (U ) used in simulations is more elaborate. It retains the feature that ∆ µ (U )ψ(x) only contains fields ψ(x ) located an odd number of lattice sites away from x in the µ direction, ensuring that the spin-diagonalization (12) still works. The full, highly improved staggered SU (3)-covariant derivative operator is [38]: Where 'symm' indicates that the product ordering is symmetrised in ρ, δ ρ approximates a covariant first derivative on the gauge links and δ (2) ρ approximates a second covariant derivative.
The third covariant derivative term removes order a 3 discretization errors coming from the approximation of the derivative. Without the epsilon term, tree level discretization errors appear going as (ap µ ) 4 . For the mesons we are interested in quarks are typically nonrelativistic, and so the error is dominated by the energy, and ultimately the mass contribution going as (am) 4 . For light quarks this is negligible, but for charm physics this must be included since current lattice spacings have am c ≈ 0.5. can be calculated straightforwardly as an expansion in (am) 2 by requiring the tree level dispersion relation lim p→0 (E 2 (p) − m 2 )/p 2 to have its correct value, 1, to a given order O(am). The expansion is [ 8 ) .
(C12) The smearings F µ remove taste changing interactions, since δ (2) ρ ≈ −4/a 2 when applied to a link carrying mo-mentum q ρ ≈ π/a The µ direction needn't be smeared as the original interaction vanishes in this case anyway. The smearing F µ , known as "Fat7" smearing [75], introduces new O(a 2 ) errors. These are removed by replacing F µ with [76] Where F ASQTAD µ is the gauge link smearing employed in the widely used a-squared tadpole improved action. Similar errors originating from the smearing on the third derivative term needn't be corrected as they go as O(a 4 ). A single smearing introduces perpendicular gauge links which are themselves unsmeared. To further suppress taste exchange we use multiple smearings. Once such smearing is: where U is a reunitarization. This combination ensures that each smearing does not introduce any additional O(a 2 ) errors, and ensures no growth in the size of two gluon vertices, since the unitarization ensures it is bounded by unity. In the HISQ operator defined in (17) we have moved the entirety of the O(a 2 ) corrections to the outermost smearing.
In order to check the taste exchange violations in HISQ one can check for taste-splittings of the pion masses. However since there are more allowed effective taste exchange vertices that there are degenerate pion multiplets this does not guarantee the theory is free of taste exchange. A better check is the explicit calculation of the couplings to taste exchange currents required to remove taste exchange. These are given in [38] in which it is clear that the HISQ action is a significant improvement over the older ASQTAD action.
where it is understood that when we add δ st it is modulo the hypercube. We have used the noise condition: to insert the random walls. Setting Now, we do not have S c ab (x, y), we have S c ab (y, x) so we can use: is the local spin-taste phase. Inserting Dirac indices: We recognise S c bc (y, x)(−1) x β M (x)ξ ca (x − σ 1 + δ st )∆ 1 (σ 1 ) as the MILC KS propagator. The naive active quark that gets made in NRQCD is then: and the contractions to do are Current ab,αβ (y) = Active * ba,κα (y)Γ κβ C 3pt = Current ab,αβ (y)Ext ba,βα (y) . (D9) Appendix E: Correlator Fits Fig. 9 shows comparison of the fit results for h A1 (1) when varying numbers of exponentials; the points are normalized by the value of taken h A1 (1) as our result for that ensemble. Plots are shown for all 8 ensembles as listed in Tab. I. In each plot, we show the full fit results to the 3 × 3 matrix of source/sink combinations (local l, or Gaussian with 2 radii, g2 and g4), as well as "diagonal" fits where only one source/sink is used. The statistical improvement of using all the data is apparent. The flatness of the curves and the constancy of the error bars shows that, for large enough N exp , the Bayesian fits are insensitive to adding further exponential terms, i.e. effects of excited states are accounted for. Our final results typically come from the N exp = 5 fits to the full 3 × 3 matrix of correlators; however, on ensembles 3 and 7, we had to include another exponential.  In each plot 4 sets of data points are shown: the full fit including all 3 × 3 source-sink combinations, and, for comparison, separate "diagonal" fits where only one type of source-sink smearing is used. (The notation is defined in Sec. III.) A significant improvement is seen in the full fit. All diagonal fits show good agreement for Nexp ≥ 4, but with the increased precision, sometimes 5 or 6 exponentials are needed to get a good 3 × 3 matrix fit.
whereF X = F [m X , −∆ mc /m X ] and The masses of the η and η are given in [77] as We take the ss pseudoscalar taste splittings equal to the pion taste splittings. This is a good approximation in the case of HISQ [36]. We can then write (to order O((a 2 δ V ) 2 ) ) from which we find The expression for h A1 (1) then reduces to where we have ignored terms expected to produce normal discretization errors and pion mass dependence, as these are included elsewhere in the fit. Following [78] we take δA ≈ δV ≈ −δt, which we implement by including δA = δV = −δt × 1.0(5) as priors. We use the pion masses computed in [78] together with the taste splittings for the pion, δt, given in [36].
Finite volume effects can be accounted for in heavy meson chiral perturbation theory [79] including tastesplitting effects in the staggered pions [60]. The functionsF X in (F6) receive a correction term corresponding to the difference between infinite volume loop integrals and finite volume discrete sums. Taste-splitting effects in the pions at non-zero lattice spacing moderate the size of the finite volume corrections because some of the pions in the loops have heavier masses than the Goldstone pion. Consequently, some of the finite volume effect appears as a lattice-spacing effect, which is dealt with by our chiralcontinuum fit.
We incorporated the finite volume corrections into our fit by subtracting from our data δ F V h A1 (1), found by adding δF X to eachF X appearing in (F6). In Fig. 10 we show δ F V h A1 (1) as a function of pion mass for the parameters appropriate for the physical pion mass lattices, Sets 3, 6, and 8 (see Table I). For the other lattices, |δ F V h A1 (1)| ≈ O(0.1%) over the M π range where we have data and is not significant.
In Table XIII we give fit results for plausible variations on our chosen fit function as a demonstration of stability under such nontrivial choices. Neglecting different powers of a 2 we see that our result is only sensitive to leading O(a 2 ) errors. The M 2 π dependence we included does not affect the central value if removed, nor do changes in the assumed correlations between NRQCD systematics between ensembles. Removing taste splitting terms in the chiral perturbation theory result down to the continuum formula results in only a small change to the central value. Adding α s Λ QCD /M B , which we have excluded from our fit due to Luke's theorem, results in a slight increase in the central value as well as the expected increase in error. Our result is also only mildly sensitive to different choices of Λ QCD which we vary by ±50%. Taken collectively we note that no tested variations result in more than a 0.25σ change to the central value.