Update of |V_{cb}| from the \bar{B}\to D^{*}\ell\bar{\nu} form factor at zero recoil with three-flavor lattice QCD

Update of from the form factor at zero recoil with three-flavor lattice QCD

Jon A. Bailey Department of Physics and Astronomy, Seoul National University, Seoul, South Korea    A. Bazavov Physics Department, Brookhaven National Laboratory, Upton, New York, USA    C. Bernard Department of Physics, Washington University, St. Louis, Missouri, USA    C. M. Bouchard Department of Physics, The Ohio State University, Columbus, Ohio, USA    C. DeTar Department Physics and Astronomy, University of Utah, Salt Lake City, Utah, USA    Daping Du Physics Department, University of Illinois, Urbana, Illinois, USA Department of Physics, Syracuse University, Syracuse, New York, USA    A. X. El-Khadra Physics Department, University of Illinois, Urbana, Illinois, USA Fermi National Accelerator Laboratory, Batavia, Illinois, USA    J. Foley Department Physics and Astronomy, University of Utah, Salt Lake City, Utah, USA    E. D. Freeland Department of Physics, Benedictine University, Lisle, Illinois, USA Liberal Arts Department, School of the Art Institute of Chicago, Chicago, Illinois, USA    E. Gámiz CAFPE and Departamento de Física Teórica y del Cosmos, Universidad de Granada, Granada, Spain    Steven Gottlieb Department of Physics, Indiana University, Bloomington, Indiana, USA    U. M. Heller American Physical Society, Ridge, New York, USA    A. S. Kronfeld ask@fnal.gov Fermi National Accelerator Laboratory, Batavia, Illinois, USA    J. Laiho jlaiho@fnal.gov SUPA, Department of Physics and Astronomy, University of Glasgow, Glasgow, United Kingdom Department of Physics, Syracuse University, Syracuse, New York, USA    L. Levkova Department Physics and Astronomy, University of Utah, Salt Lake City, Utah, USA    P. B. Mackenzie Fermi National Accelerator Laboratory, Batavia, Illinois, USA    E. T. Neil Fermi National Accelerator Laboratory, Batavia, Illinois, USA Department of Physics, University of Colorado, Boulder, Colorado, USA RIKEN-BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA    Si-Wei Qiu Department Physics and Astronomy, University of Utah, Salt Lake City, Utah, USA    J. Simone Fermi National Accelerator Laboratory, Batavia, Illinois, USA    R. Sugar Department of Physics, University of California, Santa Barbara, California, USA    D. Toussaint Department of Physics, University of Arizona, Tucson, Arizona, USA    R. S. Van de Water Fermi National Accelerator Laboratory, Batavia, Illinois, USA    Ran Zhou Department of Physics, Indiana University, Bloomington, Indiana, USA Fermi National Accelerator Laboratory, Batavia, Illinois, USA
July 16, 2019
Abstract

We compute the zero-recoil form factor for the semileptonic decay (and modes related by isospin and charge conjugation) using lattice QCD with three flavors of sea quarks. We use an improved staggered action for the light valence and sea quarks (the MILC asqtad configurations), and the Fermilab action for the heavy quarks. Our calculations incorporate higher statistics, finer lattice spacings, and lighter quark masses than our 2008 work. As a byproduct of tuning the new data set, we obtain the and hyperfine splittings with few-MeV accuracy. For the zero-recoil form factor, we obtain , where the first error is statistical and the second is the sum in quadrature of all systematic errors. With the latest HFAG average of experimental results and a cautious treatment of QED effects, we find . The QCD error is now commensurate with the experimental error.

pacs:
12.38.Gc, 13.20.He, 12.15.Hh

Fermilab Lattice and MILC Collaborations

I Introduction

The Cabibbo-Kobayashi-Maskawa (CKM) matrix element is one of the fundamental parameters of the Standard Model (SM). Together with , , and , it allows for a full SM determination of flavor and violation via processes that proceed at the tree level of the electroweak interaction. In the case of , one requires a measurement of the differential rate of mesons decaying semileptonically to a charmed final state. The hadronic part of the final state can be exclusive—e.g., a or meson—or inclusive.

