Toward accurate form factors for B -to-light meson decay from lattice QCD

We present the results of a lattice QCD calculation of the scalar and vector form factors for the unphysical B s → η s decay, over the full physical range of q 2 . This is a useful testing ground both for lattice QCD and for our wider understanding of the behavior of form factors. Calculations were performed using the highly improved staggered quark (HISQ) action on N f ¼ 2 þ 1 þ 1 gluon ensembles generated by the MILC Collaboration with an improved gluon action and HISQ sea quarks. We use three lattice spacings and a range of heavyquarkmasses from that of charm to bottom, allin the HISQ formalism. This permits an extrapolation in the heavy quark mass and lattice spacing to the physical point and nonperturbative renormalization of the vector matrix element on the lattice. We find results in good agreement with previous work using nonrelativistic QCD b quarks and with reduced errors at low q 2 , supporting the effectiveness of our heavy HISQ technique as a method for calculating form factors involving heavy quarks. A comparison with results for other decays related by SU(3) flavor symmetry shows that the impact of changing the light daughter quark is substantial but changing the spectator quark has very little effect. We also map out form factor shape parameters as a function of heavy quark mass and compare to heavy quark effective theory expectations for mass scaling at low and high recoil. This work represents an important step in the progression from previous work on heavy-to-heavy decays ( b → c ) to the numerically more challenging heavy-to-light decays. DOI:


I. INTRODUCTION
Determinations of form factors for weak semileptonic meson decays can be combined with experimental results to provide important tests of the Standard Model (SM). Decays of b quarks are of particular interest as they allow determination of some of the least well-known elements of the Cabibbo-Kobayashi-Maskawa (CKM) matrix [1,2] and tests of the unitarity of that matrix, a foundation of the weak sector of the SM. Increasingly small experimental uncertainties in CKM-dependent decay rates must be met with precise determinations of form factors from the theoretical side to pin down the CKM matrix elements (see, for example, [3,4]). The shape of the differential decay rate in q 2 , the squared momentum transfer between the initial and final states, parameterized by the form factors, provides added detail when testing the SM. Lattice quantum chromodynamics (lattice QCD) is the only modelindependent method for calculating the hadronic form factors for such decays and has been used successfully for many such calculations. For a review, see [5].
Resolving the b quark on the lattice requires a sufficiently small lattice spacing, a < 1=m b ∼ 0.05 fm. This means that lattice QCD calculations can currently only reach the b quark mass on the finest lattices available.
One approach to address this difficulty relies on the use of an effective theory description of the b quark. Examples include the relativistic heavy quark action [6], the Fermilab action [7,8], heavy quark effective theory (HQET) [9], and nonrelativistic QCD (NRQCD) [10][11][12]. Each of these must match the relevant effective theory to QCD and therefore suffer from associated matching errors. For the case of NRQCD, such matching errors are a dominant source of uncertainty.
Alternatively, with an action sufficiently improved to reduce heavy quark discretization effects, one can avoid this use of effective theory and simulate over a range of heavy quark masses m h ≲ m b and then extrapolate (or interpolate if, for example, static quark results are available) to m b . Examples of this approach include the ratio method using the twisted mass formulation [13,14], application of the Möbius domain wall formulation to the b quark [15] and our recent works using the highly improved staggered quark (HISQ) action for the b quark in several b → c decays [16][17][18].
The HISQ action [19] provides an accurate discretization of the Dirac equation for relatively heavy quarks [20]. It allows us to normalize lattice currents nonperturbatively using conserved currents, avoiding sizable systematic errors from perturbative truncation in the renormalization factors for the nonrelativistic case. This "heavy-HISQ" approach must be carried out on fine lattices, with a < 0.1 fm so that am h is not too large. On our finest lattices, am b < 1. In practice, we work at several values of a and of the heavy quark mass so that we can map out both discretization effects and physical dependence on the heavy quark mass to determine the result at m b and in the continuum. A further advantage of working on such fine lattices is that we can reach higher physical values of momentum transfer as the lattice spacing gets smaller. This is particularly important for b decays where the q 2 range for the decay is large. With the heavy-HISQ approach the range of accessible q 2 values grows on finer lattices in step with the range of heavy quark masses. This means that we can cover the full q 2 range of the heavy quark decay all the way up to that of the b [16].
The end game of this program is the determination of form factors for transitions that involve physical u and d quarks, such as B → π. In this work, we take an important step in extending our use of the HISQ action for the b quark in b → c decays toward the more demanding b → u; d decays by studying the b → s transition. As m s ≪ m c , this allows us to gauge the success of this approach for b-tolight form factors while benefiting from both a significant savings in computational cost and the typically less noisy correlators associated with s quarks. Fixing the daughter quark to the strange quark mass on each ensemble removes the need to perform a chiral extrapolation, thereby simplifying the continuum extrapolation, a key component in our study of the efficacy of the heavy-HISQ approach. Here we study the B s → η s decay, where the η s is an unphysical ss pseudoscalar meson-an easier to analyze, cheaper to compute substitute for a pion, with the same quantum numbers and no valence annihilation. For the purposes of assessing the viability of this approach, we focus on the scalar and vector form factors. 1 The form factors should not be greatly affected by changing the spectator quark from an s quark to a u=d quark, so studying this decay provides an estimate of the level of precision achievable in the computationally more expensive B → K form factor calculation. The heavy-HISQ approach allows us to extract the dependence of the form factors on the heavy quark mass as it varies from m c to m b , permitting useful tests for expectations from heavy quark symmetry.
The paper is laid out as follows. In Sec. II we set out the details of our lattice QCD calculation, including analysis of the correlation functions, normalization of the lattice currents and our fits to the form factors enabling results to be obtained for B s → η s decay in the continuum limit. Section III gives results and compares them both to expectations from heavy quark symmetry and to previous lattice QCD results for decay processes connected to B s → η s and D s → η s by SU(3) flavor symmetry, either for the active light quark in the decay or the spectator light quark. Finally, Sec. IV gives our conclusions.

