# Lattice QCD form factor for at zero recoil with non-perturbative current renormalisation

###### Abstract

We present details of a lattice QCD calculation of the axial form factor at zero recoil using the Highly Improved Staggered Quark (HISQ) formalism on the second generation MILC gluon ensembles that include up, down, strange and charm quarks in the sea. Using the HISQ action for all valence quarks means that the lattice axial vector current that couples to the can be renormalized fully non-perturbatively, giving a result free of the perturbative matching errors that previous lattice QCD calculations have had. We calculate correlation functions at three values of the lattice spacing, and multiple ‘’-quark masses, for physical and . The functional dependence on the -quark mass can be determined and compared to Heavy Quark Effective Theory expectations, and a result for the form factor obtained at the physical value of the -quark mass. We find . This is in agreement with earlier lattice QCD results, which use NRQCD quarks, with a total uncertainty reduced by more than a factor of two. We discuss implications of this result for the axial form factor at zero recoil and for determinations of .

HPQCD collaboration

## I Introduction

The study of quark flavour-changing interactions is a key component of the search for physics beyond the Standard Model (SM). There are currently a number of related tensions between experimental measurements and SM predictions Wei et al. (2009); Lees et al. (2012a, b, 2013); Aaij et al. (2014a, b); Huschle et al. (2015); Aaij et al. (2016a, 2015a, 2015b, b); Wehle et al. (2017); Sato et al. (2016); Hirose et al. (2018); Aaij et al. (2018a, b, 2017); Sirunyan et al. (2018); Aaboud et al. (2018), along with discrepancies between systematically independent determinations of Cabibbo-Kobayashi-Maskawa (CKM) matrix elements Amhis et al. (2017); Bevan et al. (2014); Alberti et al. (2015). A more precise understanding of these processes is needed to resolve these issues.

The decay (and its charge conjugate, that we simply abbreviate to from now on) supplies one of the three methods used for precisely determining the CKM element Schroder (1994); Bortoletto et al. (1990); Fulton et al. (1991); Albrecht et al. (1992); Barish et al. (1995); Buskulic et al. (1996, 1994); Abbiendi et al. (2000); Abreu et al. (2001); Adam et al. (2003); Abdallah et al. (2004); Aubert et al. (2008a, b, 2009); Dungel et al. (2010); Abdesselam et al. (2017); Bailey et al. (2014); Abdesselam et al. (2018). Measurements of branching fractions are extrapolated through to the zero recoil point to deduce , where is the value of the only form factor contributing at zero recoil. Then a determination of in the Standard Model (via Lattice QCD Bailey et al. (2014); Harrison et al. (2018)) can be divided out to infer .

is an important quantity and needs to be determined accurately. It constrains one side of the unitarity triangle via the ratio . It is also a dominant uncertainty in the determination of the -violation parameter (where there is currently tension between the SM and experiment, see for example Bailey et al. (2018a)).

Previous determinations of have shown systematic discrepancies with each other. The two competing values were those derived from exclusive decays ( and with giving the more accurare result), and inclusive (, where is any charmed hadronic state). In 2016 the Heavy Flavour Averaging Group (HFLAV) gave a value derived from exclusive decays of and from inclusive decays, using the kinetic scheme, of Amhis et al. (2017). It has since been suggested, based on unfolded Belle data Abdesselam et al. (2017), that the tension seen here arose (at least partly) from the use of a very constrained parameterization in the extrapolation of the experimental decay rates to zero recoil Bernlochner et al. (2017); Bigi et al. (2017); Grinstein and Kobach (2017). Recent exclusive determinations of have then used a less constrained parameterisation to give a larger, and less precise, result for that is no longer in tension with the inclusive result. For example, the Particle Data Group quote Tanabashi et al. (2018). Although this means that inclusive and exclusive determinations now agree, it clearly points to the need for more work to improve the accuracy of the exclusive result. On the theory side a better understanding of the form factors for from lattice QCD is required, both at zero recoil and away from zero recoil.

