B_{s}\to K\ell\nu form factors from lattice QCD

form factors from lattice QCD

C.M. Bouchard Department of Physics, The Ohio State University, Columbus, Ohio 43210, USA    G. Peter Lepage Laboratory of Elementary Particle Physics, Cornell University, Ithaca, New York 14853, USA    Christopher Monahan Physics Department, College of William and Mary, Williamsburg, Virginia 23187, USA    Heechang Na Department Physics and Astronomy, University of Utah, Salt Lake City, Utah 84112, USA    Junko Shigemitsu Department of Physics, The Ohio State University, Columbus, Ohio 43210, USA
July 12, 2019

We report the first lattice QCD calculation of the form factors for the standard model tree-level decay . In combination with future measurement, this calculation will provide an alternative exclusive semileptonic determination of . We compare our results with previous model calculations, make predictions for differential decay rates and branching fractions, and predict the ratio of differential branching fractions between and . We also present standard model predictions for differential decay rate forward-backward asymmetries and polarization fractions and calculate potentially useful ratios of form factors with those of the fictitious decay. Our lattice simulations utilize nonrelativistic QCD and highly improved staggered light quarks on a subset of the MILC Collaboration’s asqtad gauge configurations, including two lattice spacings and a range of light quark masses.

12.38.Gc, 13.20.He, 14.40.Nd, 14.40.Df
thanks: bouchard.18@osu.edu

HPQCD Collaboration

I Introduction

The decay occurs at tree-level in the standard model via the flavor-changing charged-current transition, making it an alternative to in the determination of from exclusive semileptonic decays. The difference in these processes, a spectator strange quark in vs a spectator down quark in , is beneficial for lattice QCD simulations, because it improves the ratio of signal to noise. Though this process has not yet been observed, its measurement is planned at LHCb and is possible during an run at BelleII. This provides a prediction opportunity for lattice QCD.

In addition to the calculation of form factors for , we also calculate their ratios with form factors for the fictitious decay. Such ratios are essentially free of our largest systematic error, perturbative matching. In combination with a future calculation of using a highly improved staggered (HISQ) quark, these ratios would yield a non-perturbative evaluation of the matching factor for the current with nonrelativistic QCD (NRQCD) quark. This matching factor would be applicable to and simulations using NRQCD quarks.

To include correlations among the data for both decays, correlation function fits must include vast amounts of correlated data. To make such fits feasible, we have developed a new technique, called chaining, discussed in Appendix A. In addition, the use of marginalization techniques developed in Ref. Hornbostel:2011 () significantly reduces the time required for the fits.

The chiral, continuum, and kinematic extrapolations are performed simultaneously using the modified  expansion Na:2010 (); Na:2011 () with the chiral logarithmic corrections fixed by the results of hard pion chiral perturbation theory (HPChPT) Bijnens:2010 (); Bijnens:2011 (). The factorization of chiral corrections and kinematics, as found at one-loop order by HPChPT, suggests the modified  expansion is a natural choice for carrying out this simultaneous extrapolation. We refer to the combination of HPChPT chiral logarithmic corrections and the modified  expansion as the HPChPT  expansion.

Ii Form factors and matrix elements

The vector hadronic matrix element is parametrized by the scalar and vector form factors


where and . At intermediate stages of the calculation we recast in terms of the more convenient form factors


where . In the meson rest frame, the form factors are simply related to the temporal and spatial components of the hadronic vector matrix elements,


The scalar and vector form factors are related to by


where is the kaon three-momentum. This discussion generalizes in a straightforward way for the matrix elements.

Iii Simulation

C1 2.647(3) 0.005/0.05 0.8678 1200 2 0.0070 0.0489 2.650 0.28356(15)
C2 2.618(3) 0.01/0.05 0.8677 1200 2 0.0123 0.0492 2.688 0.28323(18)
C3 2.644(3) 0.02/0.05 0.8688 600 2 0.0246 0.0491 2.650 0.27897(20)
F1 3.699(3) 0.0062/0.031 0.8782 1200 4 0.00674 0.0337 1.832 0.25653(14)
F2 3.712(4) 0.0124/0.031 0.8788 600 4 0.01350 0.0336 1.826 0.25558(28)
Table 1: Left to right: labels for the ensembles used in this analysis; lattice volume; inverse lattice spacing in -units; light/strange sea-quark masses; tadpole improvement factor ; number of configurations; number of time sources; valence -quark mass; valence -quark mass; -quark mass; and the spin-averaged ground state energies used to relate our meson simulation energies to their physical values.