The 2012 edition of the Review of Particle Physics by the Particle Data Group (PDG) Beringer et al. (2012) notes that the exclusive and inclusive values of are marginally consistent with each other. Furthermore, global fits to a comprehensive range of flavor- and -violating observables tend to prefer the inclusive value Laiho et al. (2010); Lunghi and Soni (2010); Lenz et al. (2011): when direct information on is omitted from the fit, one of the outputs of the fit is a value of that agrees better with the inclusive than the exclusive value. One should bear in mind that some tension in the global fits to the whole CKM paradigm has been seen Lunghi and Soni (2011). A full discussion of the possible resolutions of the discrepancy lies beyond the scope of this article. We conclude merely that it is important and timely to revisit the theoretical and experimental ingredients of both determinations.

In this paper, we improve the lattice-QCD calculation Hashimoto et al. (2002); Bernard et al. (2009); Bailey et al. (2010) of the zero-recoil form factor for the exclusive decay (and isopin-partner and charge-conjugate modes). Our analysis strategy is very similar to our previous work Bernard et al. (2009), but the lattice-QCD data set is much more extensive, with higher statistics on all ensembles, smaller lattice spacings (as small as  fm) and light-quark masses as small as (at lattice spacing  fm). Figure 1 provides a simple overview of the new and old data sets;

Figure 1: (color online) Range of lattice spacings and light-quark masses used here (colored or gray discs) and in Ref. Bernard et al. (2009) (black circles). The area is proportional to the size of the ensemble. The lattice spacings are , 0.12, 0.09, 0.06, and 0.045 fm. Reference Bailey et al. (2010) did not yet include the ensembles with , , , and .

further details are given in Sec. II. Our preliminary status report Bailey et al. (2010) encompassed the higher statistics but not yet four of the ensembles in the lower left-hand corner of Fig. 1.

With this work, we improve the precision of as determined from exclusive decays to that claimed for the determination from inclusive decays: 2%. Moreover, we reduce the QCD uncertainty on to the same level as the experimental uncertainty. Because normalizes the unitarity triangle, it appears throughout flavor physics. For example, the SM expressions for and for the branching ratios of the golden modes and all contain . Therefore, further improvements—beyond what is achieved here—are warranted, particularly during the course of the Belle II experiment Akeroyd et al. ().

The amplitude for semileptonic decay is expressed in terms of form factors,

(1)
(2)

where and are the (continuum QCD) electroweak currents, , , the velocity transfer , and is the polarization vector of the meson. In the SM, the differential rate for (and the charge-conjugate mode) is given by

(3)

where provides a structure-independent electroweak correction from next-to-leading-order box diagrams, in which a photon or boson is exchanged along with the boson Sirlin (1982). (See Sec. VIII for details.) The rate for (and charge conjugate) is the same as Eq. (3) but with an additional factor on the right-hand side  Ginsberg (1968); Atwood and Marciano (1990), which accounts for the Coulomb attraction of the final-state charged particles.

The notation is conventional, motivated by the heavy-quark limit. In the zero-recoil limit, , one has , and only one form factor survives:

(4)

From Eq. (1), one sees that the needed matrix element is with initial and final states both at rest.

For nonvanishing lepton mass , the rate is multiplied by , and the expressions for and receive corrections proportional to Körner and Schuler (1990). At zero recoil, these corrections reduce to an additional factor on the right-hand side of Eq. (4). Except for , lepton mass effects are not important even at the current level of accuracy.

Because precision is so crucial, the lattice-QCD calculation must be set up in a way that ensures considerable cancellation of all sources of uncertainty. The pioneering work of Hashimoto et al. Hashimoto et al. (1999, 2002) introduced several double ratios to this end. Here, we follow Ref. Bernard et al. (2009) and use a single, direct double ratio

(5)

with all states at rest and the polarization of the aligned with . In the continuum, the denominator of Eq. (5) is unity, by the definition of the flavor quantum numbers. On the lattice, however, it normalizes the flavor numbers and cancels statistical fluctuations. The main uncertainties stem, then, from the chiral extrapolation (the light-quark masses in our data exceed the up and down masses) and discretization and matching errors. In particular, we show how the discretization errors of the analogous ratio of lattice-QCD correlation functions are reduced by use of the ratio.

The rest of this paper is organized as follows. Section II describes the details of the lattice-QCD calculation. We discuss the lattice implementation of Eq. (5), the details of the numerical data, and the general structure of the computed correlation functions. Section III describes our fits to a ratio of correlation functions. Section IV discusses perturbative matching. Section V summarizes the tuning of the bottom- and charm-quark masses and presents results for the and hyperfine splittings. Our extrapolation to the continuum limit and physical light-quark mass is described in Sec. VI. Section VII gives full details of our systematic error analysis. Section VIII provides a discussion of electroweak and electromagnetic effects, which, though separate from our QCD calculation, are needed to obtain . Section IX concludes with final results for and . The appendices contain additional material, including the formulas used for the chiral extrapolation (Appendix A), an estimate of heavy-quark discretization errors (Appendix B), and a thorough discussion of our procedure for tuning the bottom- and charm-quark masses (Appendix C), which also yields the hyperfine splitting.