Another motivation for studying is the tension between SM and experimental determinations of the ratio ( or ). The combined statistical significance of the anomalies in and currently stands at Amhis et al. (2017). More precise measurements and predictions will either confirm or dismiss a new physics explanation.

The weak decay process is very similar to and could also be used to determine and test the SM. It is feasible to study this decay at LHC and from the theoretical side it is a more attractive channel than . The absence of valence light quarks means that lattice QCD results have smaller statistical errors and are less computationally expensive. Finite-volume effects and the dependence on quark masses (for quarks in the sea) are also smaller. The has no Zweig-allowed strong decay mode, unlike the , and is in fact a relatively long-lived particle Donald et al. (2014a) that can be considered ‘gold-plated’ in lattice QCD. This makes the both a useful test bed for lattice techniques (that may be later used to study decays) and a key decay process for which to make predictions ahead of experimental results.

Lattice QCD calculations have shown that several weak decay form factors are relatively insensitive to whether the spectator quark is a or quark Bailey et al. (2012); Koponen et al. (2013); Monahan et al. (2017). A combination of chiral perturbation theory and Heavy Quark Symmetry Jenkins and Savage (1992) backs up this expectation for decays. We can therefore expect the form factors to be very similar for and . A recent lattice calculation Harrison et al. (2018) found an insignificant difference at zero recoil: . Information from the study of can then be applicable to .

Lattice QCD calculations of the form factors at zero recoil have so far been performed by two collaborations using different methods. The Fermilab Lattice and MILC collaborations calculated in Bernard et al. (2009); Bailey et al. (2014) using the Fermilab action for both and quarks El-Khadra et al. (1997) and asqtad quarks Bazavov et al. (2010a). More recently the HPQCD collaboration computed both and Harrison et al. (2018) using improved NRQCD quarks Lepage et al. (1992); Dowdall et al. (2012) and Highly Improved Staggered (HISQ) and quarks Follana et al. (2007). The RBC/UKQCD Flynn et al. (2016) and LANL-SWME Bailey et al. (2018b) collaborations are also working towards these form factors using variants of the Fermilab action for heavy quarks and JLQCD has a calculation in progress using Möbius domain-wall quarks Kaneko et al. (2018).

The formalism to use for the heavy quarks is a major consideration in designing a lattice QCD calculation to determine these form factors. Most of the calculations discussed in the previous paragraph (apart from the JLQCD calculation) use approaches that make use of the nonrelativistic nature of heavy quark bound states to tune the (and in some cases also ) quark masses. This avoids potentially large discretisation effects appearing in the results in the form of a systematic error of size , where is an integer that depends on the level of improvement in the action. The absence of these discretisation errors means that quarks can be handled on relatively coarse lattices where . However the price to be paid is that the current operator that couples to the boson is also implemented within a nonrelativistic framework and must then be renormalised to match the appropriate operator in continuum QCD. This matching can be done using perturbative lattice QCD but has only be done through for these actions Harada et al. (2002); Monahan et al. (2013). This leaves a substantial source of uncertainty from missing higher-order terms in the perturbative matching that is not easily reduced. This matching uncertainty contributes 80% of the final error in the HPQCD calculation Harrison et al. (2018) and 30% in the Fermilab/MILC calculation Bailey et al. (2014) because of the differing allowances for missing higher-order terms.

Here we report details and results of a calculation of the form factor at zero recoil using an approach free of perturbative matching uncertainties. We perform our calculation on the second-generation MILC ensembles Bazavov et al. (2010b, 2013), including effects from 2+1+1 flavours in the sea using the HISQ action. We also use the HISQ action for all valence quarks. We obtain results at a number of differing masses for the (we refer to this generically as the heavy quark ), and perform an extrapolation to . By using only HISQ quarks, we can obtain the normalizations of all required currents fully non-perturbatively. We refer to this as the heavy-HISQ approach. By using many heavy masses and multiple values of the lattice spacing, including very fine lattices, we can model both the form factor dependence on the heavy mass, and the discretisation effects associated with using large values.