A. Form factors
The aim of our calculation is to determine the matrix element for the V − A electroweak current between B s and η s mesons, hB s jV μ − A μ jη s i. Here the vector current is defined as V μ ¼ψ b γ μ ψ s and the axial vector current is A μ ¼ψ b γ 5 γ μ ψ s . For pseudoscalar to pseudoscalar decays, only contributions from the vector part of the V − A current are present, as a result of QCD parity invariance.
Our heavy-HISQ approach works by determining the B s meson matrix elements from a set of matrix elements for mesons in which the b quark is replaced by a heavy quark with mass m h < m b . We denote these pseudoscalar heavy-strange mesons generically by H s . The form factors f þ ðq 2 Þ and f 0 ðq 2 Þ that are determined from the matrix elements are a function of q 2 ¼ ðp H s − p η s Þ 2 , and we compute these across the full kinematic range, As m h → m b this becomes the full range for the B s decay.
The connection between the matrix elements of the lattice temporal vector and scalar currents and the form factors of interest, f þ ðq 2 Þ and f 0 ðq 2 Þ, is Bilinears constructed from staggered quarks have a "taste" degree of freedom and, as will be discussed below, we need to arrange the tastes of mesons and lattice currents appropriately so that tastes cancel in the calculated correlation functions. Here, in spin-taste notation [19], the lattice currents are S ¼ψ s 1 ⊗ 1ψ b and V 0 ¼ψ s γ 0 ⊗ ξ 0 ψ b and H s andĤ s denote Goldstone and local non-Goldstone heavy-strange pseudoscalar mesons, respectively. Equation (2) comes from the partially conserved vector current (PCVC) relation [21], which also leads to the renormalization of the vector matrix element [21,22] (see Sec. II D). We also require that the matrix element is analytic as q 2 → 0. We can see from Eq. (1) that this demands where we will drop the superscript from now on. Both matrix elements are calculated using a Goldstone pseudoscalar strange-strange η s bilinear, η s ¼ψ s γ 5 ⊗ ξ 5 ψ s , while the scalar uses the Goldstone pseudoscalar heavystrange H s ¼ψ b γ 5 ⊗ ξ 5 ψ s , and the vector uses the non-Goldstone pseudoscalar heavy-strangeĤ s ¼ψ b γ 5 γ 0 ⊗ ξ 5 ξ 0 ψ s . All of these operators are local, giving less noisy correlation functions than their point-split counterparts.

B. Lattice details
The calculation was run on ensembles of gluon field configurations generated by MILC [23,24]. These include in the sea two degenerate light quarks, strange and charm quarks, with masses m sea l , m sea s , and m sea c , respectively, using the HISQ action. The three ensembles used have parameters listed in Table I. The gluon action is Symanzik improved to remove discretization errors through Oðα s a 2 Þ [25]. Our calculation follows the approach in the calculation of B s → D s in [16] but with a strange daughter quark in lieu of a charm. The ensembles that we use here have unphysically heavy light quark masses (of value around 1=5 of the s quark mass). In [16], little effect was seen on the form factors from the light quark mass in the sea. We similarly expect little effect here since B s → η s does not involve any valence light quarks. Our main focus here is to test the heavy quark mass dependence and so we simply address the mistuning of sea light quark masses when we extrapolate to the physical point in Sec. II E.
We denote the heavy quark h and its mass m val h and use a range of heavy masses from the physical charm to am val h ¼ 0.8, the point where discretization errors start to become significant, on each set of gluon configurations. This allows us to perform a fit to our results as a function of heavy quark mass and obtain results at the physical b mass. At the same time we determine the dependence of the form factors on the heavy mass from the charm to the bottom with D s → η s and B s → η s at the two ends of the range. On the finest lattice am val h ¼ 0.8 is close to the physical b mass, allowing good control of the subsequent extrapolation to m b .
We choose a range of daughter momenta so as to give good coverage of the full momentum transfer range of the decay (see Table II) and implement these momenta using TABLE I. Gluon field ensembles used in this work. The Wilson flow parameter w 0 ¼ 0.1715ð9Þ fm is determined in [26], following the approach outlined in [27], and is used to calculate the lattice spacing a via values for w 0 =a, in column 3, which are from [16]. Column 4 gives the spatial (N s ) and temporal (N t ) dimensions of each lattice in lattice units, while columns 5-7 give the masses of the sea quarks.