Ii Lattice setup

In this section we discuss the ingredients of our lattice-QCD calculation. We outline first the generation of ensembles of lattice gauge fields, and then the procedures for computing the three-point correlation functions needed to obtain the double ratio , which is the lattice correlation-function analog of .

ii.1 Simulation parameters

We use the MILC ensembles Bazavov et al. (2010a) of lattice gauge fields listed in Table 1. The ensembles were generated with a Symanzik-improved gauge action Weisz (1983); Curci et al. (1983); Weisz and Wohlert (1984); Lüscher and Weisz (1985a) and 2+1 flavors of sea quarks. The couplings in the gauge action include the one-loop effects of gluons Lüscher and Weisz (1985b) but not of sea quarks Hao et al. (2007); the latter were not yet available when the gauge-field generation began Aubin et al. (2004).

(fm) /  (fm) Volume Sources Configs   
/ 2.4 3.9 6.572   24 628 0.8604
/ 2.4 6.2 6.79   4 2052 0.8688
/ 2.4 4.5 6.76   4 2256 0.8677
/ 2.4 3.8 6.76   4 2108 0.8678
/ 2.9 3.8 6.76   4 2096 0.8678
/ 2.4 5.8 7.11   4 1992 0.8788
/ 2.4 4.1 7.09   4 1928 0.8782
/ 2.7 4.1 7.085   4 984 0.8781
/ 3.4 4.2 7.08   4 1012 0.8779
/ 5.5 4.8 7.075   4 788 0.877805
/ 2.9 6.3 7.48   4 576 0.8881
/ 2.9 4.5 7.47   4 672 0.88788
/ 3.4 4.4 7.465   4 800 0.88776
/ 3.8 4.3 7.46   4 824 0.88764
/ 2.9 4.6 7.81   4 800 0.89511
Table 1: Parameters of the lattice gauge fields. The columns from left to right are the approximate lattice spacing in fm, the sea-quark masses , the linear spatial dimension of the lattice ensemble in fm, the dimensionless factor (with from the Goldstone pion), the gauge coupling, the dimensions of the lattice in lattice units, the number of sources and configurations in each ensemble, and the tadpole improvement factor  (obtained from the average plaquette).

The sea-quark action is the order , tadpole-improved (asqtad) action Blum et al. (1997); Orginos and Toussaint (1999); Lagaë and Sinclair (1999); Lepage (1999); Orginos et al. (1999) for staggered quarks Susskind (1977); Sharatchandra et al. (1981). To reduce the species content from the four that come with staggered fermions, the light quarks (strange quark) are simulated with the square root (fourth root) of the determinant Hamber et al. (1983). At nonzero lattice spacing this procedure introduces small violations of unitarity Prelovšek (2006); Bernard (2006); Bernard et al. (2007); Aubin et al. (2008) and locality Bernard et al. (2006). Considerable numerical and theoretical evidence suggests that these effects go away in the continuum limit, so that the procedure yields QCD Adams (2005); Shamir (2005); Dürr (2005); Bernard (2006); Shamir (2007); Sharpe (2006); Bernard et al. (2008); Kronfeld (2007); Golterman (2008); Donald et al. (2011).

As one can see from Table 1, some ensembles contain independent gauge fields, others –1000. To increase statistics, we re-use each field four times (for  fm, 24 times) by computing quark propagators that are evenly spaced in the time direction with a spatial source origin that is chosen at random from one configuration to the next.

We also use the asqtad action for the light valence (spectator) quark. In this paper, we denote the physical quark masses by , , , and ; the variable spectator mass by ; and the sea-quark masses and , which are fixed within each ensemble. The bare spectator masses are listed in Table 2.