The heavy-HISQ approach was developed by HPQCD to compute meson masses and decay constants McNeile et al. (2012a, b) and the quark mass McNeile et al. (2010); Chakraborty et al. (2015). It is also now being used by other collaborations for these calculations Bazavov et al. (2018a); Petreczky and Weber (2019). A proof-of-principle application of heavy-HISQ to form factors was given in Lytle et al. (2016); Colquhoun et al. (2016) for decays, showing that the full range of the decay could be covered. Here we extend the approach to form factors for decays but working only at zero recoil, a straightforward extension of earlier work. Using the heavy-HISQ approach also has the added benefit of eludicating the dependence of form factors on heavy quark masses, meaning that we can test expectations from Heavy Quark Effective Theory (HQET).

This article is structured in the following way: Section II defines the form factor and gives details of the lattice calculation, including the nonperturbative normalisation and extrapolation in heavy-quark mass; Section III presents our results and compares to earlier calculations and Section IV gives our conclusions and outlook. In the appendix, we give details of a number of tests we performed on the correlator fits and the continuum, chiral and heavy-quark extrapolations.

## Ii Calculation Details

### ii.1 Form Factors

The differential decay rate for the decay is given in the SM by

(1) | ||||

where , is the 4-velocity of each meson, and is a known function of with (see, for example, appendix G of Harrison et al. (2018)). accounts for electroweak corrections from diagrams where photons or s are exchanged in addition to a , as well as the Coulomb attraction of the final-state charged particles Sirlin (1982); Ginsberg (1968); Atwood and Marciano (1990). The differential decay rate for the is identical.

The form factor is a linear combination of hadronic form factors that parameterize the vector and axial-vector matrix elements between initial and final state hadrons. A common choice of parameterization used in the context of Heavy Quark Effective Theory (HQET) is Richman and Burchat (1995)

(2) | ||||

(3) | ||||

where is the vector current and is the axial-vector current. is the polarization 4-vector of the final state.

At zero recoil (), the vector matrix element vanishes, the axial-vector element simplifies to

(4) |

and reduces to

(5) |

Our goal is to compute .

All we need to do this is the matrix element with both the and at rest, with the polarization in the same direction as the (spatial) axial-vector current.

### ii.2 Lattice Calculation

The gluon field configurations that we use were generated by the MILC collaboration Bazavov et al. (2010b, 2013). Table 1 gives the relevant parameters for the specific ensembles that we use. The gluon field is generated using a Symanzik-improved gluon action with coefficients calculated through Hart et al. (2009). The configurations include the effect of 2+1+1 flavours of dynamical quarks in the sea (,,,, with ), using the HISQ action Follana et al. (2007). In three of the four ensembles (fine, superfine and ultrafine), the bare light quark mass is set to . The fact that the value is unphysically high is expected to have only a small effect on , because there are no valence light quarks. The effect is quantified here by including a fourth ensemble (fine-physical) with (approximately) physical .

We use a number of different masses for the valence heavy quark, . This is in order to resolve the dependence of on the heavy mass, so that an extrapolation to can be performed. By varying the heavy mass on each ensemble and by using ensembles at varying small lattice spacing, we can resolve both the discretisation effects that grow with heavy quark mass () and the physical dependence of the continuum form factor on .

set | handle | T | |||||||||
---|---|---|---|---|---|---|---|---|---|---|---|

1 | fine | 1.9006(20) | 0.0074 | 0.037 | 0.440 | 0.0376 | 0.45 | 0.5, 0.65, 0.8 | 14,17,20 | ||

2 | fine-physical | 1.9518(7) | 0.0012 | 0.0363 | 0.432 | 0.036 | 0.433 | 0.5, 0.8 | 14,17,20 | ||