Set
Handle  [29]. Valence heavy quark masses am val h are chosen to span the range from the physical charm, tuned as in [29], to am val h ¼ 0.8. Simulated η s momenta a ⃗p η s are fixed using twisted boundary conditions as described in the text. On each ensemble, we use n cfg configurations and n src time sources. Data are generated for multiple temporal source-sink separations T between the η s and H s mesons. twisted boundary conditions on the daughter strange quark in the η s , as described in [28]. The heavy meson remains at rest in all stages of the calculation, meaning the strange spectator and heavy quark have no twist applied. We calculate two-point correlation functions for the Goldstone pseudoscalar (γ 5 ⊗ ξ 5 ) η s and the two heavystrange bilinears detailed above. The correlators are built using where g q ðx t ; x 0 Þ is the one-spinor component staggered propagator for a quark of flavor q, from point The twist angle θ is given by θ ¼ ja ⃗pjN s =ð ffiffi ffi 3 p πÞ, with a ⃗p in the spatial (1,1,1) direction. We sum the spatial components of x t over the lattice sites to give the two-point correlation function for each 0 ≤ t ≤ aN t . The hi denotes path integration over all fields, carried out using the averaging over ensembles, and the trace is over color. Random wall sources are used at x 0 to improve statistical precision.
The local non-Goldstone pseudoscalar (γ 5 γ 0 ⊗ ξ 5 ξ 0 ) heavy-strange meson is similarly defined, but the spintaste structure is implemented using a lattice site-dependent phase: We need to use this in the threepoint correlation function with temporal vector current in order to cancel tastes. The mass of the local non-Goldstone meson only differs from that of the Goldstone by discretization effects, which are very small and disappear in the limit of zero lattice spacing. In our case the mass splittings between H s andĤ s are so small as to only be visible above the statistical uncertainty on the fine lattice.
We also calculate three-point functions, with the scalar and temporal vector current insertions as defined in Sec. II A. We place the η s operator at x 0 , the current at x t , and the relevant heavy-strange H s orĤ s at x T ¼ ðT; ⃗ x T Þ, where we again sum over spatial components. We then need extended heavy quark propagators from x T to x t for each heavy quark mass. The three-point functions combine quark propagators as T takes several different values on each lattice, detailed in Table II, and we determine correlation functions for all x t from 0 to T. The combination of propagators needed is illustrated in Fig. 1.