(fm) /
/ 0.0097, 0.0194 1.567 0.0781 0.08354 0.1218 0.08825
/ 0.02 1.525 0.0918 0.09439 0.1259 0.07539
/ 0.01, 0.02 1.531 0.0901 0.09334 0.1254 0.07724
/ 0.007, 0.02 1.530 0.0901 0.09332 0.1254 0.07731
/ 0.005, 0.02 1.530 0.0901 0.09332 0.1254 0.07733
/ 0.0124 1.473 0.0982 0.09681 0.1277 0.06420
/ 0.0062, 0.0124 1.476 0.0979 0.09677 0.1276 0.06482
/ 0.00465 1.477 0.0977 0.09671 0.1275 0.06523
/ 0.0031, 0.0124 1.478 0.0976 0.09669 0.1275 0.06537
/ 0.00155 1.4784 0.0976 0.09669 0.1275 0.06543
/ 0.0072 1.4276 0.1048 0.09636 0.1295 0.05078
/ 0.0036, 0.0072 1.4287 0.1052 0.09631 0.1296 0.05055
/ 0.0025 1.4293 0.1052 0.09633 0.1296 0.05070
/ 0.0018 1.4298 0.1052 0.09635 0.1296 0.05076
/ 0.0028 1.3943 0.1143 0.08864 0.1310 0.03842
Table 2: Valence-quark parameters used in the simulations. The (approximate) lattice spacings  and the sea-quark masses (first two columns) identify the ensemble. Here, denotes the bare masses for the light spectator quarks, and denote the parameters in the SW action, and the rotation parameter in the current. The primes on and distinguishes the simulation from the physical values.

In every case, we compute light-quark propagators with the valence mass equal to the light mass, , and, in several cases, we also compute a partially quenched propagator with .

For the heavy and quarks we use Wilson fermions Wilson (1977) with the Sheikholeslami-Wohlert (SW) action Sheikholeslami and Wohlert (1985), adjusting the parameters in the action according to the Fermilab method El-Khadra et al. (1997). Table 2 also lists the parameters of the heavy-quark action: the hopping parameter (for each quark) and the clover coefficient of the SW action. We use and to denote the values used in the computations, reserving and for those that reproduce the and meson masses most accurately. We set to the value from tree-level tadpole-improved perturbation theory, , with from Table 1. Table 3 gives the values of where the quark mass vanishes for the SW action on each of our ensembles.

(fm) /
/ 2.2215 0.142432
/ 2.8211 0.14073
/ 2.7386 0.14091
/ 2.7386 0.14095
/ 2.7386 0.14096
/ 3.8577 0.139052
/ 3.7887 0.139119
/ 3.7716 0.139134
/ 3.7546 0.139173
/ 3.7376 0.13919
/ 5.3991 0.137582
/ 5.3531 0.137632
/ 5.3302 0.137667
/ 5.3073 0.137678
/ 7.2082 0.13664
Table 3: Derived parameters that enter the simulations. The (approximate) lattice spacings  and the sea-quark masses (first two columns) identify the ensemble. Values for are given in column three, and values for the SW action evaluated on our ensembles are given in column four. For , statistical errors are 0.1 to , and the systematic errors are comparable. For the errors are a few in the last quoted digit.

These values were determined using the methods discussed in Ref. Bernard et al. (2011); note that is only needed in the present work to fix the improvement coefficients that correct the lattice currents described below.

The relative lattice spacing is determined by calculating on each ensemble, where is related to the heavy-quark potential and is defined such that the force between static quarks, Sommer (1994); Bernard et al. (2000). A mass-independent procedure is used to set . This procedure takes the measured values and constructs a smooth interpolation/extrapolation, which we use to replace the measured values with , evaluated now at the physical masses , . Table 3 lists values for each of the ensembles that results from fitting the calculated to the smooth function and extrapolating/interpolating to physical masses. The absolute lattice spacing requires a physical quantity to set the scale. We take the absolute lattice spacing to be fm from the MILC determination of . The value used is explained and justified in Ref. Bazavov et al. (2012).

We have to adjust the light-quark bare masses and the heavy-quark hopping parameters to their physical values a posteriori. The adjustment of the light-quark masses is carried out in the chiral extrapolation, discussed in Sec. VI. For the heavy quarks, we have chosen and in Table 2 close to the physical value based on an initial set of runs that studied a range of but computed only the two-point functions for heavy-strange meson masses. After the full runs, including three-point functions, we re-analyzed the two-point functions to determine more precise values, as discussed in detail in Appendix C. Using information on the dependence, we can then fine-tune our result.

ii.2 correlation functions

To obtain the matrix elements in Eq. (5), we compute the correlation functions

(6)
(7)

and similarly and . Here, and are lattice operators with quantum numbers needed to annihilate and mesons, in the case of with polarization in the  direction; and are lattice currents for transitions. The lattices are gauge-fixed before evaluating the correlation functions so that we can use a smearing function that is extended over a spatial slice.