Ensemble averages are performed with the MILC Collaboration’s asqtad gauge configurations Bazavov:2010 () listed in Table 1. Valence quarks in our simulation are nonrelativistic QCD (NRQCD) Lepage:1992 () quarks, tuned in Ref. Na:2012 (), and highly improved staggered (HISQ) Follana:2007 () light and quarks, the propagators for which were generated in Refs. Na:2010 (); Na:2011 (). Valence quark masses for each ensemble used in the simulations are collected in Table 1 and correspond to pion masses ranging from, approximately, 260 MeV to 500 MeV.

Heavy-light meson bilinears are built from NRQCD and HISQ quarks (for details see Ref. Na:2012 ()) and light-light kaon (and similarly for the ) bilinears are built from HISQ light and quarks (for details see Ref. Na:2010 ()). From these bilinears we build two and three point correlation function data


where indices specify quark smearing. We generate data for both a local and Gaussian smeared quark, with smearing function introduced via the replacement in Eqs. (7) and (9). Three point and daughter meson two point correlation function data are generated at four daughter meson momenta, corresponding to . In three point data, these momenta are inserted at in Fig. 1. The sum over in Eqs. (8) and (9) is performed using random wall sources with U(1) phases , i.e. . In the three point correlator a meson source is inserted at timeslice , selected at random on each configuration to reduce autocorrelations. The current is inserted at timeslices such that and the daughter meson is annihilated at timeslice . Prior to performing the fits, all data are shifted to a common . This three point correlator setup is depicted in Fig. 1. Additional details regarding the two and three point correlation function generation can be found in Ref. Bouchard:2013a ().

The flavor-changing current is an effective lattice vector current corrected through . The lattice currents that contribute through this order are


Matrix elements of the continuum vector current are matched to those of the lattice vector current according to




The matching calculation is done to one loop using massless HISQ lattice perturbation theory Monahan:2013 (). In implementing the matching, we omit contributions. Ref. Gulez:2007 (), which used asqtad valence quarks, found contributions of this order to be negligible. In Ref. Bouchard:2013a (), which used HISQ valence quarks, these contributions to the temporal component of the vector current were studied and again were found to be negligible. We also omit relativistic matching corrections. These, and higher order, omitted contributions to the matching result in our leading systematic error. An estimate of this error, and its incorporation in our fit results, is discussed in the following section.

Figure 1: Setup for three point correlator data generation.

Iv Correlation function fits

Two and three point correlation function fit Ansätze, and the selection of priors, closely follows the methods of Ref. Bouchard:2013a (). Two point data are fit to


where tildes denote oscillating state contributions and is the simulated energy. The physical ground state mass is related to the simulation ground state energy by


where GeV Gregory:2011 () is adjusted from experiment to remove electromagnetic, annihilation, and charmed sea effects not present in our simulations, and is the spin-averaged energy of states calculated on the ensembles used in the simulation and listed in Table 1. The quark smearing is indicated by indices . Kaon and two point correlator data are fit to an expression of the form111The zero momentum has no oscillating state contributions due to mass degeneracy of its valence quarks.


Results of two point fits satisfy the dispersion relation and are stable with respect to variations in and the range of timeslices included in the fits, as demonstrated for kaon two point data in Ref. Bouchard:2013a ().

Three point correlation function data are described by


where the three point amplitudes , , , and are proportional to the hadronic matrix elements. The ground state hadronic matrix element is obtained from


where the factor of accounts for numerical factors introduced in the simulation and associated with taste averaging and HISQ inversion. In the correlator fits we include data for several temporal separations between the mother and daughter mesons. On the coarse ensembles we include data for while for the fine ensembles we include data.

On each ensemble we perform a simultaneous fit to two and three point correlation function data for the and decays, at all simulated momenta, including both spatial and temporal currents, and for the temporal separations listed above. This ensures correlations among these data are accounted for in the analysis. However, fits to such large data sets produce unwieldy data covariance matrices and are typically not convergent, or require a prohibitively large number of iterations. This can be partially addressed by thinning the data, e.g. by the use of singular value decomposition (SVD) cuts, but this reduces the accuracy of the fits.