C. Analysis of correlation functions
We perform a simultaneous, multiexponential fit of the two-and three-point correlation function data using a standard Bayesian approach, introduced in [30] and expanded upon in [31,32]. Further detail is available in the documentation for the Gvar [33], Lsqfit [34] and Corrfitter [35] PYTHON packages used to perform the analysis.
Bias in the small eigenvalues of a large data covariance matrix with a finite data sample is addressed with a singular value decomposition (SVD) cut. This is a conservative move which avoids underestimating errors (see Appendix D of [36]). We implement the SVD cut by replacing eigenvalues smaller than the product of the cut and the largest eigenvalue with that value. We determine an appropriate SVD cut from eigenvalues of bootstrapped data, a facility which is built into Corrfitter. We check stability against doubling and halving the SVD cut compared to the recommended value and demonstrate this in Fig. 3.
Using an SVD cut and broad priors can lead to an artificial reduction in χ 2 =d:o:f: Corrfitter has a built-in facility permitting the determination of a more realistic value (see documentation [33][34][35] and Appendix D of [36]) by adding SVD and prior noise. We have checked that the fits reported below give values of χ 2 =d:o:f: close to 1 with this augmented noise. We report the raw χ 2 =d:o:f: values in Fig. 3 since they still provide a useful comparison between fits.
Bayesian fits provide an additional fit statistic, the Bayes factor, which penalizes overfitting, thereby providing a measure of fit quality complementary to χ 2 . For each fit, Corrfitter calculates the Gaussian Bayes factor (GBF), the Bayes factor under assumed Gaussian probability distributions. When evaluated together, GBF and χ 2 =d:o:f: provide a useful diagnostic for evaluating the ability of a fit to describe the data while not overfitting.
We aim to extract the ground state energies from the twopoint functions and use these, combined with ground state amplitudes, to extract ground state to ground state matrix elements from the three-point correlators.
We fit two-point correlators for a meson M to where a tower of excited states of energy E M;n i and amplitude a M;n i above the ground state (i ¼ 0) are generated by our lattice operator. Discarding data for t < t min allows us to fit a finite number N 2pt exp of these states, and t min =a takes values in the range 3-9 for different correlators and different lattice spacings. As detailed in [19], HISQ twopoint correlators also produce states which oscillate in time from lattice site to lattice site, with the exception of the zero momentum η s , where the quark and antiquark of the same mass prevent this effect from being exhibited. We determine priors for the ground state energies and amplitudes using the effective mass and effective amplitude, defined as We give each a broad uncertainty, checking that the final result of the fit is much more precisely determined than this prior. The lowest oscillating state prior is taken to be the nonoscillating ground state plus Λ QCD (which we take to be 0.5 GeV), with an error around 50%. The energy differences between all excited states are taken to be Λ QCD with an error of 50%. We use log-normal priors throughout to enforce positive values on energy splittings and amplitudes. Priors for excited state nonoscillating and all oscillating amplitudes are based on previous experience of amplitude sizes, and some are slightly adjusted to maximize the GBF; these are listed in Table III. In all cases, priors are many times broader than the final fit uncertainties, as demonstrated in Fig. 2.  We perform three-point fits to for different masses of H s (for the scalar current insertion) orĤ s (for the temporal vector current insertion) mesons and different twists of η s meson. J no ij represents the amplitude for the ith nonoscillating state of the η s and the jth oscillating state of the heavy meson. J ¼ S, V, for our scalar and vector current insertions. We create the η s at t ¼ 0, insert the current at t and annihilate the H s (Ĥ s ) at T.
Priors for J nn 00 are based on the effective three-point amplitudes, which can be determined from Priors for all other J kl ij values are listed in Table III. Figure 2 shows representative plots of the two-point and three-point correlator data, illustrating prior selection and providing a comparison of fit results with both the prior and data. The effect of doubling and halving the standard deviation given to all priors on the overall results of the fit are shown in Fig. 3.
On each ensemble, we perform a chained, marginalized fit to our two-and three-point correlator data. For detailed descriptions of chaining and marginalization, see [31,32] and the Corrfitter documentation [35].
The chained fit begins with a simultaneous fit to all of the two-point correlators (H s andĤ s for each m h and η s for each a ⃗p), fixing N 2pt exp in Eq. (9) for each lattice spacing such that it gives an acceptable χ 2 and maximizes the GBF. We take N 2pt exp ¼ 5 in the case of set 1 and N 2pt exp ¼ 6 in the case of sets 2 and 3. The next step in the chained fit is a simultaneous fit to all three-point correlators. This includes both S and V current insertions and data at the values for T chosen for each ensemble (listed in Table II). The chained fit prescription uses two-point correlator fit posteriors as priors for the twopoint fit parameters that appear in the subsequent three-point correlator fit, accounting for correlations between these posteriors and the three-point correlator data.
In the three-point correlator fits, the number of states N 3pt exp in Eq. (12) must be understood in terms of marginalization. Marginalization [31] subtracts fit function terms, evaluated using priors, from the data before performing the fit. In this way, effects from these terms are accounted for while the fit function used by the minimizer is simplified. For sets 1, 2 and 3, we choose N 3pt exp ¼ 2, 3 and 2, respectively, such that we achieve an acceptable fit (χ 2 per degree of freedom of 0.342, 0.079 and 0.047, respectively). On each set, the total number of states accounted for, either explicitly fit using Eq. (12) or subtracted from the data, is equal to N 2pt exp . For example, on set 1 we fit two-point correlators with N 2pt exp ¼ 5. For the fit to the three-point correlators, we first subtract from the data contributions from terms in Eq. (12) with i or j equal to 3, 4 or 5. We then fit this data using Eq. (12) with N 3pt exp ¼ 2. This is useful here because our threepoint data are noisier than our two-point data, so fewer states are required in their fits. Marginalization allows us to include information about higher states obtained from two-point fits.
We also check that the momentum dispersion relation for our η s fit results agrees with the momenta given in the FIG. 3. Stability tests of the chained, marginalised fit used on a typical three-point correlator. Test 0, the final result, shows the value of V nn 00 for am h ¼ 0.45, aj ⃗pj ¼ 0.1430 on set 2, with N 3pts exp ¼ 3 exponential terms and three additional states marginalised (as discussed in the text), with t min =a ¼ 2, the number of data points removed from the fit at the start and end of the data. Tests 1 and 2 show the effects of increasing and decreasing the number of fitted exponentials by 1, tests 3 and 4 show the effect of doubling and halving the SVD cut, and 5 and 6 show the effect of doubling and halving the error on all priors. Test 7 shows the effect of an increase on t min =a by 1, and test 8 shows the reduction of the marginalised exponentials from 6 to 5. Finally, test 9 shows the result of just fitting the vector 3 point correlator for this mass and twist, and the relevant 2 points; this gives a reduced error since the smaller fit requires a smaller SVD cut. Fitting like this does not preserve correlations, however, so we use a global fit. Other two and three-point correlators behaved similarly well under the same tests. The χ 2 =dof values (purple × s) are also plotted for reference. Note that these are the raw values and hence artificially small (see text) and the degrees of freedom are not the same across all tests.
lattice calculation. The two should differ by discretization effects only, which are small for the η s as it contains only s quarks but grow with momentum. This is displayed in Fig. 4.
Fit results are converted according to to extract the matrix elements which appear in the definition of the form factors [Eqs. (1) and (2)]. We always use the mass of the Goldstone heavy-strange pseudoscalar for M H s as the non-Goldstone mass is the same in the continuum limit. Tests showed that changing this to the non-Goldstone mass in the case of the vector matrix element made no difference at all to our continuum form factors, as discretization errors are accounted for in our extrapolation to the physical point. The results from the fits for each of the three lattice spacings are summarized in the Appendix in Tables VI-VIII.