We form the interpolating operators from a staggered fermion field  and heavy-quark field  in the SW action:

(8)
(9)
(10)

where is a spatial smearing function. . The free Dirac index on can be interpreted as a taste index (in which case we average over taste) Kronfeld (2007), or one can promote to a four-component field Wingate et al. (2003), which leads to the same results for the correlation functions of bilinear operators.

We employ two smearing functions. One is the local . The other is the ground-state 1S wavefunction of the Richardson potential. See Ref. Bazavov et al. (2012) for details.

We define the lattice vector and axial-vector currents to be

(11)
(12)

where , are flavor indices. The fermion field includes a correction factor to reduce discretization effects El-Khadra et al. (1997),

(13)

where is a nearest-neighbor covariant difference operator. Its coefficient is set to its value in tree-level tadpole-improved perturbation theory, where it does not depend on the other quark in the current. The matrix elements of the lattice currents satisfy ( means “has the same matrix elements as”) Kronfeld (2000); Harada et al. (2002)

(14)

where is the continuum current corresponding to the lattice current and the matching factors are defined such that Eq. (14) holds. In practice, can be determined only approximately, via either perturbative or nonperturbative methods. Thus, for tree-level matching of or , for one-loop matching, etc. Nonperturbative matching schemes could be set up, which would remove all powers of . Here we implicitly use nonperturbative matching for flavor-diagonal , one-loop matching for suitable ratios of factors (see below), and tree-level matching for . Higher-loop and nonperturbative calculations, except for , are not available.

In the double ratio like Eq. (5) but with matrix elements of lattice currents, the following ratio of matching factors remains:

(15)

In this ratio, all corrections associated with wave-function renormalization cancel out, leaving only vertex diagrams. Each contains the difference between continuum and lattice vertex diagrams, and the ratio introduces further cancellations. It is not surprising, then, that one-loop calculations of yield very small coefficients of  Harada et al. (2002).

With the Fermilab method applied to the SW action, the Lagrangian also leads to discretization effects of order , where counts, as above, the matching of the SW (clover) term. Again, one-loop matching is not completely available (see Ref. Nobes and Trottier (2006)), so we use tree-level matching. Table 2 lists the values of used in this work. Appendix B discusses the discretization effects in (as extracted here) in detail.

For large enough time separations and , the correlation function

(16)

where and are the masses of the and mesons and . The omitted terms from higher-mass states are discussed in Sec. III. The other correlation functions , , have analogous large-time behavior. Therefore, the ratio of correlation functions

(17)

where

(18)

is a lattice version of , up to the matching factor  and discretization errors. The analysis of to extract is discussed in Sec. III, the calculation of is discussed in Sec. IV, the light-quark discretization errors are analyzed in Sec. VI, and the heavy-quark discretization errors are derived in Appendix B.

Above we mentioned that we increase statistics by choosing four (24 at  fm) sources. This means we choose four (24) origins in Eqs. (6) and (7). We do so by picking at random four (24) equally separated timeslices for . On each timeslice, we choose a completely random point for .

Starting at each origin [ in Eqs. (6) and (7)], we construct the three-point correlation functions as follows. We compute the parent heavy-quark propagator from smeared to all points, in particular . We also compute the spectator staggered-quark propagator from to all points. At time , we convolve this propagator with the Dirac matrix and smearing function of the sink, projecting onto a fixed momentum (here, ). This combination is used for a further inversion for the daughter heavy quark; this inversion yields a sequential propagator encoding the propagation of the spectator quark, a flavor change at the sink, and (reverse) propagation of the daughter quark back to the decay. This sequential propagator and the parent propagator are then inserted into the appropriate trace over color and Dirac indices.

Iii Analysis of correlation functions

To obtain from with sufficient accuracy, we have to treat the excited states [denoted by in Eq. (16)] carefully. From the transfer-matrix formalism, one finds

(19)

where even and label excitations of desired parity, and odd and label excitations of opposite parity. The appearance of the opposite-parity states and their oscillating time dependence are consequences of using staggered fermions for the spectator quark. The are transition matrix elements, multiplied by uninteresting factors. For the desired , these factors cancel in .

In practice, we can choose the time separations such that only the lowest-lying states of each parity make a significant contribution. As discussed in detail in Ref. Bernard et al. (2009), it is advantageous to smear over time in a way that suppresses the opposite-parity state, and define

(20)

which is very close to , with small time-dependent effects that one can disentangle via a fit to the dependence.