3 | superfine | 2.896(6) | 0.0048 | 0.024 | 0.286 | 0.0234 | 0.274 | 0.427, 0.525, 0.65, 0.8 | 22,25,28 | ||

4 | ultrafine | 3.892(12) | 0.00316 | 0.0158 | 0.188 | 0.0165 | 0.194 | 0.5, 0.65, 0.8 | 31,36,41 |

Staggered quarks have no spin degrees of freedom, so that solution of the Dirac equation on each gluon field is numerically fast. The remnant of the doubling problem means that quark bilinears of specific spin-parity have multiple copies, called ‘tastes’ Follana et al. (2007). They differ in the amount of point-splitting between the fields and the space-time dependent phase needed to substitute for the appropriate matrix. In this calculation we can use only local (non point-split) bilinears, which is an advantage in terms of statistical noise, since no gluon fields are included in the current operator. In the standard staggered spin-taste notation, the operators that we use are: pseudoscalar, ; vector, and axial-vector, .

We compute several ‘two-point’ correlation functions on the ensembles detailed in table 1, combining HISQ propagators from solving the Dirac equation for each random wall time source. These correlation functions take the form

(6) | ||||

where represents a functional integral, are valence quark fields of the flavours the meson is charged under, is the spin-taste structure of and is the staggered quark normalisation for closed loops. The random-wall source and the sum over at the sink project onto zero spatial momentum. We compute the correlation functions for all values, i.e. .

The correlation function for the heavy-strange pseudoscalar meson, , with valence quark content and spin-taste structure is constructed from HISQ propagators as:

(7) |

Here is a HISQ propagator for flavour , the trace is over color and is the staggered quark normalisation. and and the sum is over spatial sites , . We also compute correlators for a charm-strange vector meson , with structure , using

(8) |

We average over polarisations, .

We also compute correlation functions for two tastes of pseudoscalar heavy-charm mesons denoted and respectively. has spin-taste structure , and has structure . correlators are computed using Eq. (7) (with replaced with ), while correlators are given by

(9) |

where we use the notation . These correlators will be used to normalise the axial vector current as discussed in Section II.4.

A useful physical proxy (that does not run) for the quark mass is that of the pseudoscalar meson made from that flavour of quark. It is therefore also useful, for our heavy quark mass extrapolation, to calculate correlation functions for heavy-heavy pseudoscalars, denoted , with spin-taste structure using Eq. (7). Likewise, to test the impact of any mistuning of the charm and strange quark masses, we also determine and correlators analogously. We can tune the and masses using the experimental values for the and masses, allowing for slight shifts from missing QED effects and the fact that we do not allow these mesons to annihilate to gluons McNeile et al. (2010). The mass of the meson (which is not a physical state) can be fixed in lattice QCD from the and meson masses Davies et al. (2010a); Dowdall et al. (2013).

We then generate the ‘three-point’ correlation functions that contain the to transition.

(10) | ||||

Our source is given spin-taste , the sink, , and the current insertion . This gives the required cancellation of tastes within the three-point function Donald et al. (2014b). In terms of HISQ propagators

(11) |

where we fix , and . We compute the three-point correlation functions for all values within , and 3 values that vary between ensembles and are given in Table 1. We average over the 3 directions for for increased statistical precision.

### ii.3 Analysis of Correlation Functions

We use simultaneous Bayesian fits Lepage et al. (2002); cor (2018) to extract the axial vector matrix element and meson masses from the two- and three-point correlation functions. This allows us to include the covariance between results at different heavy quark masses on a given ensemble into our subsequent fits in Section II.5.

We fit the two-point correlation functions using the functional form

(12) | ||||

where is the temporal extent of the lattice, and , are fit parameters, with the excited-state energy parameters implemented as energy differences to the state below Lepage et al. (2002). The second term accounts for the presence of opposite-parity states that contribute an oscillating term to the correlation function when using staggered quarks Follana et al. (2007). These terms do not appear when is a pseudoscalar with a quark and antiquark of the same mass, so in the and cases the second term is not required. For all correlator fits we set ; this allows the impact of systematic effects from excited states to be included in the ground-state parameters that we are interested in.

