B s → D s l ν form factors for the full q 2 range from lattice QCD with nonperturbatively normalized currents

We present a lattice QCD determination of the B s → D s l ν scalar and vector form factors over the full physical range of momentum transfer. The result is derived from correlation functions computed using the highly improved staggered quark (HISQ) formalism, on the second generation MILC gluon ensembles accounting for up, down, strange and charm contributions from the sea. We calculate correlation functions for three lattice spacing values and an array of unphysically light b -quark masses, and extrapolate to the physical value. Using the HISQ formalism for all quarks means that the lattice current coupling to the W can be renormalized nonperturbatively, giving a result free from perturbative matching errors for the first time. Our results are in agreement with, and more accurate than, previous determinations of these form factors. From the form factors we also determine the ratio of branching fractions that is sensitive to violation of lepton universality: R ð D s Þ ¼ B ð B s → D s τν τ Þ = B ð B s → D s l ν l Þ , where l is an electron or a muon. We find R ð D s Þ ¼ 0 . 2993 ð 46 Þ , which is also more accurate than previous lattice QCD results. Combined with a future measurement of R ð D s Þ , this could supply a new test of the Standard Model. We also compare the dependence on heavy quark mass of our form factors to expectations from heavy quark effective theory.


I. INTRODUCTION
The weak decay processes of mesons such as the B and B s , containing b quarks, are a key potential source of insights into physics beyond the Standard Model (SM). Flavor-changing B decays have gained a lot of interest because of a number of related tensions between experimental measurements and SM predictions [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. These tensions drive the need for improved theoretical calculations in the SM using methods and studying processes where we have good control of the uncertainties.
Lattice QCD is the method of choice for providing the hadronic input known as form factors that determine, up to a normalization factor, the differential branching fraction for exclusive decay processes such B → Dlν (we suppress all electric charge and particle-antiparticle labels here in referring to decay processes). The normalization factor that can then be extracted by comparison of theory with experiment is the Cabibbo-Kobayashi-Maskawa matrix element, in this case jV cb j [20][21][22][23][24]. Determination of jV cb j then feeds into constraints on new physics through, for example, tests of the unitarity triangle.
There has been a longstanding tension in determinations of jV cb j between exclusive (from B → Dlν and B → D Ã lν decays), and inclusive (from B → X c lν, where X c is any charmed hadronic state) processes. The most accurate exclusive results came from studies of the B → D Ã lν decay at zero recoil. It now seems likely that the uncertainties there were being underestimated because of the use of a very constrained parametrization in the extrapolation of the experimental B → D Ã data to the zero recoil limit [25][26][27][28], but see also [29,30]. This underlines the importance in future of comparing theory and experiment across the full range of squared 4-momentum transfer (q 2 ) (a point emphasized for D → K in [31]). It also demonstrates the need for comparison of accurate results from multiple decay processes for a more complete picture. Improved methods for producing the theoretical input to jV cb j excl , namely lattice QCD determinations of form factors, are clearly necessary.
Here we provide improved accuracy for the form factors for the B s → D s lν decay using a new lattice QCD method that covers the full q 2 range of the decay for the first time. Preliminary results appeared in [32]. The B s → D s form factors are more attractive than B → D for a first calculation to test methodology. They are numerically faster to compute and have higher statistical accuracy and smaller finite-volume effects because no valence u=d quarks are present. Chiral perturbation theory [33] expects that the B → D form factors should be relatively insensitive to the spectator quark mass and hence should be very similar between B s → D s and B → D. This is confirmed at the 5% level by lattice QCD calculations [34,35]. Hence improved calculations of B s → D s form factors can also offer information on B → D.
Given an experimental determination, the B s → D s lν decay can supply a new method for precisely determining the Cabibbo-Kobayashi-Maskawa (CKM) element jV cb j. It can also supply a new test of the SM through quantities sensitive to lepton universality violation. We give the SM result for RðD s Þ ¼ BðB s → D s τν τ Þ=BðB s → D s lν l Þ, where l ¼ e or μ. An experimental value for comparison to this would help to clarify the tension found between the SM and experiment in the related ratios RðD ðÃÞ Þ ¼ BðB → D ðÃÞ τν τ Þ=BðB → D ðÃÞ lν l Þ [36] (see also a preliminary new analysis by Belle [37]).
Three lattice QCD calculations of B ðsÞ → D ðsÞ form factors have already been performed. The FNAL/MILC collaboration [22,34] used the Fermilab action for the b and c quarks and the asqtad action for the light quarks on MILC gluon field ensembles that include 2 þ 1 flavors of asqtad sea quarks. On the same gluon field ensembles the HPQCD collaboration has calculated the form factors using nonrelativistic QCD (NRQCD) for the valence b and the highly improved staggered quark (HISQ) action for the other valence quarks [23,35]. A further calculation has been done using maximally twisted Wilson quarks on n f ¼ 2 gluon field ensembles [38]. Preliminary results using domainwall quarks are given in [39].
A considerable limitation in the FNAL/MILC and HPQCD/NRQCD studies is the requirement for normalization of the lattice QCD b → c current. The matching between this current and that of continuum QCD is done in lattice QCD perturbation theory through Oðα s Þ, giving a Oðα 2 s Þ systematic error which can be sizeable. Systematic errors coming from the truncation of the nonrelativistic expansion of the current are also a problem. In the Fermilab formalism the missing terms become Oðα s aÞ discretization effects on fine enough lattices; in the NRQCD formalism they mix discretization effects and Oðα s =m b Þ (where m b is the b quark mass) relativistic corrections. Here we dispense with both of these problems by using a relativistic formalism with absolutely normalized lattice QCD currents.
Another limitation present in each of the previous studies is that the lattice QCD results are limited to a region of high q 2 , close to zero recoil. The reason for this is mainly to avoid large statistical errors. The signal/noise degrades exponentially as the spatial momentum of the meson in the final state grows. For b quark decays the maximum spatial momentum of the final state meson can be large (tending to m B =2 for light mesons, where m B is the B meson mass). Systematic errors from missing discretization (and relativistic) corrections also grow away from zero recoil. This is particularly problematic if discretization effects are OðaÞ as above and relatively coarse lattices are used (to reduce numerical cost). Working close to zero recoil means that the lattice results then have to be extrapolated from the high q 2 region into the rest of the physical q 2 range. Here we also overcome this problem by working with a highly improved quark action in which even Oða 2 Þ errors have been eliminated at tree-level [40]. We cover a range of values of the lattice spacing that includes very fine lattices and include results from lighter than physical b quarks and this enables us to cover the full q 2 range in our lattice calculation.
We perform our calculation on the second-generation MILC gluon ensembles [41], including effects from 2 þ 1 þ 1 flavors in the sea using the HISQ action [40]. We also use the HISQ action for all valence quarks. Our calculation employs HPQCD's heavy-HISQ approach. In this we obtain lattice results at a number of unphysically light masses for the b (we refer to this generically as the heavy quark h), reaching the b quark mass on the finest lattices. This allows us to perform a combined fit in m h and lattice spacing that we can evaluate in the continuum limit at m h ¼ m b [and as a function of m h to compare, for example, to expectations from heavy quark effective theory (HQET)]. By using only HISQ quarks, we can normalize all the lattice currents fully nonperturbatively and avoid systematic errors from current matching.
This calculation adds to a growing number of successful demonstrations of the heavy-HISQ approach. The method was developed for determination of the b quark mass, B meson masses and decay constants [42][43][44] and is now also being used by other groups for these calculations [45,46]. A proof-of-principle application of heavy-HISQ to form factors was given for B c → η c and B c → J=ψ in [47,48], covering the full q 2 range for these decays and this work builds on those results. The B s → D Ã s axial form factor at zero recoil was calculated using heavy-HISQ in [49].
This article is structured in the following way: Section II lays out our lattice QCD approach to calculating the form factors and then Sec. III presents our results, along with several consistency checks and our determination of RðD s Þ. We also give curves showing the heavy-quark mass dependence of some features of the form factors that can be compared to HQET. For those simply hoping to use our calculated B s → D s form factors, Appendix A gives the parameters and covariance matrix required to reconstruct them.

A. Form factors
In this section we specify our notation for the form factors and matrix elements. The differential decay rate for B s → D s lν decays is given in the SM by where m l is the mass of the lepton, η EW is the electroweak correction, q 2 ¼ ðp B s − p D s Þ 2 is the momentum transfer and f s 0 ðq 2 Þ, f s þ ðq 2 Þ are the scalar and vector form factors that parametrize the fact that the decay process involves hadrons. We use superscript "s" to denote the strange spectator valence quark. The allowed range of q 2 values if the final states are on-shell is The form factors are determined from matrix elements of the electroweak current between B s and D s states, hD s jðV − AÞ μ jB s i where V μ ¼bγ μ c is the vector component and A μ ¼bγ 5 γ μ c is the axial vector component. In a pseudoscalar-to-pseudoscalar amplitude, only V μ contributes, since hD s jA μ jB s i does not satisfy the parity invariance of QCD. In terms of form factors, the vector current matrix element is given by Analyticity of this matrix element demands that Via the partially conserved vector current relation (PCVC), the form factor f s 0 ðq 2 Þ is also directly related to the matrix element of the scalar current S ¼bc; In our calculation we determine the form factors by computing matrix elements of the temporal vector current V 0 and the scalar current S. The form factors can be extracted from this combination using expressions derived from Eqs. (3) and (5) (once the currents have the correct continuum normalization-see Sec. II D): Our goal is to compute f s 0 ðq 2 Þ and f s þ ðq 2 Þ throughout the range of q 2 values 0 ≤ q 2 ≤ ðM B s − M D s Þ 2 ≡ q 2 max . We extend the range to q 2 ¼ 0 in order to take advantage of the constraint from Eq. (4).

B. Lattice calculation
This calculation closely follows the approach employed in our calculation of the B s → D Ã s axial form factor at zero recoil [49]. Here, however, we must give spatial momentum to the charm quark in the final state so that we can cover the full q 2 range of the decay.  [41,50]. a is the lattice spacing, determined from the Wilson flow parameter [52]. The physical value of w 0 was determined using f π to be 0.1715(9) fm in [53]. This allows the determination of w 0 =a on individual sets of gluon field configurations to be converted into a value of a. These are given in column 3 with the uncertainty quoted being the statistical uncertainty in w 0 =a on that ensemble [54][55][56]. There is an additional uncertainty in a, correlated across all ensembles, from the physical value of w 0 . N x is the spatial extent and N t the temporal extent of the lattice in lattice units. Light (m u ¼ m d ), strange and charm quarks are included in the sea, their masses are given in columns 5-7.

Set
Handle The gluon field configurations used in this calculation were generated by the MILC collaboration [41,50]. The relevant parameters for the specific ensembles we use are given in Table I. The gluon fields are generated using a Symanzik-improved gluon action with coefficients matched to continuum QCD through Oðα s a 2 ; n f α s a 2 Þ [51]. The gluon fields include the effect of 2 þ 1 þ 1 flavors of quarks in the sea (u, d, s, c, where m u0 ¼ m d0 ≡ m l0 ) using the HISQ action [40]. In three of the four ensembles (sets 1, 3 and 4), 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 the form factors here, since we have no valence light quarks. We quantify this small effect by including a fourth ensemble (set 2) with roughly physical m l0 .
We use a number of different masses for the valence heavy quark am val h0 . This allows us to resolve the dependence of the form factors on the heavy quark mass, so that a fit in m h can be performed and the results of the fit evaluated at m h ¼ m b . With a heavy quark mass varying both on a given ensemble and between ensembles, we can resolve both the discretization effects that grow with large (am val h0 ≲ 1) masses and the physical dependence of the continuum form factors on m h . Using unphysically light hquarks also reduces the q 2 range, meaning that we can obtain lattice results across the full range while the statistical noise remains under control.
Staggered quarks have no spin degrees of freedom. Spinparity quantum numbers are accounted for by construction of appropriate fermion bilinears and including an appropriate space-time dependent phase with each operator in the path integral. We categorize these phases according to the standard spin-taste notation, ðγ n ⊗ γ s Þ, where γ n is the spin structure of the operator in the continuum limit, and γ s is the "taste" structure which accounts for the multiple possible copies of the operator constructed from staggered quark fields.
We have designed this calculation to use only local operators (combining fields at the same space-time point and having ðγ n ⊗ γ n Þ spin-taste) for the calculation of the current matrix elements that we require. This is an advantage since point-split operators can lead to noisier correlation functions. The spin-taste operators we use are: scalar ð1 ⊗ 1Þ, pseudoscalar ðγ 5 ⊗ γ 5 Þ, vector ðγ μ ⊗ γ μ Þ, and temporal axial-vector ðγ 0 γ 5 ⊗ γ 0 γ 5 Þ.
We compute a number of correlation functions on the ensembles detailed in Table I. Valence quark masses, momenta and other inputs to the calculation are given in Table II. We use random wall sources to generate all staggered propagators from the source since this gives improved statistical errors [57]. First we compute two-point correlation functions between meson eigenstates of momentum ap, give the s and c valence quark masses in lattice units, which were tuned in [54]. In column 4 we give the heavy quark masses that we used in lattice units. We use a number of heavy quark masses to enable the heavyquark mass dependence to be determined in our fit. Column 5 gives the ratio of heavy and charm valence quark masses, to give physical intuition for the heavy masses we use. The physical value for this is given by 4.528 (54) (for N f ¼ 4) [54]. Column 6 gives the absolute value of the spatial momentum (in lattice units) given to the D s meson using a momentum twist on the charm quark propagator. These values are chosen with the following rationale: when only two values are given, these correspond to the q 2 ¼ 0 and q 2 max points (except on the fine-physical ensemble, where we use the points q 2 max and q 2 max =2); when three values are given, the momenta correspond to q 2 ¼ 0, q 2 ¼ q 2 max =2, and q 2 max ; when four values are given, these are points corresponding to q 2 max , 3q 2 max =4, q 2 max =2, q 2 max =4 and q 2 ¼ 0. We used twisted momenta in the (1,1,1) direction to minimize discretization effects. Column 7 gives the number of gluon field configurations used for that ensemble, n cfg , and the number of different time sources used per configuration to increase statistics, n src . Column 8 gives the temporal separations between source and sink, T, of the three-point correlation functions computed on each ensemble. Columns 9-13 give the t cut values chosen for the corresponding correlator fits (see Sec. II C); η q denotes the correlators for η s , η c and η h mesons and J denotes all 3-point correlators. Column 8 gives the t cut used for both tastes of H s meson that we use. where hi represents a functional integral over all fields, q, q 0 are valence quark fields of the flavors the M meson is charged under, Γ is the spin-taste structure of M and the division by the number of tastes is required to normalize closed loops made from staggered quarks [40]. We compute these correlation functions for all t values, i.e., 0 ≤ t ≤ N t . We compute correlation functions for a heavy-strange pseudoscalar, H s , with spin-taste structure ðγ 5 ⊗ γ 5 Þ, at rest. In terms of staggered quark propagators this takes the form where g q ðx; yÞ is a staggered propagator for flavor q, and the trace is over color. Here x 0 ¼ 0 and y 0 ¼ t, and the sum is over spatial sites labelled x, y. We also compute correlators for a charm-strange pseudoscalar meson D s , with structure ðγ 5 ⊗ γ 5 Þ. For these correlators we need both zero and non-zero spatial momentum. Nonzero spatial momentum is given to the D s by imposing twisted boundary conditions on the gluon fields when computing the charm quark propagators [58]. Then where g θ q ðx; yÞ denotes a propagator with momentum twist θ. We compute these correlation functions using several different twists to produce the range of momenta given in Table II. We design the c propagators to have momentum ap ¼ japjð1; 1; 1Þ, by imposing a twist θ ¼ N x japj=π ffiffi ffi 3 p in each spatial direction. Necessary for extracting the vector current matrix element, we also compute correlation functions for a non-goldstone pseudoscalar heavy-strange mesons at rest, denotedĤ s . This has spin-taste structure ðγ 0 γ 5 ⊗ γ 0 γ 5 Þ.Ĥ s correlators are computed using where we use the notationz μ ¼ P ν≠μ z ν . We also compute correlators for H c mesons, heavycharmed pseudoscalars, using the same form as those for H s , Equation (8). These are used to find H c decay constants, which are useful in some of our continuum and m h fits. In our fits to heavy-quark mass dependence we will use the mass of the heavy-heavy pseudoscalar meson, η h as a physical proxy for the quark mass. To quantify mistuning of the charm and strange quark masses, we also require masses for η c and η s mesons, identical to η h with h replaced c and s quarks respectively. We compute correlators for each of these at rest, using a spin-taste structure ðγ 5 ⊗ γ 5 Þ, taking the same form as those of the H s , Equation (8). Note that all of the η mesons discussed here are artificially forbidden to annihilate in our lattice QCD calculation. We expect this to have negligible effect, for the purposes of this calculation, on the masses of the η c and the η b [54]; the η s is an unphysical meson that can be defined in this limit in a lattice QCD calculation and is convenient for tuning the s quark mass [53,59].
Three-point correlation functions are needed to allow determination of the current matrix elements for B s → D s decay. We require two sets of such correlation functions, one with a scalar and one with a temporal vector current insertion. The first takes the form Sðy; tÞ ¼cðy; tÞhðy; tÞ: ð11Þ In terms of the staggered quark formalism, both the H s source and D s sink are given structure ðγ 5 ⊗ γ 5 Þ, and the current insertion ð1 ⊗ 1Þ. We combine staggered propagators to construct these correlation functions as: where we fix x 0 ¼ 0, y 0 ¼ t and z 0 ¼ T, and once again the charm propagator is given the appropriate twist θ. We compute these correlation functions for all t values within 0 ≤ t ≤ T, using 3 T values to make sure that excited state effects are accounted for. The T values vary with lattice spacing to give approximately the same physical range and always include both even and odd values. The values are given in Table II. The three-point correlation function with temporal vector current insertion is given by V 0 ðy; tÞ ¼cðy; tÞγ 0 hðy; tÞ: ð13Þ This is generated using spin-taste ðγ 0 γ 5 ⊗ γ 0 γ 5 Þ at theĤ s source, ðγ 5 ⊗ γ 5 Þ at the D s sink, and ðγ 0 ⊗ γ 0 Þ at the current insertion. To achieve this we compute The non-goldstoneĤ s is required here to ensure that taste cancels in the correlation function. The difference between theĤ s and the H s , for example in their masses, is generated by taste-exchange discretization effects. In practice it is very small for heavy mesons [40], being suppressed by the heavy meson mass.

C. Analysis of correlation functions
We now describe our simultaneous multiexponential fits to the correlation functions using a standard Bayesian approach [60,61]. The parameters that we wish to determine are ground-state energies, two-point amplitudes and ground-state to ground-state matrix elements. Our correlation functions, however, are contaminated by contributions from excited states. These excited states must be included in our fits so that the systematic error on the ground-state parameters from the presence of the excited states is fully taken into account. Multiexponential fits are then mandatory, guided by Bayesian priors for the parameters, discussed below. To reduce the number of exponentials needed by the fits, we drop values of the correlation functions when they are within t cut of the endpoints (where excited states contribute most). We use values of t cut varying from 2 to 10 throughout the correlator fits; these values are given in Table II. We take results from fits using 5 exponentials (N exp ¼ 5 in the fit forms below), where good χ 2 values are obtained and the ground-state parameters and their uncertainties have stabilized.
Two-point correlation functions are fit to the form are fit parameters. The second term in Eq. (15) accounts for the opposite-parity states that arise from the staggered quark time doublers and are known as oscillating states (see Appendix G of [40]). These oscillating states do not appear when M is a Goldstone-taste 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. Figure 1 shows the quality of our results. We plot effective energies and amplitudes for the D s and H s mesons on the fine ensemble. For the effective energy E eff ðtÞ we use where CðtÞ is defined in Eqs. (8) and (9). By comparing CðtÞ to Cðt þ 2Þ and Cðt − 2Þ, we avoid contamination by the oscillating states. Similarly an effective amplitude can be computed according to As can be seen from Fig. 1, these quantities agree well with the results for the ground-state energies and amplitudes found from our multiexponential Bayesian fits. For the three-point correlation functions we use the fit form This includes fit parameters common to the fits of H s (when J ¼ S),Ĥ s (when J ¼ V 0 ) and D s two-point correlators, along with new fit parameters J jk that are related to the current matrix elements. We perform a single simultaneous fit containing each correlator computed (C H s ; CĤ s ; C D s ; C η h ; C η c ; C η s ; C H c ; C S ; C V 0 ) at every am h and every jap D s j, for each ensemble.
These simultaneous fits are very large and this causes problems for the covariance matrix which must be inverted to determine χ 2 . We take two steps toward mitigating this. The first is to impose an svd (singular value decomposition) cut c svd . This replaces any eigenvalue of the covariance matrix smaller than c svd x with c svd x, where x is the largest eigenvalue of the matrix. 1 The small eigenvalues are driven to zero if the statistics available are not high enough [62]. The application of the svd cut makes the matrix less singular, and can be considered a conservative move since the only possible effect on the error of the final results is to inflate them. An appropriate value for c svd is found by comparing estimates of covariance matrix eigenvalues between different bootstrap samples of the data using the Corrfitter package [61]. The resulting c svd varies between ensembles since it depends on the statistical quality of the dataset, but we find them to be of order 10 −3 . For more details of this approach, see Appendix D of [62].
The other step we take toward a stable fit is employing a chained-fitting approach. We first perform an array of smaller fits, each fitting the correlators relevant only to one m h and one jap D s j value. In the case of set 4, for example, this results in 11 separate fits. Then, a full simultaneous fit of all of the correlators is carried out, using as priors the results of the smaller fits. This both speeds up the full fit and improves stability of the results.
The priors for the fits were set up as follows. We set gaussian priors for the parameters J jk , and log-normal priors for amplitudes a M i , ground-state energies E M 0 , and excited-state energy differences E M iþ1 − E M i . Using lognormal distributions forbids ground-state energies, excited state energy differences and amplitudes moving too close to zero or becoming negative, improving stability of the fit.
Priors for ground state energies E M 0 and amplitudes a M 0 are set according to an empirical-Bayes approach, plots of the effective amplitude of the correlation functions are inspected to deduce reasonable priors. The ground-state oscillating parameters a M;o 0 , E M;o 0 , are given the same priors as the non-oscillating states, with uncertainties inflated by 50%. The resulting priors always have a standard deviation at least 10 times that of the final result. The logs of the excited-state energy differences are given prior values 2aΛ QCD AE aΛ QCD where Λ QCD was taken as 0.5 GeV. The log of oscillating and non-oscillating excited state amplitudes are given priors of −1.9 AE 3.3. The ground-state nonoscillating to nonoscillating three-point parameter, J nn 00 is given a prior of 1 AE 0.5, and the rest of the three-point parameters J nn jk are given 0 AE 1. The physical quantities that we need here are extracted from the ground-state fit parameters and given in Tables III  and IV (20) and (21). D s energies at non-zero spatial momentum are given in IfΦ M is a γ 5 ⊗ γ 5 pseudoscalar operator P, the decay constant can be found from this via where q, q 0 are the quark flavors that M is charged under. We use this to determine the H c meson decay constant in Table III. The current matrix elements that we are focused on here can be extracted from the fit parameters via These can be converted into values for the form factors once the currents have been normalized (Sect. II D). Figure 2 shows the results of a number of tests we performed on the fits to correlators on the fine ensemble. Each test modifies one of the features of the fits and we then plot the resultant value of the key output parameter J nn 00 . The robustness of the fits can be gauged by the effect of these changes, which are all small. Figure 3 shows a comparison of the D s meson dispersion relation on the fine (set 1) and superfine (set 3) lattices. The dispersion relation is sensitive to discretization effects in our quark action. The figure shows them to be small (see [63] for more discussion of discretization effects in dispersion relations for mesons using HISQ quarks). These discretization effects, along with all other discretization effects, are dealt with systematically via the continuum extrapolation detailed in Sec. II E.  (22) and (24). The final two columns give the ratio of these form factors to f H c ffiffiffiffiffiffiffiffiffi M H c p . We show statistical/fit errors on each of the quantities, including the value of q 2 that is derived from M H s and E D s . q 2 and R s uncertainties also include those from the determination of the lattice spacing.

D. Current normalization
In the HISQ formalism, the local scalar current ð1 ⊗ 1Þ (multiplied by the mass difference of flavors it is charged under) is conserved, and hence requires no renormalization. This is not the case for the local temporal vector current ðγ 0 ⊗ γ 0 Þ. We use this instead of the conserved vector current because it is much simpler, but we then require a renormalization factor to match to the continuum current. This is simple to obtain fully nonperturbatively within this calculation [31,64], at no additional cost.
When both meson states in the matrix elements are at rest (the zero recoil point), the scalar and local vector matrix elements are related via the PCVC relation: Z V can be extracted from this relation using the matrix elements we have computed. The Z V values found on each ensemble and for each am val h0 are given in Table V. We also remove Oðam 4 h Þ tree-level mass-dependent discretization effects from the current using a normalization constant, Z disc derived in [65] and discussed in detail in [49]. Z disc values are also tabulated in Table V; they have only a very small effect.
Combining these normalizations with the lattice current from the simultaneous correlation function fits, we find values for the form factors at a given heavy mass, lattice spacing, and q 2 : Tests on the correlator fits on the fine ensemble. The y-axis shows the best fit result for J nn 00 , with the appropriate current, heavy mass and q 2 specified. At N test ¼ 1 we give our final result, reproduced by the light grey band for ease of comparison. N test ¼ 2 and 3 give the results of setting N exp ¼ 4 and 6 respectively. N test ¼ 4 gives the result of setting t cut ¼ 3 for all 2-point correlators (in the final fit t cut ¼ 2 for all correlators). N test ¼ 5 gives the value when the prior width on the J nn 00 parameters is doubled. N test ¼ 6 gives the results from the output from a fit to the appropriate correlators from that heavy mass and q 2 value only, and therefore not including correlation with results from other masses and momentum values. 074513-9 where δ M is defined in Eq. (6) and we have made the dependence of the matrix elements on q 2 explicit.

E. Obtaining a result at the physical point
We now discuss how we fit our results for f s 0 ðq 2 Þ and f s þ ðq 2 Þ as a function of valence heavy quark mass, sea light quark mass and lattice spacing. Evaluating these fits at the mass of the b, with physical l, s and c masses and zero lattice spacing will then give us the physical form factor curves from which to determine the differential decay rate, using Eq. (2).
Following [49] we use two methods; one a direct approach to fitting the form factors and the other in which we fit the ratio in which discretization effects are somewhat reduced. We will take our final result from the direct approach and we describe that here. We use the ratio approach as a test of uncertainties and we describe that in more detail in Appendix B. We use identical fit functions for both approaches. We feed into the fit our results from Tables III and IV, retaining  the correlations (not shown in the Tables) between values for different heavy quark masses and q 2 values on a given gluon field ensemble that we are able to capture in our simultaneous fits (Sec. II C). We also include, where needed, correlated lattice spacing uncertainties.

Kinematic behavior
Our fit form is a modified version of the Bourrely-Caprini-Lellouch (BCL) parametrization for pseudoscalarto-pseudoscalar form factors [66]: The function zðq 2 Þ maps q 2 to a small region inside the unit circle on the complex q 2 plane, defined by Here t þ ¼ ðM H s þ M D s Þ 2 and we choose t 0 to be t 0 ¼ 0. This t 0 choice means that q 2 ¼ 0 maps to z ¼ 0 and the fit functions simplify to f s 0;þ ð0Þ ¼ a 0;þ 0 . For the physical range of q 2 for B s to D s decay, the range covered by z is jzj < 0.06, resulting in a rapidly converging series in powers of z. We truncate at N ¼ 3; adding further powers of z n does not effect the results of the fit.
The factors in front of the sums in the BCL parametrization account for lowest mass pole expected in the full q 2 plane for each form factor coming from the production of on-shell H c0 and H Ã c states in the crossed channel of the semileptonic decay. Note that these poles, even though they are below the cut (for H s þ D s production) that begins at t þ , are at much higher q 2 values than those covered by the semileptonic decay here (with maximum q 2 given by ðM H s − M D s Þ 2 ).
We must estimate M H c0 , the scalar heavy-charm meson mass, at each of the heavy masses we use. For this we use the fact that the splitting Δ 0 ¼ M H c0 − M H c is an orbital excitation and therefore largely independent of the heavy quark mass. The splitting has been calculated in [67] to be Δ 0 ¼ 0.429ð13Þ GeV at the b quark mass. Combined with an H c mass from our lattice results, we construct the H c0 mass as M H c0 ¼ M H c þ Δ 0 . We do not include the uncertainty on Δ 0 in the fit, since any shift in the precise position of the pole will be absorbed into the other fit parameters.
To estimate M H Ã c , the vector heavy-charm meson mass, we use the fact that the hyperfine splitting M H Ã c − M H c should vanish in the infinite m h limit. M H Ã c then takes the

Heavy quark mass and discretization effects
To account for dependence on the heavy quark mass and discretization effects in a general way, we use the following form for each of the a 0;þ n coefficients: To understand this form, focus first on the terms inside the sum. Powers of ð2Λ QCD =M η h Þ allow for variation of the coefficients as the heavy quark mass changes, using an HQET-inspired form since this is a heavy-light to heavylight meson transition. M η h =2 is proportional to m h at leading order in HQET, so is a suitable physical proxy for the heavy quark mass. We take Λ QCD here to be 0.5 GeV. The other two terms in the sum allow for discretization effects. These can be set by two scales. One is the variable heavy quark mass am val h0 and the other is the charm quark mass, am val c0 , constant on a given ensemble. Adding further discretization effects set by smaller scales such as aΛ QCD had no impact on the results since such effects are subsumed into the larger am val c0 terms. The coefficients d 0;þ ijkn are fit parameters given Gaussian prior distributions of 0 AE 2.
To account for any possible logarithmic dependence on m h , arising from, for example, an ultraviolet matching between HQET and QCD, we include a log term in front of the sum. ρ 0;þ n are fit parameters with prior distribution 0 AE 1. The fact that f s þ ð0Þ ¼ f s 0 ð0Þð⇒ a þ 0 ¼ a 0 0 Þ is a powerful constraint within the heavy-HISQ approach. Since this relation must be true at all m h , it translates to constraints on the fit parameters; d þ i000 ¼ d 0 i000 ∀ i and ρ þ 0 ¼ ρ 0 0 .We impose these constraints in the fit.

Quark mass mistuning
To account for any possible mistunings in the c, s and l quark masses, we include the terms N 0;þ mistuning;n in each a 0;þ M η s is the mass of an unphysical pseudoscalar ss meson where the valence quarks are artificially not allowed to annihilate; its mass (in the continuum and chiral limits) can be determined in lattice QCD calculations in terms of the masses of the π and K mesons. In [54] this analysis was done to yield M η s ¼ 0.6885ð20Þ GeV. This is the value that we use here for M phys η s to tune the s quark mass. M phys η c is taken to be 2.986(3) GeV from [54]. This allows for a small (0.1%) adjustment for cc annihilation and QED effects, neither of which are included in our calculation.
We similarly account for (sea) light quark mass mistuning by defining δ l ¼ m l0 − m tuned l . We find m tuned l from m tuned s , using the fact that the ratio of quark masses is regularization independent, and was determined in [45]: We set m tuned l to m tuned s divided by this ratio. All higher order contributions, such as δ 2 s;l , ðM η c − M phys η c Þ 2 , or ðΛ QCD =M η h Þ 2 are too small to be resolved by our lattice data, so are not included in the fit.
In our lattice QCD calculation we set m u ¼ m d ≡ m l ; this means that our results do not allow for strong-isospin breaking in the sea quarks. By moving the m tuned l value up and down by the PDG value for m d − m u [28], we found that any impact of strong-isospin breaking on our results was negligible.

Finite-volume and topology freezing effects
We expect finite-volume effects to be negligible in our calculation and we do not include any associated error. Finite-volume corrections to the B → Dlν form factors were calculated in chiral perturbation theory in [33] and found to be very small, less than one part in 10 4 , for typical lattice QCD calculations. For B s → D s lν form factors we expect finite-volume effects to be smaller than this because there are no valence u=d quarks.
The finest lattices that we use here have been shown to have only a slow variation of topological charge in Monte Carlo time. This means that averaging results over the ensemble could introduce a bias if the quantities we are studying are sensitive to topological charge. A study calculating the adjustment needed to allow for this gives only a 0.002% effect on the D s decay constant [68]. For B s to D s form factors we might expect an effect of similar relative size. This is negligible compared to our other uncertainties.

Uncertainties in the physical point
Once we have fit our lattice results as described above, we can determine the physical form factors by We take the experimental value for M η b but allow for an additional AE10 MeV uncertainty beyond the experimental uncertainty, since our lattice QCD results do not allow for QED effects or for η b annihilation to gluons [44]. This additional uncertainty has no effect, however, because the heavy quark mass dependence is mild.

III. RESULTS AND DISCUSSION
In Tables III and IV, we give our results for the form factors on each ensemble along with the meson masses needed for the fits of the form factors as a function of m h and a discussed in Sec. II E.
We first show results from simplified fits to zero recoil data to find f s 0;þ ðq 2 max Þ. This allows us to test the behavior in m h . We then perform the larger fit, described in Sec. II E, taking into account all the lattice data throughout the q 2 range.

A. Zero recoil
We performed a fit to f s 0 ðq 2 max Þ as a function of m h and a using the fit form This is the same fit function as described earlier for the individual z-space coefficients in Eq. (28) and we take the same priors for the corresponding coefficients as discussed there. The fit has χ 2 =N dof ¼ 0.21, for 12 degrees of freedom. Evaluating the result at a ¼ 0 and physical b quark mass, we find We show the dependence on M η h of our results and the fit in Fig. 4. The error budget corresponding to Eq. (33) is given in Table VI. Note that we do not impose the constraint that f s 0 ðq 2 max Þ ¼ 1.0 when m h ¼ m c . If we do this, we reduce the uncertainty in Eq. (33) by 25%.
We include in Fig. 4 a previous lattice determination of f s 0 ðq 2 max Þ [35], shown as a red triangle. Our result, containing independent uncertainties, is in agreement with this earlier value but much more accurate. The older study used the n f ¼ 2 þ 1 MILC asqtad gluon ensembles, with HISQ s and c valence quarks, and an NRQCD b quark. Using NRQCD meant that the calculation could be performed directly at the physical b mass. However, the matching of lattice NRQCD currents to continuum QCD, performed at Oðα s Þ, is a significant source of systematic error absent in our calculation.
We perform a number of tests of the fit at zero recoil, and present results in Fig. 5. The tests show that the fits are robust.

B. Full q 2 range
We now proceed to fit our full set of data including zero and nonzero recoil points to the fit form given in Eqs. (26) and (28) and discussed in Sec. II E. We include the Notice that the y-axis scale does not begin at zero. We also include the result from a previous lattice calculation, which used the NRQCD discretization for the b quark with a nonrelativistic expansion of the current through OðΛ=m b Þ and Oðα s Þ matching to continuum QCD [35]. Sets of gluon field configurations listed in the legend follow the order of sets in Table I. covariance matrix between the form factor values obtained from correlator fits on a given ensemble along with correlated lattice spacing uncertainties. The goodness of fit obtained is χ 2 =N dof ¼ 0.51, with 58 degrees of freedom. The fit parameters that are constrained by the fit are listed in Table X of Appendix C.
In Fig. 6 we show our results and fit function in z-space for the product P × f of each form factor and its pole factor, P, given by 1 − q 2 =M H 0;Ã c in Eq. (26). This shows that the zdependence of P × f is relatively benign for both form factors and the main m h -dependent effect is the smooth reduction in value of P × f as m h increases. The final result at the physical b quark mass is given by the grey band.
In Fig. 7, we show the results and fit function in q 2space. The form factors for the physical b quark mass (i.e., those corresponding to B s → D s decay) are given by the grey band.
From Fig. 6 it is clear that the largest effect in the lattice results is a z-independent shift as a function of heavy quark mass. Not surprisingly then, the parameters that are best constrained by the fit are d 0 i000 ¼ d þ i000 and ρ 0 0 ¼ ρ þ 0 , i.e., the parameters that control this heavy-quark mass extrapolation. We find that d 0;þ 1000 ¼ 1.397ð82Þ and ρ 0;þ 0 ¼ 0.419ð20Þ. Other coefficients are not as well constrained by the fit, including those that allow for discretization effects. Allowing for such terms in the fit, however, means that their impact on the final uncertainty is included. Figure 8 shows the physical f þ and f 0 form factors on the same plot and covering the full q 2 range for the B s → D s decay. Figure 9 plots the associated error budget for the two form factors throughout the q 2 range. The dominant uncertainty comes from statistical errors. There are also significant uncertainties from the q 2 and m h dependence for FIG. 6. Pf s 0;þ in z-space, where P is the appropriate pole function for each form factor given in Eq. (26). The colored points show lattice results, i.e., outputs from the correlator fits. The colors correspond to the legend given in Fig. 7. Sets listed in this legend follow the order of sets in Table I. The lowest grey band shows the result of our fit at a ¼ 0 and physical l,s, c and b masses. Each of the higher grey bands show the fit form evaluated at the heavy quark masses, lattice spacings and l,s and c masses of each of our sets of lattice data. f þ at larger values of q 2 . This is because there are no lattice QCD results at q 2 max for f þ . The impact of uncertainties in the lattice spacing (both in w 0 =a and in w 0 ) are smaller than the errors shown in the figure and so not plotted there. This is because the form factors themselves are dimensionless and lattice spacing effects in the determination of q 2 largely cancel as, for example, the pole masses are given in terms of lattice masses.
As discussed in Sec. II E an alternative approach to the fit is to take ratios of the form factors to the H c decay constant and fit the ratios to the fit form of Eqs. (26) and (28). This fit is described in Appendix B. It has the advantage of smaller discretization effects but the disadvantage of larger lattice spacing uncertainties because the ratios being fit are dimensionful. In the end the ratio method has larger uncertainty for the final physical form factors. We therefore take the results from the direct method as our final result, and use the ratio method results as a consistency test. Since the two approaches have quite different systematic errors, their comparison supplies a strong consistency check. In Fig. 10, we plot the form factors from the two methods on top of each other. As is clear from this plot, the results are in good agreement. The direct method gives a more accurate result for both form factors and at all q 2 .
We compare the coefficients from our fits to unitarity bounds in Appendix D as a further test.  Table I. In Fig. 11, we compare our final form factors to those determined from the lattice QCD calculation using the NRQCD approach for the b quark already used as a comparison at q 2 max in Fig. 4 [35]. The NRQCD calculation works directly at the b quark mass but on relatively coarse lattices and hence is unable to obtain results at large physical momenta for the D s meson. The results close to zero-recoil are extrapolated to q 2 ¼ 0 using a z-space parametrization. As the figure shows, our results are in excellent agreement with the NRQCD calculation but are more precise for both f s 0 ðq 2 Þ and f s þ ðq 2 Þ throughout all q 2 . This is because we can avoid the significant systematic uncertainty that the NRQCD calculation has from the perturbative matching to continuum QCD of the NRQCD current that couples to the W.
The B s → D s lν form factors have also recently been obtained by the Fermilab Lattice/MILC collaborations from their earlier B → D form factors and ratios of B s → D s and B → D form factors using the Fermilab formalism for b and c quarks in [69]. Our results are consistent with theirs (shown in Fig. 20 of [69]) but we have a smaller uncertainty at q 2 ¼ 0.
C. RðD s Þ Using our calculated form factors f s 0;þ ðq 2 Þ, we can calculate the differential rate for B s → D s lν decay from Equation (1). This is a function of the lepton mass and so differs between the heavy τ and the light e, μ leptons. The differential rate for μ and τ is compared in Fig. 12. We take the meson and lepton masses needed for Eq. (1) from [28] and η EW ¼ 1.011ð5Þ [23]. The distribution in the τ case is cut off at q 2 ¼ m 2 τ and so, although there is enhancement from m 2 l =q 2 terms in Eq. (1) that reflect reduced helicity suppression, the integrated branching fraction for the τ case is smaller than for the μ.
The ratio of branching fractions for semileptonic B decays to τ and to e=μ is being used as a probe of lepton universality with an interesting picture emerging [36,37]. Here we provide a new SM prediction for the quantity where l ¼ e or μ [the difference between e and μ is negligible in comparison to our precision on RðD s Þ].
Our result is FIG. 10. Results for f s 0;þ ðq 2 Þ against q 2 at the physical point, comparing the ratio method (from Appendix B) and the direct method (from Sec. III B).
FIG. 11. Our final result for f s 0;þ ðq 2 Þ compared to form factors calculated using an NRQCD action for the b quark [35]. Part of the NRQCD band is shaded darker than the rest (q 2 ⪆9.5 GeV 2 ) to signify the region where lattice results were directly calculated. The NRQCD form factors in the rest of the q 2 range are the result of an extrapolation using a BCL parametrization.
in which we averaged over the l ¼ e and l ¼ μ cases. Note that jV cb j and η EW cancel in this ratio. We give an error budget for this result in terms of the uncertainties from our lattice QCD calculation in Table VII. Our result agrees with, but is more accurate than, the previous lattice QCD value of RðD s Þ [0.301 (6)] from [35]. An experimental result for RðD s Þ would allow a new test of lepton universality.
We expect very little difference between RðD s Þ and the analogous quantity RðDÞ because the mass of the spectator quark has little effect on the form factors [34]. Lattice QCD calculations that involve light spectator quarks have larger statistical errors, however, which is why the process B s → D s is under better control. Previous lattice QCD results for RðDÞ are 0.300(8) [23] and 0.299(11) [22], in which any difference with our result for RðD s Þ is too small to be visible with these uncertainties.

IV. COMPARISON TO HQET
In Fig. 13 we show our form factor results at two key values of q 2 , the zero recoil point and q 2 ¼ 0, as a function of heavy quark mass, given by M η h . The plot demonstrates how f þ at zero-recoil increases as the heavy quark mass increases, but f 0 changes very little. The value at q 2 ¼ 0, where the form factors are equal, falls with growing heavy quark mass, as the q 2 range opens up.
Knowledge of the functional form of f s 0 ðq 2 Þ and f s þ ðq 2 Þ, along with h s A 1 ð1Þ from [49], against m h gives us access to the functional form in m h of a number of quantities of interest in HQET.
In HQET the vector current matrix element is parametrized with a different set of form factors, h s þ ðwÞ and h s − ðwÞ according to where v μ ¼ p μ H s =M H s and v 0μ ¼ p μ D s =M D s are the 4-velocities of the initial and final state mesons, and w ¼ v · v 0 is an alternative parameter to q 2 used in the context of HQET.
As a test of HQET, one can construct ratios of form factors that should become unity in the m c ; m h → ∞ limit. Following the notation defined in [70], one can redefine the form factors such that each of them reduce to the Isgur-Wise function ξðwÞ in the m c ; m h → ∞ limit. In the B s → D s case these new form factors are where r ¼ M D s =M H s . h s A 1 also reduces to ξ in the infinite mass limit. Hence any ratio between S s 1 ,V s 1 and h s A 1 should FIG. 14. Form factor ratios against M η h , a proxy for the heavy quark mass. S s 1 and V s 1 are defined in Eqs. (37) and (38). The colorful points are NLO HQET expectations from [25], derived with input from QCD sum rules. become unity in this limit. From our results at zero recoil and for m h ¼ m b we find Figure 14 illustrates how these ratios vary with m h , and gives the NLO HQET expectation for these values for comparison [25]. As can be seen here, our results show large deviations from the HQET expectation, implying that NLO HQET is missing significant higher order contributions. As discussed in [71], previous lattice results have shown similarly large deviations. These are also shown in preliminary results from the JLQCD collaboration [72].
Another set of quantities of interest in HQET are the slopes of the form factors at q 2 ¼ 0 [73,74]: ; ð42Þ To obtain these values from our results for the form factors, we take the derivative of the fit function [Eq. (28)] evaluated at continuum and physical l, s and c masses.
In Fig. 15 we show how these quantities vary with m h . By rewriting Eq. (43) in terms of h s þ;− , and recognising that in the heavy quark limit h s þ ≈ ξ and h s − ≈ 0, one can find a leading order HQET expectation for δ [74]: This leading order expression is shown in Fig. 15 as a blue band. Our results (grey band) are in good agreement with this, up to the uncertainties from higher order terms shown in Eq. (45).

V. CONCLUSIONS
We have calculated the scalar and vector form factors for the B s → D s lν decay for the full q 2 range using lattice QCD with a fully nonperturbative normalization of the current operators for the first time. Our calculation used correlation functions at three values of the lattice spacing, including an ensemble with an approximately physical light quark mass. We used the relativistic HISQ action with a range of values for the heavy valence quark and by fitting this dependence, along with the lattice spacing dependence, we are able to determine results at the b quark mass. The valence c and s quark masses are accurately tuned to their physical values. By working on very fine lattices we are able both to reach a heavy quark mass close to the b but also to cover the full q 2 range of the decay.
Our results for the form factors are given in Fig. 8 and the differential rate that this implies for B s → D s lν in Fig. 12. This will allow a determination of jV cb j from future experimental data from this semileptonic process. Instructions on how to reproduce our form factors are given in Appendix A. Our error budget is given in Fig. 9.
Our results are more accurate than previous lattice QCD determinations using a nonrelativistic b quark formalism because we do not have a systematic uncertainty from the perturbative matching of the lattice current to continuum QCD.
From our results we can predict the ratio RðD s Þ of the branching fraction to a τ lepton compared to that to e=μ (see Sec. III C). We are also able to compare functions of the form factors and their slopes to HQET expectations (see Sec. IV).
Our calculation shows that a heavy-HISQ determination of the B → Dlν form factors is feasible. This would allow direct comparison to existing experimental data. Such a FIG. 15. We show two quantities derived from the form factor slopes as a function of M η h . 1=βðm h Þ is defined in Eq. (42) and δ in Eq. (43). Our results are shown by the grey bands. The blue band shows the leading order HQET expectation for δ given in Eq. (45). calculation could use an identical strategy to the one demonstrated here, with the strange valence quark replaced with a light one and additional calculations on ensembles spanning a range of light quark masses. Higher statistics would be needed since statistical uncertainties will be larger than in this calculation, and the issue of topology freezing on fine lattices will be more significant. Our calculation demonstrates, however, that lattice QCD is no longer limited for these form factors by the systematic uncertainties coming from current matching and, with sufficient computer time, a 1% accurate result for B → Dlν form factors is achievable.    Table I. from the amplitude for the ground-state in the two-point correlation functions as given in Eq. (21). Our results for R s 0;þ ðq 2 Þ are given in Table IV. We fit R s 0;þ ðq 2 Þ using an identical fit function to that of the direct approach, given in Eqs. (26) and (28). The fit parameters are tabulated in Table XI in Appendix C. Discretization effects change in the ratio given in Eq. (B1), compared to those from the form factors themselves, changing the continuum extrapolation. The dependence on heavy quark mass of the ratio is also very different. The value of the ratio at the physical point (where m h ¼ m b and a ¼ 0) can then multiplied by f B c ffiffiffiffiffiffiffiffi ffi M B c p to obtain the form factors. We find f B c via a separate calculation (detailed in Appendix A of [49]) and take the experimental value for the B c meson mass, M B c ¼ 6.2756ð11Þ GeV [28]. This approach has the disadvantage of introducing larger uncertainties from scale-setting since R s 0;þ ðq 2 Þ are dimensionful quantities (as opposed to f s 0;þ ðq 2 Þ which are dimensionless). Hence we do not use it for our final value. It provides a useful test of our uncertainties, however.

Zero recoil
We first show results from a fit at the zero recoil point, R s 0 ðq 2 max Þ, as a function of heavy quark mass and lattice spacing. To do this we use the same fit form as for our fits to f s 0;þ ðq 2 max Þ, given in Eq. (32), with the same priors. We find, with physical parameters for all quark masses and a ¼ 0 The fit has χ 2 =N dof ¼ 0.63 for N dof ¼ 16. The lattice QCD results and fit are shown in a plot against M η h in Fig. 16. As can be seen from this plot, the lattice results have stronger dependence on the heavy quark mass and somewhat less on the lattice spacing compared to that seen in Fig. 4. The error budget for our physical value in Eq. (B2) is given in Table IX. Notice that, compared to Table VI it now has a significant scale-setting uncertainty.

Full q 2 range
In Fig. 17, we show the complete set of lattice results along with the results of the full fit given by the fit form in Eq. (28). As for Fig. 7 we see that as the lattice spacing decreases, the range of heavy quark masses increases and the q 2 range, 0 < ðM H s − M D s Þ 2 expands. The goodness of fit here was χ 2 =N dof ¼ 0.57, N dof ¼ 58.
To obtain the form factors, the resulting functions R s 0;þ ðq 2 Þ were multiplied by ffiffiffiffiffiffiffiffi ffi M B c p (using the experimental value) and f B c from our determination detailed in Appendix A of [49]. The resulting form factors are shown in Figure 10 in a comparison to those found by our direct method of Sec. III B.

APPENDIX C: FIT OUTPUTS
Here we give the results for the free parameters from our fit described in Sec. II E, with Eqs. (28) and (29). Most of the parameters in these equations are not well determined by the fit. The point of including them in the fit is then to make sure that the uncertainty from not knowing them is Þ plotted against q 2 . The grey band shows the result of our fit at a ¼ 0 and physical l,s, c and b masses. Sets listed in the legend follow the order of sets in Table I. fed into our final results. In these cases the fit simply returns the prior value for that parameter, and we do not list these.
The parameters that can be determined from the fit are mainly those associated with the leading n ¼ 0 coefficient in the z-expansion. Some information about the n ¼ 1 term is also obtained, but none of the n ¼ 2 coefficients are constrained. Table X gives the mean values and uncertainties of the coefficients d 0 ijkn and d þ ijkn from Eq. (28) for n ¼ 0 and n ¼ 1 that are constrained from the fit. The d þ coefficients for n ¼ 0 and j ¼ k ¼ 0 are the same as those for d 0 , to implement the constraint that f þ ð0Þ ¼ f 0 ð0Þ. The values of the ρ parameters obtained for n ¼ 0 and n ¼ 1 (n ¼ 2 is not constrained) are also given in Table X. For n ¼ 0 ρ 0 0 ¼ ρ þ 0 , again to enforce f 0 ð0Þ ¼ f þ ð0Þ. The mass mistuning parameters given in Eq. (29) are not constrained by the fit. Table XI gives the same information for the alternative ratio fit described in Appendix B.

APPENDIX D: TESTS OF UNITARITY BOUNDS
Unitarity and crossing symmetry imposes bounds on the coefficients of the BCL parametrization of f 0;þ ðq 2 Þ, fa n g [75,76]. As another consistency test, we show here that the coefficients for f þ found in our fit satisfy these bounds.
To obtain bounds on the BCL coefficients [66], one must relate them to those of a different parametrization, that of Boyd, Grinstein and Lebed (BGL) [77]: BðzÞ is known as the Blashke factor: where z Ã ¼ zðM 2 ðD3Þ In the f s þ case, κ ¼ 6πM 2 B s χ V , p ¼ 3, s ¼ 2. The quantity χ V is the once-subtracted dispersion relation at q 2 ¼ 0 for the vector b → c current, computed in [77] to be χ V ¼ 5.7 × 10 −3 =m 2 b . The BGL coefficients, fb n g, obey the unitarity constraint by construction of the parametrization. To see how this applies to the BCL coefficients fa n g, one must relate them to fb m g by equating the two parametrizations to find [66] X M m¼0 b m z m ¼ ψðzÞ X N n¼0 a n z n ; ðD5Þ TABLE X. Fit results for the parameters d 0 ijkn , d þ ijkn ρ 0 n and ρ þ n defined in Eq. (28) and used in our preferred fit described in Sec. II E. We give the mean values and uncertainties, but do not include the correlation matrix. The n ¼ 0 parameters for j ¼ k ¼ 0 are the same for d 0 and d þ , as are the n ¼ 0 parameters for ρ 0 and ρ þ . We restrict the list to those parameters that are constrained by the fit beyond the prior value of 0 AE 2, but include values that are not constrained for n ¼ 1 where the n ¼ 0 result is constrained. n ¼ 2 parameters are also included in the fit but they are not constrained.  (20) −0.23ð96Þ 0.419 (20) 0.08 (25) TABLE XI. Fit results for the parameters d 0 ijkn , d þ ijkn , ρ 0 n and ρ þ n defined in Eq. (28) but used with the alternative ratio fit described in Appendix B. We give the mean values and uncertainties, but do not include the correlation matrix. The n ¼ 0 parameters for j ¼ k ¼ 0 are the same for d 0 and d þ , as are the n ¼ 0 parameters for ρ 0 and ρ þ . We restrict the list to those parameters that are constrained by the fit beyond the prior value of 0 AE 2, but include values that are not constrained for n ¼ 1 where the n ¼ 0 result is constrained. n ¼ 2 parameters are also included in the fit but they are not constrained.  (21) where ψðzÞ is given by and M pole ¼ M B Ã c in the f s þ case. Expanding ψðzÞ around z ¼ 0, comparing coefficients of z in Eq. (D5), and imposing the constraint of Eq. (D4), we arrive at a constraint for the BCL coefficients B ≡ X L;L j;k¼0 B jk a j a k ≤ 1; ðD7Þ where fη n g are the Taylor coefficients of ψðzÞ. ψðzÞ is bounded on the closed disk jzj < 1, so its Taylor coefficients are rapidly decreasing. We computed values for B jk by truncating the sum in its definition [Eq. (D8)] at 100.
These values are given in Table XII. With these B jk values, and the a n coefficients from our fit (via the direct method), we find B þ ¼ 0.0210ð54Þ: This comfortably satisfies the unitarity bound. Additionally, as discussed in [78], the leading contributions to B þ are of order ðΛ QCD =m b Þ 3 ≃ 10 −3 in the heavy quark expansion. This expectation is approximately fulfilled by our result.