The average in Eq. (20) is designed to suppress the contribution from oscillating states that changes sign only when the total source-sink separation is varied (the “same-sign” oscillating-state contributions). The double ratio, including the leading effects of the wrong-parity states, is

(21)

where the function contains the oscillating-state contributions, and is given by

(22)

The terms in square brackets in Eq. (22) are the suppression factors for the oscillating state contributions. The are the splittings between the ground-state masses and the opposite-parity masses, and their values can be computed precisely from fits to two-point correlators. We find values for these splittings in the range between about 0.1 and 0.4 in lattice units. With these values of the parameters the “same-sign” contributions [the third term in Eq. (22)] are suppressed by a factor of –20 by Eq. (20), where the suppression is greater at finer lattice spacings. The other oscillating-state contributions change sign as a function of and are given by the first two terms in Eq. (22). These contributions are very small for our double ratio, and they are further suppressed by a factor of 2 by the average in Eq. (20).

Figure 2: at on 0.12 fm (left) and on 0.09 fm (right) lattice spacings.
Figure 3: at on 0.12 fm for two different combinations of source-sink separations.

Figure 2 shows and for two different, representative ensembles. One can see that the plateau is lower for odd total source-sink separation than for even total source-sink separation, whether the odd source-sink separation is larger or smaller than the even source-sink separation. This feature holds for all ensembles. It suggests that the “same-sign” oscillating states are visible in our data, and are comparable to, but somewhat larger than, our current statistical errors. The average of Eq. (20) suppresses this effect to around on our coarser lattices and around on our finest lattices. This effect is negligible compared to other errors. The fact that this effect is visible independently of whether the odd source-sink separation is larger or smaller than the even source-sink separation indicates that this effect is larger than other excited-state contributions and that within the current statistical precision of our data, these can also be neglected. That this is the case is verified by a calculation at a larger source-sink separation on the 0.12 fm ensemble. Figure 3 shows a comparison between the square-root of the average Eq. (20) for two different combinations of source-sink separations. The larger source-sink separation is computed with 16 time sources on 2256 configurations, compared with four time sources on the same configurations for the smaller separation. The source was moved around the lattice randomly with a different seed for the two calculations, so we expect the ratios to be less correlated than is typical for quantities computed on the same configurations. The agreement between the best fits to the different source-sink separations is good to the 1 level, as expected if residual excited state contamination is small. Since ordinary excited state contamination would tend to cause the plateau fit to be too high, as can be seen by the higher values of near the source and sink, this contamination must be negligible within our current statistics because the fit with the larger separation and smaller contamination gives a slightly higher plateau value.

Figure 4: at on 0.12 fm (left), 0.09 fm (right), and 0.06 fm (bottom). The plateau fits are shown with 1 error bands.

The square-root of the average Eq. (20) is shown in Fig. 4 for 0.12 fm, 0.09 fm, and 0.06 fm lattice spacings. These plots show data at unitary (full QCD) points, with valence spectator- and light sea-quark masses equal to . The square-root of is fit to a constant in the identified plateau region, including the full covariance matrix to determine the correlated and to ensure that the fits yielded acceptable  values. The fits are shown in Fig. 4 superimposed over the data with error bands. Source-sink separations and plateau ranges are approximately the same in physical units for all lattice spacings. Time ranges for fits, their  values, and the raw values for are given in Table 4.

(fm) /    fit range   value
/ 10, 11  5-7 0.85 0.9141(51)
/ 12, 13  5-8 0.80 0.9035(28)
/ 12, 13  5-8 0.97 0.9052(44)
/ 12, 13  5-8 0.63 0.9160(53)
/ 12, 13  5-8 0.68 0.9143(55)
/ 17, 18  7-11 0.63 0.9162(31)
/ 17, 18  7-11 0.54 0.9135(45)
/ 17, 18  7-11 0.78 0.9212(73)
/ 17, 18  7-11 0.95 0.9092(68)
/ 17, 18  7-11 0.79 0.9208(90)
/ 24, 25  8-14 0.84 0.9126(50)
/ 24, 25  8-14 0.93 0.9097(64)
/ 24, 25  8-14 0.13 0.9073(67)
/ 24, 25  8-14 0.55 0.9147(64)
/ 32, 33  7-14 0.87 0.9029(45)
Table 4: Fit results for double ratios at the full QCD points. The (approximate) lattice spacings  and the sea-quark masses (first two columns) identify the ensemble. The third column is the pair of spectator quark source-sink separations, the fourth is the time-slice fit range, the fifth is the  value of the fit, and the sixth is the value of determined from the fit.