D. Current normalization
The PCVC relation, applied at zero spatial momentum for the daughter meson, allows us to normalize the vector matrix element nonperturbatively using the scalar matrix element [16,22,37]. This uses the fact that this current is conserved in the HISQ formalism, that is to say that the product of the bare mass and the scalar matrix element does not require renormalization. We also make the small correction Z disc to account for small tree-level mass-dependent discretization effects beginning at order ðam h Þ 4 . For the determination of Z disc see [38,39]. Values for these normalizations can be found in Table IV.

E. Continuum and quark mass extrapolation
Having calculated f 0 ðq 2 Þ and f þ ðq 2 Þ for the three lattice spacings and at a range of heavy quark masses and q 2 values on each lattice, we now perform a fit in heavy quark mass, sea quark mass and lattice spacing. We can then evaluate our form factors at the physical quark masses and zero lattice spacing. Our fits also allow us to examine the heavy quark mass dependence of the form factors, all the way down to the charm mass.

Fit ansatz and priors
Following the method successfully employed in [16], we fit the form factors on the lattice using the Bourrely-Caprini-Lellouch parameterization [40]: where we use a mapping of q 2 to z, a region inside the unit circle of the z plane,   (1) and (2). Z V is calculated as in Eq. (15) and Z disc is defined in [38].
These have the correct limit at physical quark mass values. We choose to take t 0 ¼ 0. To fit the data for f 0 ðq 2 Þ and f þ ðq 2 Þ, tabulated in the Appendix, we calculate for each quark mass and momentum simulated the corresponding value of z, using the associated meson masses and values of q 2 .
The poles in Eq. (16) account for the production of on shell H s0 and H Ã s states for q 2 > q 2 max , and the mass M H s0 is taken as M H s þ 0.4 GeV, which is consistent with lattice results in [41] and experimental results [42] for the axial vector-vector splitting, M B s ð1 þ Þ − M B s ð1 − Þ. We do not need to know this number precisely as we are simply removing a pole in the data to ease the fitting process and then replacing it later. Indeed, excluding the pole from the f 0 fit function completely leads to fit results which are consistent with those from including the pole. The position of the M H Ã s pole can be estimated, as in [16], using the fact We go one step further to ensure that this ansatz also gives the correct PDG value for M phys with splittings ΔðB s Þ ¼ 0.0489ð15Þ GeV and ΔðD s Þ ¼ 0.14386ð41Þ GeV, from the PDG. We find no significant difference in the final form factors from the change of ansatz, supporting our assertion that the exact pole position is not crucial, as any small errors here are accounted for by higher orders of the z expansion. We use N ¼ 3 in Eq. (16) for our final results. We fit coefficients a 0;þ n to a general fit form, accounting for heavy quark mass dependence and discretization effects: We use M H s as a physical proxy for the heavy quark mass, as the two are equal at leading order in HQET. Terms in Λ QCD =M H s (with Λ QCD ¼ 0.5 GeV) parameterize the effect of changing heavy mass, while the other terms in the sum allow for discretization effects, which for the HISQ action appear as even powers of energy scales. In this case the two relevant energies are the heavy quark mass and Λ QCD . The log term comes from the matching of our HQET-inspired fit function to QCD [43,44]. From [43], we expect the coefficient of the log term to be of order unity, so we use a prior of 0 AE 1.
We find M η c on the three sets from [16] and take M phys GeV. This value differs from the experimental η c mass [42] by 7 MeV to allow for the effect determined in [46] of missing quark-line disconnected diagrams in the lattice calculation of the η c mass.
We give all d coefficients a prior of 0 AE 1, with the exception of d i10n , which multiply terms with ðam h Þ 2 in them. Since the HISQ action is improved up to second order in the lattice spacing, we expect these terms to be small, giving them a prior of 0.0 AE 0.5. We set d þ i000 ¼ d 0 i000 and ρ þ 0 ¼ ρ 0 0 to enforce f 0 ð0Þ ¼ f þ ð0Þ on the fit, in the continuum and in the absence of quark mistuning, although relaxing this constraint still leaves the two values agreeing within errors, giving f þ ð0Þ=f 0 ð0Þ ¼ 0.95ð11Þ. We take c val s ¼ 0 AE 1 based on a study of s quark mistuning. In the case of maximum mistuning, where m s ¼ m l and we have the B → π decay, we can compare our form factors with those from [7] and find that this gives an upper bound on our valence quark mistuning of c val s ≈ 2. This is a very extreme case of quark mistuning, so we take the prior width at half of this. Sea quark mistunings, as well as those of the valence charm quark, make less of a contribution so we give all other c coefficients a prior of 0.0 AE 0.3. In Eq. (19) we take N ijk ¼ 3.
In our fit we also include a data point corresponding to the B s → η s scalar form factor in the continuum, f 0 ðq 2 max Þ ¼ 0.811ð17Þ from previous work by the HPQCD Collaboration [47]. This data point was obtained in a calculation using NRQCD b quarks, working directly at the tuned b quark mass. A ratio was constructed to remove the systematic errors from renormalization of the NRQCD-HISQ current that would otherwise reduce the accuracy of the result. For this reason, this point can be included alongside our HISQ data, without introducing additional errors associated with NRQCD. This result is included as a statistically independent data point for the f 0 fit function in the continuum and physical quark mass limits and reduces our error at f 0 ðq 2 max Þ. The effect of its removal is demonstrated by test 4 in Fig. 7.