The three-point correlation functions have the fit form

(13) | ||||

This includes fit parameters common to the fits of the and two-point correlators, along with new fit parameters .

We perform a single simultaneous fit containing each correlator computed (, and three-point) for each ensemble. We set gaussian priors for the parameters , and log-normal priors for all other parameters. Using log-normal distributions forbids energy differences and amplitudes (which can be taken to be positive here) from moving too close to zero or changing sign, improving stability of the fit.

Ground state energies are given priors of , where and are the masses of the appropriate quarks, and is the confinement scale, which we set to 0.5GeV. For or , this corresponds to the leading order HQET expression for a heavy meson mass. Ground-state energies of oscillating states, , are given priors of . Excited state energy differences, , are given prior values . Priors for ground state amplitudes , are set from plots of effective amplitudes. The resulting priors always have a variance at least 10 times that of the final result for the ground-state. We use log(amplitude) priors of -1.20(67) for non-oscillating excited states and -3.0(2.0) for oscillating excited states. The ground-state non-oscillating to non-oscillating 3-point parameter, is given a prior of , and the rest of the 3-point parameters are given .

is the mass of the ground-state meson in lattice units. The masses and can both be used as proxys for in the extrapolation to . The annihilation amplitude for an -meson is given (in lattice units) by

(14) |

The (as yet unnormalised) matrix element that we need to obtain is given by

(15) |

To ensure that truncating the sum over states at accounts for the full systematic error from excited states, we cut out some data very close to the sources and sinks, where even higher excited states might have some effect. To do this we only include data with and in the two-point case and in the three-point case. We can in principle use a different for every correlation function included in our fit, but we do not use a big range of values. They range from 1 to 3 for the three-point functions and up to 8 for the two-point functions.

The determination and minimisation of the function in our fit procedure requires the inversion of the covariance matrix that captures correlations between the different pieces of ‘data’ (correlation functions) in our fit. The low eigenmodes of the correlation matrix are not well determined with the statistics that we have and so we implement an SVD (singular value decomposition) cut in the inversion of the correlation matrix to avoid underestimating the uncertainty in the parameters of the fit cor (2018). This replaces correlation matrix eigenvalues below , equal to svdcut times the largest eigenvalue, with . is estimated using the diagnosis tools in the Corrfitter package cor (2018) and corresponds typically to an svdcut of here.

Figure 1 summarises stability tests of our fits, focussing on the key parameter that is converted to the ground-state to ground-state transition amplitude using Eq. (15).

The fit parameters determined by our fits that we will use to calculate the physical value for are given in Table 2. Notice that the statistical errors on the results grow with the heavy quark mass. This is a well understood problem in lattice heavy-light meson physics (see, for example Davies et al. (2010b)). Our method here has the advantage of including information from lighter-than- heavy quarks with improved statistical precision.

Set | |||||||||
---|---|---|---|---|---|---|---|---|---|

1 | 0.5 | 0.9255(20) | 0.95972(12) | 0.96616(44) | 1.419515(41) | 0.186299(70) | 1.471675(38) | 1.367014(40) | 0.313886(75) |

0.65 | 0.9321(22) | 1.12511(16) | 1.573302(40) | 0.197220(77) | 1.775155(34) | ||||

0.8 | 0.9434(24) | 1.28128(21) | 1.721226(39) | 0.207068(78) | 2.064153(30) | ||||

2 | 0.5 | 0.9231(21) | 0.95462(12) | 0.93976(42) | 1.400034(28) | 0.183472(62) | 1.470095(25) | 1.329291(27) | 0.304826(52) |

0.8 | 0.9402(27) | 1.27577(22) | 1.702456(23) | 0.203407(45) | 2.062957(19) | ||||