Iv Perturbation theory for

As discussed in Sec. II.2, we need the ratio of matching factors, , defined in Eq. (15). This ratio has been calculated in one-loop perturbation theory, which will be discussed in detail in another publication. The perturbative expansion for is

(23)

where we make explicit a choice of scheme and scale for the perturbative series. The calculation of is a straightforward extension of the work in Ref. Harada et al. (2002), modified to use the improved gluon propagator.

For the expansion parameter , we would like to make a choice that prevents large logarithms associated with the function from making the neglected terms unnecessarily large. Brodsky, Lepage, and Mackenzie Brodsky et al. (1983) discussed how to do so by exploiting the dependence of the second order in , and Lepage and Mackenzie Lepage and Mackenzie (1993) explained how to define an equivalent scale choice when the second order is not yet available. The Lepage-Mackenzie version requires a coefficient defined by weighting the Feynman integral for with an additional factor of , where is the gluon momentum in the one-loop diagram(s). Then the recommended (and empirically successful Lepage and Mackenzie (1993); Harada et al. (2003)) scale  is given through

(24)

when the scheme is the  scheme, such that the interquark potential in momentum space is .

Unfortunately, as the heavy-quark masses vary over the range of interest, nearby zeroes of the numerator and denominator in Eq. (24) lead to physically unreasonable values for . Fortunately, the way to deal with such cases has been spelled out by Hornbostel, Lepage, and Morningstar (HLM) Hornbostel et al. (2003). The HLM method requires integrals weighted by higher powers of . This prescription results in values for that are close to . We therefore use to obtain the listed in Table 5.

(fm) /       
0.0097/ 0.0484 0.3589 0.99422(4)
0.02/ 0.05 0.3047 0.99650(5)
0.01/ 0.05 0.3108 0.99623(5)
0.007/ 0.05 0.3102 0.99618(5)
0.005/ 0.05 0.3102 0.99617(5)
0.0124/ 0.031 0.2582 0.99978(4)
0.0062/ 0.031 0.2607 0.99963(4)
0.00465/ 0.031 0.2611 0.99957(4)
0.0031/ 0.031 0.2619 0.99950(4)
0.00155/ 0.031 0.2623 0.99946(4)
0.0072/ 0.018 0.2238 1.00334(3)
0.0036/ 0.018 0.2245 1.00323(3)
0.0025/ 0.018 0.2249 1.00317(3)
0.0018/ 0.018 0.2253 1.00312(3)
0.0028/ 0.014 0.2013 1.00608(2)
Table 5: One-loop estimate of . The first two columns label each ensemble with the approximate lattice spacing in fm and the sea simulation light- and strange-quark masses. The third column is with . The fourth column is on that ensemble with statistical errors from the VEGAS evaluation of the one-loop coefficients.

As expected, varies somewhat as a function of lattice spacing. It is even slightly different from ensemble to ensemble at the same nominal lattice spacing, because these ensembles have slightly different lattice spacings.

V Heavy-quark mass tuning

Our approach to tuning is similar to that described in Ref. Bernard et al. (2011), and a detailed description of the current approach is given in Appendix C. We start with the lattice dispersion relation

(25)

where defines the meson rest mass and the kinetic mass is given by

(26)

The meson masses differ from the corresponding quark masses, and , by binding-energy effects. In the Fermilab method, the lattice pole energy is fit to the dispersion relation Eq. (25), and is adjusted so that the kinetic mass agrees with experiment. We tune to the experimental and meson masses to obtain and , respectively.

The simulation values differ from our current best estimates of these parameters because of improvements in statistics and methodology since the initial tuning runs. Table 6 shows our best estimates of , along with errors. The first error is a combination of statistical and fitting systematics, and the second error is that due to fixing the lattice scale. For comparison, Table 6 also shows the values used in the runs.