Continuum and physical quark mass limit
The fit outlined in the previous section has a χ 2 value of 0.16 per degree of freedom, for 109 degrees of freedom. It produces best-fit results for the coefficients in Eqs. (19) and (20), from which we construct the z-expansion coefficients of Eq. (16).
By evaluating Eq. (19) at a; N 0;þ n ¼ 0, we obtain the z-expansion coefficients, and therefore the form factors from Eq. (16), in the continuum limit and at physical light, strange and charm quark masses. In Table V

Fit analysis and stability check
In Figs. 5 and 6 we show our lattice data in z space, as well as the results of the fit at the physical point for each form factor. In both cases these are plotted with their respective poles removed. We see very little z dependence in the f 0 case, which we can also infer from our a 0 1 and a 0 2 values (Table V), both of which are consistent with zero. In contrast, f þ displays a negative linear z dependence, again clear in the expansion coefficients. Both of these trends are similar to the findings in [16]. Both cases have large errors in some ultrafine data, which simply arises from lack of statistics on the very computationally expensive ultrafine configurations.
The lowest masses on each set correspond approximately to the physical charm mass, and we can see in Fig. 5 that these points lie on top of each other, indicating that lattice artifacts such as discretization errors are small at this mass. Other masses that are approximately equal are the ultrafine am h ¼ 0.45 and superfine am h ¼ 0.6, the ultrafine am h ¼ 0.6 and superfine am h ¼ 0.8, and the superfine am h ¼ 0.45 and fine am h ¼ 0.683. By comparing these values in Fig. 5 we can see that, while lattice artifacts become slightly more significant above the charm mass, they are still small, and that the heavy mass dependence itself is what dominates this plot. The picture is less clear in Fig. 6 because of larger statistical errors, but it appears to be similarly dominated by heavy quark mass dependence.
We verify our results for the form factors at the physical point are stable with respect to reasonable variations of the fit by modifying the fit as illustrated in Fig. 7 and discussed in the caption. The fit is stable under these variations, suggesting associated systematic uncertainties are adequately accounted for. This is plotted over the full q 2 range. We note that the error in the f 0 form factor shrinks with q 2 , while that in f þ grows. This is true even without the continuum data point TABLE V. Values of fit coefficients a 0;þ n and pole masses at the physical point for the B s → η s decay with correlation matrix are given below. Form factors can be reconstructed by evaluating Eq. (16) using these coefficients and pole masses. Note that M B s0 is set to M B s þ 0.4 GeV. Masses are in GeV. The pole masses are very slightly correlated due to the way the fit function is constructed. These correlations are too small to have any meaningful effect on the fit, but we include them for completeness in reconstructing our results.  respectively, where e.g., m0.8 indicates am h ¼ 0.8 on that ensemble. Lines between data points of a given heavy mass over the full z range are there to guide the eye. The additional continuum data point from [47] is shown as a purple diamond and helps to pin down the form factor in the high q 2 limit. from [47], so that statistical errors completely dominate f 0 ðq 2 max Þ. The vector form factor has a minimum error somewhere in between 0 and q 2 max , where our data are most densely distributed. This trend is similar in the scalar form factor if we remove the continuum data point which dominates the error at high q 2 . We also note that the quark mistuning and input errors for both cases are small and almost independent of q 2 , as we would expect. It is clear that our error is statistics dominated, which is a strong affirmation of the heavy HISQ method and nonperturbative current renormalization, as well as of the suitability of our z expansion. This also leaves the door open to a significant reduction in error, simply by increasing our statistics, particularly on the finest ensemble, a costly but straightforward exercise. We can see that, with sufficient computing time, errors could be reduced to 2%-3% across the full q 2 range for both the scalar and vector form factors.