Figure 2: (color online). Chained and marginalized fit results for the ground state amplitude of the decay on ensemble F2. Fit results are shown as a function of the number of (top) states explicitly included in the fit and (bottom) total states accounted for in the fit. Final results are taken from (6,1)/(8,8) fits, represented by gray bands.

To address this problem we introduce a technique, which we refer to as chaining, to simplify fits to very large data sets. Consider a data set consisting of correlators, . Before the fit, all fit parameters are assigned priors. Chaining first fits then uses the best fit mean values and covariances to replace the corresponding priors in subsequent fits. The updated set of priors is then used in the fit to . In this and all subsequent fits, correlations are accounted for between the data being fit and those priors which are best fit results from previous fits — this is an important step as it prevents “double counting” data. After this second fit, the priors are again updated according to the best fit mean values and covariances. This process is repeated for all correlators. The collection of best fit mean values and covariances following the fit to are the final fit results. Chaining is described in greater detail in Appendix A.

We combine the use of Bayesian Lepage:2002 (), marginalized Hornbostel:2011 (), and chained fitting techniques. Our final fit results use marginalization with a total of states accounted for, of which are explicitly fit. We refer to such fits with the shorthand notation, . States accounted for but not explicitly fit are marginalized in that their contributions are subtracted from the data prior to the fit. This technique reduces significantly the time required to perform the fits. In Fig. 2 we show the stability of the fits under variations in the numbers of states explicitly included and the total number of states accounted for in the fit.

C1 0.8244(23) 0.7081(27) 0.6383(30) 0.5938(41)
C2 0.8427(25) 0.6927(35) 0.6036(49) 0.536(12)
C3 0.8313(29) 0.6953(33) 0.6309(30) 0.5844(46)
F1 0.8322(25) 0.6844(35) 0.5994(43) 0.5551(56)
F2 0.8316(27) 0.6915(38) 0.6119(43) 0.5563(61)
C1 2.087(16) 1.657(14) 1.378(13)
C2 1.880(12) 1.412(16) 1.142(33)
C3 1.773(11) 1.4212(84) 1.184(10)
F1 1.878(13) 1.385(12) 1.158(13)
F2 1.834(14) 1.396(10) 1.163(14)
Table 2: Fit results for the scalar and vector form factors on each ensemble and for each simulated momentum.

The form factor results from the correlation function fits are tabulated in Table 2 and additional details are given in Appendix B.

Figure 3: (color online). (top) Heat map of the correlation matrix for ensemble C1. (bottom) Distribution of correlations among the form factors for and for all ensembles.

The form factors obtained from these fits preserve correlations resulting from shared gauge field configurations and quark propagators used in data generation. The preservation of correlations is demonstrated in the top panel of Fig. 3 where, e.g., significant correlations among the form factor fit results are seen at common momenta and nonzero correlations among form factors for the two decays is suggested. The bottom panel of Fig. 3 shows the distribution over all ensembles of correlations among form factors for the two decays. Accounting for these correlations is useful in our determination of the ratio of form factors for the two decays. Fit results for , and the resulting form factor ratios, are presented in Appendix D.

V Chiral, continuum, and kinematic extrapolation

The results of HPChPT Bijnens:2010 (); Bijnens:2011 () suggest a factorization, to at least one-loop order, of the soft physics of logarithmic chiral corrections and the physics associated with kinematics in the form factors describing semileptonic decays of heavy mesons,


The logarithmic chiral corrections, calculated in Ref. Bijnens:2011 () for several decays, are independent of . An unspecified function characterizes the kinematics.

To obtain results over the full kinematic range one must include lattice simulation data over a range of energies. However, for any relevant physical scale (e.g. , ), at nominal lattice momenta and there is no convergent expansion of the unknown function in powers of . This is an inherent limitation of characterizing the kinematics in terms of energy. The energy of the daughter meson is a poor variable with which to describe the kinematics.