(fm) /
/ 0.0775(16)(3) 0.12237(26)(20) 0.0781 0.1218
/ 0.0879(9)(3) 0.12452(15)(16) 0.0918 0.1259
/ 0.0868(9)(3) 0.12423(15)(16) 0.0901 0.1254
/ 0.0868(9)(3) 0.12423(15)(16) 0.0901 0.1254
/ 0.0868(9)(3) 0.12423(15)(16) 0.0901 0.1254
/ 0.0972(7)(3) 0.12737(9)(14) 0.0982 0.1277
/ 0.0967(7)(3) 0.12722(9)(14) 0.0979 0.1276
/ 0.0966(7)(3) 0.12718(9)(14) 0.0977 0.1275
/ 0.0965(7)(3) 0.12714(9)(14) 0.0976 0.1275
/ 0.0964(7)(3) 0.12710(9)(14) 0.0976 0.1275
/ 0.1054(5)(2) 0.12964(4)(11) 0.1048 0.1295
/ 0.1052(5)(2) 0.12960(4)(11) 0.1052 0.1296
/ 0.1051(5)(2) 0.12957(4)(11) 0.1052 0.1296
/ 0.1050(5)(2) 0.12955(4)(11) 0.1052 0.1296
/ 0.1116(3)(2) 0.130921(16)(7) 0.1143 0.1310
Table 6: Errors in the tuned parameters. The (approximate) lattice spacings  and the sea-quark masses (first two columns) identify the ensemble. The third and fourth columns are the tuned values for the and quarks, respectively. The first error is the statistics plus fitting error, and the second is an error due to the uncertainty in the lattice scale. The fifth and six columns are the values used in the simulations.

A detailed discussion of how the tuned values of are obtained is given in Appendix C. As a cross-check of our tuning procedure, we calculate the hyperfine splittings and . In Appendix C.4 we find

(27)

where the error includes statistics and the sum of all systematic errors in quadrature. These are in good agreement with the experimental values  MeV and  MeV.

We correct our values of for the mistuning of using information on the heavy-quark mass dependence from an additional run with nearer their physical values on the coarse ensemble with . To apply the correction we exploit information from heavy-quark effective theory (HQET); the form factor at zero-recoil has the heavy-quark expansion Falk and Neubert (1993); Mannel (1994)

(28)

up to order , where is a factor that matches HQET to QCD and the ’s are long-distance matrix elements of the HQET. Heavy-quark symmetry forbids terms of order at zero-recoil Luke (1990). The form factor depends on both the bottom quark mass and the charm quark mass; we correct for this dependence and propagate the uncertainty due to the error in to the form factor before performing the chiral/continuum extrapolation. The leading dependence is given by the term that is inversely proportional to in brackets in Eq. (28), and this dependence, inversely proportional to for fixed charm-quark mass, is the one used to correct the form factor for the mistuning in . The leading charm-quark mass dependence is, however, given by the term that is inversely proportional to the charm quark mass squared. Thus, we determine the adjustment that must be made from the simulated form factor to the tuned value using

(29)

where is the kinetic or quark mass, and sets the relative lattice spacing on different ensembles. The slope parameters are determined by a linear interpolation between the two sets of points shown in Fig. 5.

Figure 5: at different values close to the tuned and quark masses. Each is plotted as a function of the leading (assuming is sufficiently heavier than ) heavy-quark mass dependence in Eq. (28), and for and , respectively.

One of these points in each of these plots is from our original production run, while the other points are from runs where were separately varied and chosen to be closer to their tuned values.

The slopes are also used to propagate the errors in the tuned kappa values due to “statistics and fitting” to the errors in each individual data point before performing the chiral/continuum extrapolation. This is done by inflating the jackknife error of on each data point by adding to it in quadrature the parametric error in due to the “statistics and fitting” part of the tuning error. We make the assumption that the statistics and fitting errors in the tuned values on different ensembles are independent of one another, though we also test the size of the additional error induced if this assumption is not true and find that it is small. The tuning “statistics and fitting” error is thus directly incorporated into the statistical error of . The scale error in the tuned values, however, is correlated across ensembles, and is therefore treated as a separate systematic error.

Vi Chiral-continuum extrapolation

Because the light and -quark masses used in the calculation are heavier than the physical ones, an extrapolation in quark mass is necessary. This extrapolation can be controlled using an appropriate chiral effective theory, where one can also incorporate discretization effects particular to staggered quarks. The chiral effective theory that incorporates these effects is rooted staggered chiral perturbation theory (rSPT), which was extended to include heavy-light quantities in Ref. Aubin and Bernard (2006).

There are discretization effects that are particular to staggered quark actions. The staggered quark discretization only partially solves the fermion doubling problem, reducing the number of species from 16 to 4. There remain unphysical species of quarks, commonly referred to as tastes. Quarks of different tastes can exchange high momentum gluons with momenta of order the lattice cutoff, and this exchange breaks the degeneracy in the pion spectrum for pions made of quarks of different tastes. This taste-symmetry breaking leads to the staggered theory having 16 light pseudoscalar mesons instead of 1.

The tree-level relation in the chiral theory between the pseudoscalar meson masses and the quark masses is given by