3 | 0.427 | 0.9107(46) | 0.77453(24) | 0.63589(49) | 1.067224(46) | 0.126564(70) | 1.233585(41) | 0.896806(48) | 0.207073(96) |

0.525 | 0.9165(49) | 0.88487(31) | 1.172556(46) | 0.130182(72) | 1.439515(37) | ||||

0.65 | 0.9246(65) | 1.02008(39) | 1.303144(46) | 0.133684(75) | 1.693895(33) | ||||

0.8 | 0.9394(66) | 1.17487(54) | 1.454205(46) | 0.137277(79) | 1.987540(30) | ||||

4 | 0.5 | 0.9143(51) | 0.80245(24) | 0.47164(39) | 1.011660(32) | 0.098970(52) | 1.342639(65) | 0.666586(89) | 0.15412(17) |

0.65 | 0.9273(62) | 0.96386(33) | 1.169761(34) | 0.100531(60) | 1.650180(56) | ||||

0.8 | 0.9422(72) | 1.11787(43) | 1.321647(37) | 0.101714(70) | 1.945698(48) |

### ii.4 Normalisation of the Axial Current

The partially-conserved axial-vector current for the HISQ action is a complicated linear combination of one-link and three-link lattice currents. In this study we use only local axial vector currents. This simplifies the lattice QCD calculation but creates the need for our resulting current matrix element to be multiplied by a matching factor to produce the appropriate continuum current. We determine via a fully non-perturbative method Donald et al. (2014b).

We use the fact that the staggered local pseudoscalar current of spin-taste , multiplied by the sum of its valence quark masses, is absolutely normalized via the PCAC relation. From the two-point and correlator fits we can extract the decay amplitudes: and as in Eq. (14). Then, the normalization for the local current (common to that of the local spatial axial-vector current up to discretisation effects), , is fixed by demanding that

(16) |

The values found on each ensemble and are given in Table 3.

There is an ambiguity in what mass to use on the right hand side of Eq. (16). We use the non-goldstone mass , but one could just as well replace this with since the difference is a discretisation effect. The meson mass difference is very small for heavy mesons Follana et al. (2007) and so we find the effect of changing the taste of meson mass used never exceeds 0.15% of throughout the range of ensembles and heavy masses that we use and has no impact on the continuum result.

We also remove tree-level mass-dependent discretisation effects coming from the wavefunction renormalisation Follana et al. (2007) by multiplying by a factor . This is derived in Monahan et al. (2013) as:

(17) | ||||

See also Bazavov et al. (2018a). is the tree-level pole mass in the HISQ action. It has an expansion in terms of the bare mass Follana et al. (2007)

(18) | ||||

fixes the Naik parameter Naik (1989) ( is the coefficient of the tree-level improvement term for the derivative) in the HISQ action when it is being used for heavy quarks Follana et al. (2007). is set to its tree-level value, removing the leading tree-level errors from the dispersion relation. As an expansion in it begins at Follana et al. (2007). To determine we use the closed form expression for it given in Monahan et al. (2013) and this can also be used along with Eq. (18) to evaluate . The pole condition can be used to show that the expansion of begins at as . The effect of is then very small, never exceeding . values on each ensemble for each are given in table 3.

Combining these normalizations with the lattice current from the 3-point fits, we find a value for the form factor at a given heavy mass and lattice spacing:

(19) |

Set | |||
---|---|---|---|

1 | 0.5 | 1.03178(57) | 0.99819 |

0.65 | 1.03740(58) | 0.99635 | |

0.8 | 1.04368(56) | 0.99305 | |

2 | 0.5 | 1.03184(47) | 0.99829 |

0.8 | 1.04390(39) | 0.99315 | |

3 | 0.427 | 1.0141(12) | 0.99931 |

0.525 | 1.0172(12) | 0.99859 | |

0.65 | 1.0214(12) | 0.99697 | |

0.8 | 1.0275(12) | 0.99367 | |