In contrast, the  expansion Boyd:1996 (); Arnesen:2005 (); Bourrely:2010 () provides a convergent, model-independent characterization of the kinematics over the entire kinematically accessible range. Combining a  expansion on each ensemble222This assumes the general arguments on which the  expansion is based hold for heavier than physical quark masses and at finite lattice spacing. with the HPChPT inspired factorization of Eq. (19) allows a simultaneous chiral, continuum, and kinematic extrapolation of lattice data at arbitrary energies. Because the chiral logs are the same for and , linear combinations (i.e. and ) factorize in the same way and have the same chiral logs. Motivated by these observations, we construct a HPChPT-motivated modified  expansion, which we call the “HPChPT  expansion”, and fit the lattice data of Tables 2 and 9, with accompanying covariance matrix, to fit functions of the form


where [logs] are the continuum HPChPT logs of Ref. Bijnens:2011 (), and generic analytic chiral and discretization effects are accounted for by . Resonances above but below the production threshold, i.e. those in the range , are accounted for via the Blaschke factor, . Though not observed, we allow for the possibility of a state in , with choice of mass guided by Ref. Gregory:2011 (). Our fit results are insensitive to the presence of this state. The factorization suggested by HPChPT may not hold at higher order Colangelo:2012 () so we allow chiral analytic terms, which help parametrize effects from omitted higher order chiral logs, to have energy dependence (i.e. to vary with ).

We note that Eq. (20) is the modified  expansion introduced in Refs. Na:2010 (); Na:2011 (), with the coefficients of the chiral logarithmic corrections fixed by the results of HPChPT. In the chiral and continuum limits


and Eq. (20) is equivalent to the Bourrely-Caprini-Lellouch parametrization Bourrely:2010 () of the form factors.

Following Ref. Bourrely:2010 () we impose a constraint on from the expected scaling behavior of in the neighborhood of . The resulting fit function for is


We write , , and , explicitly exposing the dependence on and . This is useful in explaining the implementation of a second kinematic constraint we impose on the form factors. At the kinematic endpoint , the continuum extrapolated form factors and are equal, i.e. . We impose this constraint by fixing the coefficient ,


Imposing this constraint results in the fit function for :


In the fit functions for and , Eqs. (22) and (24), and [logs] are given by,


with implicit indices in Eq. (25) specifying the scalar or vector form factor. We account for momentum-independent and momentum-dependent discretization effects in . The values of that enter the fit are the values from the simulation and are, of course, small. Finite volume effects in the simulation are included via a shift in the pion log Bernard:2002 (). The infinite volume limit is taken by setting this shift to zero. Eq. (26) gives the HPChPT Bijnens:2011 () result for the chiral logarithmic correction to form factors. These expressions make use of the dimensionless quantities


where . We determine and on each ensemble using correlator fit results for meson masses and simulation momenta. Light and heavy quark discretization effects are accommodated for by making the mild functions of the masses, accomplished by the replacements

where is chosen so that as varies over the coarse and fine ensembles .

Lastly, we account for uncertainty associated with the perturbative matching of Sec. III. With the matching coefficients calculated in Ref. Monahan:2013 (), we find contributions to be of the total contribution to . Of this the majority, , comes from the one loop correction and from the NRQCD matching via . For we find contributions at this order to be , with coming from the correction and from the NRQCD matching. The matching error results from omitted higher order corrections, the size of which we estimate from observed leading order effects, where we conservatively use the larger 4%. Following the arguments outlined in Ref. Bouchard:2013a () we estimate the matching error to be the same size as the observed contributions and take the matching error to be 4%. This is equivalent to taking the matching coefficient to be four times larger than the matching coefficient (13 times larger than ). This uncertainty is associated with the hadronic matrix elements and therefore, by Eqs. (3) and (4), with and . To correctly incorporate it in the results for and we convert our fit functions for into , multiply by , where is a coefficient representing the matching error with a prior central value of zero and width 0.04, then convert back to before performing the fit. Schematically, we modify the fit functions, defined in Eqs. (22) and (24), by


then we use to fit the results of the correlation function fits of Sec. IV. Conversions between the form factors and are performed using Eqs. (5) and (6).

Figure 4: (color online). form factor results from a simultaneous chiral, continuum, and kinematic extrapolation via the HPChPT  expansion are shown (top) relative to coarse ensemble data (C1, C2, and C3) and (bottom) relative to fine ensemble data (F1 and F2).