III. FORM FACTOR RESULTS AND COMPARISONS
Our form factors at zero lattice spacing and physical quark mass are shown over the full physical q 2 range in Fig. 9. We can compare these with B s → η s results from a lattice calculation that used NRQCD b quarks given in the Appendix of [32]. We find the results to be in good agreement with an improvement in uncertainty across the q 2 range in the case of the f 0 form factor and an improvement by a factor of 2 at q 2 ¼ 0. The systematic uncertainties in the NRQCD calculation are dominated by the extrapolation to q 2 ¼ 0 from high q 2 values close to zero recoil and the associated discretization errors. The use of relatively coarse lattices in the NRQCD approach means that results are restricted to small daughter meson momentum. There is also a sizable systematic uncertainty from current renormalization present in the NRQCD results. We do not have these sources of error here. Our result for f þ ðq 2 max Þ agrees to 1σ with the NRQCD value, but with significantly larger uncertainty. This is a region of q 2 space Test 4 is the fit without the data point from [47]. Test 5 adds a cubic term in the z expansion [Eq. (16)]. Test 6 shows the effect of extending the i, j, k sum in Eq. (19). Tests 7 and 8 remove the highest masses and momenta for all lattice spacings, respectively. Test 9 is without the log term in Eq. (19); here we find that d i000 terms change to mimic the Taylor expansion of the log, and we require much larger priors (0 AE 5) to account for this. Test 10 shows the effect of doubling the width of all d ijkn priors. We see that our extrapolation is stable to all of the above modifications. Increasing the prior widths decreases the GBF, giving us confidence our priors are chosen conservatively.  [42] and used in the fit as described above. The purple dotted line ("q mistunings") adds, negligibly, to the inputs the error contribution from the quark mistunings associated with c fit parameters, while the solid green line ("statistics") further adds the error from our correlator fits. The blue dot-dashed line ("HQET") includes the contribution from the expansion in the heavy quark mass, and, finally, the thick black line ("Discretization"), the total error on the form factor, also includes the discretization errors. The percentage variance adds linearly and the scale for this is given on the lefthand axis. The percentage standard deviation, the square root of this, can be read from the scale on the right-hand side. where our data have large statistical errors because of the way that f þ is constructed from a temporal vector current in that limit. The differential rate for the decay vanishes rapidly toward q 2 max so it is the smaller values of q 2 at which we want to improve lattice QCD determination of the form factors and we have succeeded in doing this.

A. Comparisons testing SU(3) flavor and heavy quark symmetries
While the B s → η s decay does not correspond to a physical process, it is related to a host of physical decays via combinations of SU(3) flavor and heavy quark symmetry. In this section, we evaluate these symmetries by comparing to published results for symmetry-related decays. Figure 10 shows the effect of changing heavy quark mass over the full range of M H s from the physical M D s to the physical M B s , for both form factors at q 2 ¼ 0 [recall that f þ ð0Þ ¼ f 0 ð0Þ] and at maximum physical q 2 . Our use of a range of heavy masses from the physical charm to the physical bottom allows for good control of this heavy mass dependence. The uncertainty at the lighter end is particularly small, as all three ensembles had a physical charm mass data point, whereas only set 3 was fine enough to give data at the physical bottom mass. f 0;þ ð0Þ, f 0 ðq 2 max Þ and f þ ðq 2 max Þ are converging as M H s is reduced and one can imagine them meeting if extrapolated in mass below M D s to M η s . That point would correspond to the η s → η s decay, where only q 2 ¼ 0 is kinematically allowed and we expect f þ ¼ f 0 ¼ 1. A similar effect was seen in [16].
Previous lattice QCD results for other decay processes related by SU(3) flavor symmetry are included in Fig. 10 in the same color labeling system. We see very good agreement with the D → K and B → K decays for both form factors at both ends of the q 2 range, suggesting that the mass of the spectator quark has almost no effect on the form factors and supporting our use of B s → η s to test the viability of a B → K calculation. B s → K data show good agreement for f 0 but f þ ðq 2 max Þ is in slight tension. This suggests, as expected, that the form factors are much more sensitive to SU(3) flavor symmetry breaking in the daughter quark in the transition than in the spectator quark. This is further supported by the D → π results, which are in poor agreement with our D s → η s form factors across the board. B → π results are in even worse agreement and are not included in the plot. This implies that symmetry breaking in the light daughter quark becomes even more important as the heavy parent quark becomes heavier.