4 | 0.5 | 1.00896(44) | 0.99889 |

0.65 | 1.01363(49) | 0.99704 | |

0.8 | 1.01968(55) | 0.99375 |

### ii.5 Obtaining a Result at the Physical Point

We now discuss how we fit our results for the zero recoil form factor, , as a function of valence heavy quark mass, sea light quark mass and lattice spacing to obtain a result at the physical point where the heavy quark mass is that of the , the sea quark masses are physical and the lattice spacing is zero.

In summary, we fit our results for to the following form

(20) | |||||

The terms in the first line allow for dependence on the valence heavy quark and charm quark masses (with ) using input from HQET, to be discussed below. and account for discretisation and mass mistuning effects, also discussed below. The physical result is then .

#### ii.5.1 Dependence on the heavy valence quark mass

Our fit of the dependence is guided by HQET, which considers both the quark and the heavy quark of mass to be heavy here. In particular, for the parameter , HQET forbids terms of where can be or Luke (1990). The HQET expression for is then given by Falk and Neubert (1993); Mannel (1994):

(21) | ||||

where , and are (with possible mild dependence on whether the spectator quark is or ). accounts for ultraviolet matching between HQET and QCD, and has been computed to 2-loops in perturbative QCD Czarnecki (1996). It has mild dependence on through logarithms of ; at one-loop has explicit form Close et al. (1984)

(22) |

The coefficient of then varies from -0.66 to -0.29 across the range of from to , taking Bazavov et al. (2018b). The two-loop correction is small Czarnecki (1996). is then close to 1 and differs by a few percent across our range of .

Our calculation has results at multiple values of , and could therefore in principle provide information on the coefficients and of the -dependent terms in the HQET expansion. The charm quark mass is fixed to its physical value and so we cannot access the value of independent of a choice of at . The terms in round brackets in Eq. (21), multiplying , are all very small because of the suppression by heavy-quark masses. To constrain them tightly requires very precise data and, as we will see, we are not able to determine , or accurately with our results. It therefore does not make sense to attempt to compare them accurately to HQET expectations. To do so would require using an appropriate quark mass definition (since different definitions will move quark mass dependence between the term and the others in Eq. (21) ) and the two-loop expression for with appropriate value for (since logarithmic dependence of can be misinterpreted as part of a polynomial in ).

Instead we simply take an HQET-inspired form for the -dependence and set to 1, resulting in the first line of our fit form, Eq. (20). This is sufficient to test, through the results we obtain for , and using this expression, that the HQET expectation for the approximate size of these coefficients is fulfilled. We take priors on in our fit of .

We have several different proxies, derived from heavy meson masses, that we can take for the heavy quark mass that appears in in Eq. (20). We do not expect our physical result for to vary significantly depending on which meson mass we use, but the results for , and will vary because of different sub-leading terms in the relationship between meson and quark mass. The most obvious substitutions to use for the heavy quark mass are the mass of the pseudoscalar heavy-strange meson, , and half the mass of the pseudoscalar heavyonium meson, . We also tested using the quark mass in the minimal renormalon subtracted (MRS) scheme suggested in Brambilla et al. (2018). This takes

(23) |

where with for pseudoscalar mesons and for vectors. For this case we use parameters determined in Bazavov et al. (2018b) for the MRS scheme: , and . We take from Eq. (23) using our results for the mass of the pseudoscalar heavy-strange meson and from our results for the mass of the meson.

We take our central fit, for simplicity, from the result of using half the pseudoscalar heavyonium mass for and half the pseudoscalar charmonium mass for i.e. taking

(24) |

We test the stability of the fit results under the different choices discussed above in Section III.2.

#### ii.5.2 Mistuning of other quark masses