The results of a simultaneous fit to the data for and , in which the maximum order of [specified by in Eqs. (22) and (24)] is 3 and , are shown relative to the data in Fig. 4 for . Details of prior choices and fit results are given in Appendix C.

We test the stability of this fit to the following modifications of the fit Ansätze:

  1. Truncate the  expansion at .

  2. Truncate the  expansion at .

  3. Truncate the  expansion at .

  4. Drop momentum-dependent and momentum-independent discretization terms in Eq. (25).

  5. Drop the -dependent discretization terms in Eq. (LABEL:eq-disc).

  6. Drop the light-quark mass-dependent discretization terms in Eq. (LABEL:eq-disc).

  7. Add the following next-to-next-to-leading-order (NNLO) chiral analytic terms to as defined in Eq. (25):

  8. Drop the sea- and valence-quark mass difference term from Eq. (25).

  9. Drop the strange quark mistuning term from Eq. (25).

  10. Drop finite volume effects, i.e. set in Eq. (26).

The stability of the fit results to these modifications is shown in Fig. 5, where results are shown at the extrapolated point. This point is furthest from the data region where simulations are performed and therefore is particularly sensitive to changes in the fit function. In Fig. 5 our final fit result, as defined by Eqs. (22) and (24) with and by Eqs. (25)–(LABEL:eq-disc), is indicated by the dashed line and gray band.

Modifications 1, 2, and 3 vary the order of the truncation in and demonstrate that by fit results have stabilized and errors have saturated. We therefore conclude that the error of the fit adequately accounts for the systematic error due to truncating the  expansion.

Momentum-dependent and momentum-independent discretization effects proportional to are removed in modification 4. This results in a modest increase in and a negligible shift in the fit result. This suggests our final fit, which includes the effects, adequately accounts for all discretization effects observed in the data.

In modifications 5 and 6 we remove heavy- and light-quark mass-dependent discretization effects with essentially no impact on the fit. That our results are independent of light-quark mass dependent discretization effects suggests that staggered taste violating effects are accommodated for by a generic dependence.

Modification 7 tests the truncation of chiral analytic terms after next-to-leading-order (NLO) by adding the NNLO terms listed in Eq. (34). This results in a slight decrease in but has no noticeable effect on the fit central value or error. From this we conclude that errors associated with omitted higher order chiral terms are negligible.

Differences in sea and valence quark masses, due in part to our use of HISQ valence- and asqtad sea-quarks, are neglected in modification 8. This results in a small increase in and negligible change in the fit results. We account for these small mass differences in our final fit, though this test suggests they are unimportant in the fit.

Effects due to strange quark mass mistuning on the ensembles are omitted in modification 9, resulting in a modest increase in and no change in the fit central value and error. We include these effects in our final fit.

Modification 10 results in nearly identical fit results, suggesting that finite volume effects are negligible in our data. We include these effects in our final fit results.

Figure 5: (color online). The stability of the HPChPT  expansion is demonstrated by studying the fit results under various modifications, discussed in the text. The top panel shows with 70 degrees of freedom (d.o.f.) for each test fit and the bottom panel shows form factors extrapolated to .

Vi Form Factor Results

Figure 6: (color online). form factor results from a simultaneous chiral, continuum, and kinematic extrapolation via the HPChPT  expansion. The region for which lattice simulation data exist is indicated by the shaded region.

In this section we present final results, with a complete error budget, for the form factors. We provide the needed information to reconstruct the form factors and compare our results with previous model calculations.

Fig. 6 shows the results of the chiral, continuum, and kinematic extrapolation of Sec. V, plotted over the entire kinematic range of . The form factors, extrapolated to , have the value .

vi.1 Fit errors for the HPChPT  expansion