B. Tests of HQET
That we are able to evaluate our form factors over the full range m c ≤ m h ≤ m b means we are in a unique position to test predictions of HQET. One such set of predictions relates to the characterization of form factor shape. The quantities α, δ and β −1 are used to describe the shape of the form factors in HQET [49,50]. The latter two of these are related to the slope of the form factors at q 2 ¼ 0 and the first to the value at high q 2 : : ð25Þ Figure 11 shows our results for these quantities, plotted across the full range of heavy masses from c to b using as the x axis the mass of the heavy-strange pseudoscalar meson. Our results for α and β are qualitatively in agreement with expectations from HQET [49] with α and β close to one at the heaviest masses and differing further from one as the heavy quark mass falls. Our results are accurate enough that they could be used to constrain scaling laws in the mass from other theoretical approaches. We see that δ is close to zero at the B s end of the plot but clearly nonzero at the D s end. We find values of α M B s ¼ 0.698ð56Þ, The form factor ratio f 0 Fig. 12, where it is compared with the HQET expectation [51] lim This is included in [51] as a B → π expectation; to test it here in B s → η s we replace B with B s . We take the ratio of decay [52]. No difference is visible in this ratio between B s and B in [52]. We take the coupling g B Ã s B s η s ≈ g B Ã Bπ ¼ 0.56ð8Þ [53], because again the light quark mass dependence seen in [53] is mild. This leads us to expect little impact from SU(3) flavor symmetry breaking in our test of Eq. (26). This is also consistent with our observation in Fig. 10 that SU(3) flavor symmetry breaking effects in the daughter quark affect both f 0 and f þ at large q 2 , and so there will be some cancellation of the effects in their ratio. Figure 12 shows reasonable agreement with Eq. (26) in the limit q 2 → M 2 B s , as is found for B → π in [7]. Figure 13 tests the relationships between form factors for a changing initial state but fixed final state with a fixed energy. In [49] it is shown that the f 0 form factor for a pseudoscalar heavy meson decay to a pseudoscalar light meson at fixed energy is inversely proportional to the square root of the heavy meson mass. This scaling should work both at small energy, close to zero recoil, and also at large energy, high recoil. In [49] this is used to compare B → π and D → π decay. Here we compare B s → η s to H s → η s for variable H s mass from D s upward.
, as compared with the HQET expectation in the limit q 2 → M 2 B s (red band), defined in Eq. (26).
Our results at both energies are flatter than the ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1=M H s p expectation, indicating that sizable corrections are needed to this expectation to describe the physical behavior. This is reminiscent of results for the decay constant of heavystrange pseudoscalar mesons in that it does not vary so strongly with mass as predicted; [20] shows that this decay constant only changes by 9.4(1.4)% over the range from c to b when the leading-order HQET behavior is as ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1=M H s p , i.e., a 65% change.

IV. CONCLUSIONS
We have performed the first calculation of form factors for a b → light quark transition in which we use our heavy-HISQ technique. This requires results at multiple values of the heavy quark mass on multiple sets of gluon field configurations with fine lattice spacing (going down to 0.045 fm here) so that we can map out the heavy quark mass dependence of the form factors and obtain physical results for a heavy quark mass equal to that of the b. One advantage of this technique over previous calculations is that we can normalize the lattice currents completely nonperturbatively. Here we do this for the vector and scalar currents that give the vector and scalar form factors. This means that we can avoid sizable systematic errors from the one-loop matching of lattice currents to continuum currents that is done, for example, for NRQCD b quarks. A second advantage of the heavy-HISQ technique is that it enables us to cover the full range in q 2 of the decay rather than just values of q 2 close to zero recoil (low momentum for the daughter meson). This is possible because the accessible range in q 2 grows as the accessible range in heavy quark mass grows on finer lattices.
As a stepping-stone toward a variety of physical decay processes, we have chosen to study first the unphysical process B s → η s here because this does not involve valence u or d quarks and the s quark mass can be accurately tuned to its physical value on all of our gluon field configurations. We present our final form factor results in Fig. 9. The form factor values at the end points of the q 2 range are Our uncertainty for the form factor at the kinematically important point (for the differential rate) q 2 ¼ 0 is 8%. This is an improvement by a factor of 2 over earlier results that used NRQCD b quarks and coarser lattices. The uncertainties of the NRQCD result were dominated by the extrapolation of lattice results from relatively high q 2 values to q 2 ¼ 0, along with the associated discretization effects, statistical errors and a current matching uncertainty of 3%. Our error budget as a function of q 2 is given in Fig. 8  and is dominated by statistical errors that can be improved at the cost of additional computing resource, to 2%-3% over the full q 2 range.
Although our results correspond to an unphysical process, B s → η s is related to physical processes through SU(3) flavor symmetry for the light quark. Because we have results for the range of heavy quark masses from c to b we can study this SU(3) symmetry breaking through comparison to previous lattice QCD results for the physical processes for both B and D decay. This is shown in Fig. 10. We find that SU(3) flavor symmetry breaking in the daughter quark in the transition affects the form factors increasingly as the parent quark gets lighter. In contrast, symmetry breaking in the spectator quark has very little effect.
HQET expectations for the mass scaling behavior of form factors for h → l decay should hold for B s → η s up to effects from the s quark mass, which should be small. We show comparison to such expectations in Figs. 12-14. The latter two show substantial corrections to the leading-order HQET behavior are present.
Our results provide further evidence that the heavy-HISQ approach is an improved method for calculating hadronic form factors for semileptonic decays involving heavy quarks. This leads us to conclude that a heavy HISQ calculation of form factors for a physical b → s process, B → Kl þ l − will be able to improve upon the previous errors in [10,56]. An accurate determination of the renormalization of the lattice tensor current [57], possible with HISQ quarks, will allow us to improve the determination of the tensor form factor for that process as well. Our results are also encouraging for similar calculations involving b → l decays, such as B → π and B s → K, enabling improvement in the determination of the CKM element V ub when combined with experimental results.

ACKNOWLEDGMENTS
We are grateful to the MILC Collaboration for the use of their configurations and their code. We would also like to thank J. Harrison (24) 0.54 (35)