Our calculation has results for multiple different heavy quark masses on each gluon field configuration. The valence charm and strange quark masses, however, are tuned to their physical values. This is done by fixing the and meson masses to their physical values in a pure QCD world allowing, for example, for annihilation as discussed in Chakraborty et al. (2015). Any possible mistuning of the charm quark mass is accounted for in our fit function by the dependence on the charm quark mass that is included in the first line of Eq. (20). When the fit function is evaluated at the physical point we set from the physical mass.

The strange (valence and sea) and light (sea) mass mistunings are accounted for using the tuning in Chakraborty et al. (2015). For the strange quark, we define , where is given by

(25) |

is determined in lattice simulations from the masses of the pion and kaon Dowdall et al. (2013). The ratio then gives the fractional mistuning. The valence strange quark masses are very well tuned, but the sea strange quark masses less so.

We similarly account for mistuning of the masses of the (sea) light quarks by defining . We find from , leveraging the fact that the ratio of quark masses is regularization independent, and the ratio was calculated in Bazavov et al. (2018a):

(26) |

We set to divided by this ratio.

The full term we include to account for mistuning is then given by

(27) |

where , and are fit parameters with prior distributions . We neglect contributions since these are an order of magnitude smaller and are not resolved in the results of our lattice calculation.

The gluon field configurations that we use have in the sea. In the real world this is not true. We test the impact of possible isospin-breaking on our fits by testing for sensitivity to the sea light quark masses. We do this by changing the value up and down by the expected value for Tanabashi et al. (2018). We find the effect to be completely negligible in comparison to the other sources of error.

#### ii.5.3 Discretisation Effects

Discretisation effects in our lattice QCD results are accounted for following the methodology of McNeile et al. (2012b). We take

(28) |

The leading terms, with , allow for discretisation effects that are set by the heavy quark mass and also discretisation effects that are set by the charm quark mass (or indeed any other lighter scale that is independent of heavy quark mass). The terms allow for discretisation effects to vary as the heavy quark mass is varied, with being used here as a proxy for the heavy quark mass. are fit parameters with prior distributions . All discretisation effects are of even order by construction of the HISQ action Follana et al. (2007).

We tested the impact on the fit of including extra discretisation effects set by the scale but this made no difference (since such effects are much smaller than those already included by the terms). We also tested the effects of increasing the number of terms in each sum, but the final result remained unchanged.

#### ii.5.4 Finite-volume Effects

The finite volume effects in our lattice results are expected to be negligible, because we are working with heavy mesons that have no valence light quarks and no Zweig-allowed strong decay modes. Coupling to chiral loops or decay channels with pions that could produce significant finite-volume effects Laiho and Van de Water (2006) is therefore absent and we can safely ignore finite-volume effects here.

In section III.2 we detail the results of several tests of the stability of our final result under changes to the details of the fit.

## Iii Results and Discussion

### iii.1 Result for

The results of our correlation function fits (discussed in Section II.3) are given in Table 2. We tabulate values for at each heavy quark mass that we have used on each gluon field ensemble from Table 1. We also tabulate the meson masses needed to allow determination of at the physical point, using the fit form of Eq. 20.

The fit function of Eq. (20) is readily applied, giving a of 0.21 for 12 degrees of freedom. Figure 2 shows our results for along with the fit function at zero lattice spacing and physical , and quark masses as the grey band. Evaluating the fit function at the physical mass, as determined by , gives our final result

(29) |

Adding the statistical and systematic errors in quadrature, we find a total fractional error of . The error budget for this result is given in table 4. Note that we allow for an additional uncertainty in the physical value of the mass beyond the experimental uncertainty, since our lattice QCD results do not include the effect of annihilation and QED McNeile et al. (2012b). This has no effect, however, since the heavy quark mass dependence is so mild.

Source | % Fractional Error |
---|---|

Statistics | 1.06 |

0.73 | |

0.69 | |

mass mistuning | 0.20 |

Total | 1.45 |

Our total uncertainty is dominated by the statistical errors in our lattice results. The systematic error is dominated by that from the continuum extrapolation.

We include in Figure 2 the value from the only other lattice determination of