The inputs in our chiral, continuum, and kinematic extrapolation fits are data (the correlator fit results for and in Tables 2 and 9 with the accompanying covariance matrix) and priors. The total hessian error of the fit can be described in terms of contributions from these inputs, as described in detail in Appendix A. We group priors in a meaningful, though not unique, way and discuss the error associated with the chiral, continuum, and kinematic extrapolation based on these groupings. As the priors are, by construction, uncorrelated with one another, we can group them together in any way we find meaningful. The resulting error groupings are uncorrelated and add in quadrature to the total error. In Fig. 7 we plot the following relative error components as functions of :

  1. experiment: This is the error in the fit due to uncertainty of experimentally determined, and other, input parameters. It is the sum in quadrature of the errors due to priors for the “Group I” fit parameters listed in Table 7. This error is independent of and subdominant.

  2. kinematic: This error component is due to the priors for the coefficients in Eqs. (22) and (24). A comparison of the fit results from modifications 1, 2, and 3 in Fig. 5 shows that by the fit results have stabilized and errors have saturated. The kinematic error therefore includes the error associated with truncating the  expansion. The extrapolation to values of for which we have no simulation data is controlled by the  expansion. As a result, the growth in form factor errors away from the simulation region is due almost entirely to kinematic and statistical errors.

  3. chiral: This error component is the sum in quadrature of errors associated with priors for in Eq. (25). These terms are responsible for extrapolating to the physical light quark mass and for accommodating for the slight strange quark mistuning and the small mismatch in sea and valence quark masses due to the mixed action used in the simulation. As shown in Fig. 7, these errors are subdominant and do not vary significantly with .

  4. discretization: We account for momentum-dependent discretization effects via the , and momentum-independent discretization effects via the , terms of Eq. (25). In addition we allow for heavy- and light-quark mass-dependent discretization effects via the and terms in Eq. (LABEL:eq-disc). The discretization error component, which is essentially independent of , is the sum in quadrature of the error due to the priors for these fit parameters.

  5. statistical: The statistical component of the error is due to uncertainty in the data, i.e. the errors from form factor fit results of Table 2. Simulation data exist for for and over the range for . Extrapolation beyond these regions leads to increasing errors.

  6. matching: The matching error is due to the uncertainty associated with the priors for introduced in Eq. (32) and discussed in the surrounding text.

Figure 7: (color online). (top) and (bottom) relative error components. The total error (solid line) is the sum in quadrature of the components.

In addition to the largest sources of error, which we account for directly in the fit, there are remaining systematic uncertainties.

We simulate with degenerate light quarks and neglect electromagnetism. By adjusting the physical kaon mass () used in the chiral, continuum and kinematic extrapolation, we estimate the “kinematic” effects of omitting electromagnetic and isospin symmetry breaking in our simulation to be . It is more difficult to determine the size of the full effects. However, in general electromagnetic and isospin effects are expected to be sub-percent. We assume the error in our form factor calculation due to these effects is negligible relative to other sources of uncertainty.

Our simulations include up, down, and strange sea quarks and we assume omitted charm sea quark effects are negligible. This has been the case for processes in which it has been possible and appropriate to perturbatively estimate effects of charm quarks in the sea Davies:2010 ().

Our final form factor results, multiplied by the Blaschke factor , are shown in Fig. 8 where they are compared with results from a model calculation using perturbative QCD (pQCD) Wang:2012 () and a relativistic quark model (RQM) Faustov:2013 (). Our results provide significant clarification on the form factors at large .

Figure 8: (color online). Comparison of our (top) and (bottom) form factors with those from a pQCD model Wang:2012 () and the RQM Faustov:2013 (). To ease comparison the vertical scale is reduced by multiplying the form factors by the Blaschke factor .

vi.2 Reconstructing Form Factors

Coefficient Value
Table 3: (Top) Physical extrapolated coefficients of the HPChPT  expansion for the form factors, defined in Eqs. (35) and (36) and (bottom) the associated covariance matrix.

In the physical limit our form factor results are parametrized in a BCL Bourrely:2010 () form with coefficients [see Eq. (21)]. Including the kinematic constraint and terms through order , we have




and the resonance masses are and . The values of the coefficients , derived from the extrapolation fit results of Sec. V, and the associated covariance matrix, are given in Table 3. Note that it is necessary to take into account the correlations among the coefficients to correctly reproduce the form factor errors.

Vii Phenomenology

With the benefit of ab initio form factors from lattice QCD, we explore the standard model implications of our results. In this section we make standard model predictions for several observables related to the decay for and .

The standard model differential decay rate is related to the form factors by


In Fig. 9 we plot predicted differential decay rates for and , divided by , over the full kinematic range of .