B-meson decay constants: a more complete picture from full lattice QCD

We extend the picture of $B$-meson decay constants obtained in lattice QCD beyond those of the $B$, $B_s$ and $B_c$ to give the first full lattice QCD results for the $B^*$, $B^*_s$ and $B^*_c$. We use improved NonRelativistic QCD for the valence $b$ quark and the Highly Improved Staggered Quark (HISQ) action for the lighter quarks on gluon field configurations that include the effect of $u/d$, $s$ and $c$ quarks in the sea with $u/d$ quark masses going down to physical values. For the ratio of vector to pseudoscalar decay constants, we find $f_{B^*}/f_B$ = 0.941(26), $f_{B^*_s}/f_{B_s}$ = 0.953(23) (both $2\sigma$ less than 1.0) and $f_{B^*_c}/f_{B_c}$ = 0.988(27). Taking correlated uncertainties into account we see clear indications that the ratio increases as the mass of the lighter quark increases. We compare our results to those using the HISQ formalism for all quarks and find good agreement both on decay constant values when the heaviest quark is a $b$ and on the dependence on the mass of the heaviest quark in the region of the $b$. Finally, we give an overview plot of decay constants for gold-plated mesons, the most complete picture of these hadronic parameters to date.


I. INTRODUCTION
Lattice QCD calculations are now an essential part of B physics phenomenology (see for example [1]), providing increasingly precise determinations of decay constants, form factors and mixing parameters needed, along with experiment, in the determination of Cabibbo-Kobayashi-Maskawa (CKM) matrix elements. As the constraints being provided by lattice QCD become more stringent it is increasingly important to expand the range of hadronic matrix elements being calculated to allow tests both against experiment where possible and/or against expectations from other approaches. Decay constants are particularly useful in this respect because they are single numbers expressing the amplitude for a meson to annihilate to a single particle (for example a W boson or a photon), encapsulating information about its internal structure. They are straightforwardly calculated in lattice QCD from the same hadron correlation functions being used to determine the hadron masses. The only additional complication is that normalisation of the appropriate operator for the meson creation/annihilation is required. In this way we can build up a tested and consistent 'big picture' of meson decay constants within which sit the results being used for CKM element determination.
To this end we determine here the decay constants that parameterise the amplitude to annihilate for the vector * christine.davies@glasgow.ac.uk † http://www.physics.gla.ac.uk/HPQCD mesons B * and B * s . These mesons are the partners of the B and B s whose weak decay matrix elements are critical to understanding heavy flavour physics. Decay modes of the B * and B * s are dominated by electromagnetic radiative decays [2] to B and B s , however, and so it is unlikely that processes in which the decay constant is the key hadronic parameter will be measured experimentally. The determination of the vector decay constants is nevertheless useful because the relationship with that of the pseudoscalar decay constant can be understood within the framework of Heavy Quark Effective Theory (HQET) and the decay constants appear in phenomenological analyses of the vector form factor for semileptonic decay processes for the pseudoscalar mesons (see [3] for a recent discussion of this).
Since vector and pseudoscalar heavy-light mesons differ only in their internal spin configuration, their decay constants might be expected to have rather similar values. The key question is then: by how much do they differ and which is larger? A recent review [4] showed tension between the results for the ratio of the B * to B decay constants from QCD sum rules and from lattice QCD. The lattice QCD results used u/d quarks (only) in the sea and obtained results for mesons containing b quarks from an interpolation between results for quarks close to the c mass and the static (infinite mass) limit [3]. This gave a result for the ratio greater than 1 whereas the QCD sum rules approach quoted preferred a value less than 1.
The results we give here build on our state-of-the-art calculation of the B and B s decay constants [5] using an improved NonRelativistic QCD (NRQCD) formulation [6] that allows us to work to high accuracy directly at the b quark mass. We also use lattice QCD gluon field configurations that have the most realistic QCD vacuum to date, include u/d, s and c quarks in the sea (using the Highly Improved Staggered Quark formalism [7]) with the u/d quark mass taking values down to the physical value. We therefore avoid significant systematic errors from extrapolations in the u/d quark mass. We are able to give results for the ratio of vector to pseudoscalar decay constants for both the B s and the B and the SU(3)breaking ratio of these ratios. We find clearly that the vector decay constant is smaller than the pseudoscalar decay constant in both cases.
We also give results for the decay constant of the B c meson and its vector partner the B * c . The B c has been seen experimentally only relatively recently [2] and is interesting because it can be viewed both as a heavy-heavy meson and as a heavy-light meson. Here we compare its ratio of vector to pseudoscalar decay constants to that of the B and B s , and find the ratio is significantly larger, now being very close to 1.
The decay constant of the B c is quite different from that of the B and B s , being nearly double their size. The B c decay constant can be used to predict its partial width for leptonic decay that may be observed in the future. We determine this decay constant here using NRQCD b quarks and HISQ c quarks and compare to our previous result [8] that used the HISQ action for both quarks and mapped out the behaviour of a range of decay constants for valence heavy quarks in the region between the c quark mass and the b quark mass. Since the HISQ action is fully relativistic this is a good test of our understanding of systematic errors in lattice QCD, and confirmation of how well improved actions work.
In a further study of this point we go on to look at the dependence of decay constants on the valence heavy quark mass using quark masses lighter than that of the b in the NRQCD action. This enables us to compare both the value of specific decay constants and the dependence on the heavy quark mass with that from using the relativistic HISQ action. We also demonstrate the consistency of our results for the ratio of vector to pseudoscalar decay constants for the B s meson here to our earlier result for the same ratio for the D s meson [9] using HISQ quarks.
A very consistent picture thus emerges from both a nonrelativistic and a relativistic approach to heavy quarks within lattice QCD. Both approaches are the result of several stages of improvement to reduce discretisation errors and other systematic uncertainties to a low level, important for making a detailed comparison.
We begin by outlining the methods used in our lattice calculation, which follow [5,10]. Section III A gives results for the decay constants of the B * s and B * and their comparison, and then section III D gives results for the decay constants of the B c and the B * c . Section III E works with quarks lighter than b to demonstrate the heavy quark mass dependence of the decay constants and compare to our earlier results using the HISQ formalism for b quarks. Section IV compares our results for vector meson decay constants to those of earlier determinations using other methods, including HQET arguments, and shows how the B c fits in between results for heavyonium and heavy-light mesons. Section V gives our conclusions, including the promised 'big picture' for the decay constants of gold-plated mesons from lattice QCD, the most complete picture of these hadronic parameters to date.

II. LATTICE CALCULATION
Since the first lattice NRQCD calculations were done for heavy-light mesons [11], huge improvements have been made. The current state-of-the-art [5,10] uses an improved NRQCD action for the heavy quark coupled to a HISQ light quark on gluon configurations that include an improved gluon action and HISQ sea quarks. Here we extend these calculations to include the decay constants of vector heavy-light mesons.
The gluon field configurations that we use were generated by the MILC collaboration [12,13]. These are n f = 2+1+1 configurations that include the effect of light (up/down), strange and charm quarks in the sea with the HISQ action [7,14] and a Symanzik improved gluon action with coefficients correct through O(α s a 2 , n f α s a 2 ) [15]. The lattice spacing values that we use range from a = 0.15 fm to a = 0.09 fm. The configurations have well-tuned sea strange quark masses and sea light quark masses (m u = m d = m l ) with ratios to the strange mass from m l /m s = 0.2 down to the value that corresponds to the experimental π meson mass of m l /m s = 1/27.4 [16].
In [6] we accurately determined the lattice spacings using the mass difference of the Υ and Υ mesons using the same NRQCD action for the b quark as we use here. The details of each ensemble, including the lattice spacing, sea quark masses and spatial volumes, are given in table I. All ensembles were fixed to Coulomb gauge.

A. NRQCD valence quarks
We use improved NRQCD for the b quark, which takes advantage of the nonrelativistic nature of the b quark within its bound states for very good control of discretisation uncertainties. It has the advantage that the same action can be used for both bottomonium and B-meson calculations so that tuning of the b-quark mass and determination of the lattice spacing can be done using bottomonium and there are no new parameters to be tuned at all for B-mesons. b-quark propagators are calculated with a very fast forward evolution using NRQCD and this also means that high statistics can be accumulated readily for precise results. The action used here builds on the standard NRQCD action [19] accurate through v 4 in the heavy quark velocity v (using power-counting terminology for bottomonium) by including one loop radiative corrections to many of the v 4 coefficients [6,20].  [12,13]. β is the bare gauge coupling, aΥ is the lattice spacing as determined by the Υ(2S −1S) splitting in [6], where the three errors are statistics, NRQCD systematics and experiment. Column 4 gives the corresponding values of αs used in renormalisation factors. This is taken as αV (n f = 4, q = 2/a) and determined from [17,18]. am l , ams and amc are the sea quark masses in lattice units. We also give, in column 8, values for δxsea, the fractional difference in the sum of the light quark masses from their physical values. δxsea is defined as (2m l + ms)/(2m l,phys + m s,phys ) − 1, using values of m s,phys from [6] and ms/m l = 27.4 [16]. L/a × T /a gives the spatial and temporal extent of the lattices and n cfg is the number of configurations in each ensemble. 16 time sources were used for the valence quark propagators on each configuration for increased statistics. Sets 1, 2 and 3 will be referred to in the text as "very coarse", 4, 5 and 6 as "coarse" and 7 as "fine". Sets 3 and 6 include light sea quarks with their physical masses. We studied the effect of these improvements on the bottomonium spectrum in [6,21,22] and in B, B s and B c meson masses in [10]. The NRQCD Hamiltonian we use is given by [19]: Here ∇ is the symmetric lattice derivative and ∆ (2) and ∆ (4) the lattice discretization of the continuum i D 2 i and i D 4 i respectively. am b is the bare b quark mass. The parameter n has no physical significance, but is included for numerical stability of high momentum modes. We take the value n = 4 here in all cases.Ẽ andB are the chromoelectric and chromomagnetic fields calculated from an improved clover term [23]. TheB andẼ are made anti-hermitian but not explicitly traceless, to match the perturbative calculations done using this action.
The coefficients c i in the action are unity at tree level but radiative corrections cause them to depend on am b at higher orders in α s . These were calculated for the relevant b quark masses using lattice perturbation theory in [6,20] and the values used in this paper are given in Table II. Including the one-loop radiative corrections to c 4 is particularly important here, since this coefficient controls the hyperfine splitting between the vector and pseudoscalar states. We showed in [10] that improving c 4 leads to accurate results for b-light hyperfine splittings in keeping with the results of [6] for bottomonium. Most of the correlators we use here for determining the vector heavy-light meson decay constants come from the same calculation as that of [10]. The tuning of the b quark mass on these ensembles was discussed in [6]. We use the spin-averaged kinetic mass of the Υ and η b and tune this to an experimental value of 9.445(2) GeV. This allows for electromagnetism and η b annihilation effects missing from our calculation [24]. Note that we no longer have to apply a shift for missing charm quarks in the sea [24]. The values used in this calculation are the same as those in [5,10] and given in table III along with other parameters.

B. HISQ valence quarks
For the u/d, s and c valence quarks in our calculation we use the same HISQ action as for the sea quarks. The advantage of using HISQ is that am q discretisation errors are under sufficient control that it can be used both for light and for c quarks [7,14,26]. The HISQ action is also numerically inexpensive which means we are able TABLE III. Parameters used in the NRQCD action. am b is the bare b quark mass and u0L the Landau link tadpoleimprovement factor used in the NRQCD action [25]. δx b gives the fractional mistuning in the b quark mass ((am b − am b,phys )/am b,phys ) obtained from the determination of the spin-averaged kinetic mass of the Υ and η b [6], when this has a magnitude larger than 0.5%. The column asm gives the size parameters of the quark smearing functions (see section II C), which take the form exp(−r/asm). asm is kept approximately constant in physical units a .
Set to perform a very high statistics calculation to combat the signal to noise ratio problems that arise in simulating B-mesons. For example, we use 16 time sources for both NRQCD and HISQ propagators on each configuration, so we are typically generating 16,000 correlators per ensemble.
The masses used on each ensemble are given in table IV. Again these are the same as in [5,10]. In [6] accurate strange quark masses were determined for each ensemble, tuned from the mass of the η s meson, a pseudoscalar ss which can be prevented from mixing with other states on the lattice so that its mass can be determined very accurately [27]. Using experimental K and π meson masses we found M ηs = 0.6893 (12) GeV (see also [28]). The values of am val s in table IV correspond to these tuned values. The light valence quarks are taken to have the same masses as in the sea.
Charm quark masses are tuned by matching the mass of the η c to experiment. The experimental value is shifted by 2.6 MeV for missing electromagnetic effects and 2.4 MeV for not allowing it to annihilate to gluons, giving 2.985(3) GeV [27]. The Naik term in the action is not negligible for charm quarks and we use the tree level formula given in [26]; the values appropriate to our masses are given in table IV.

C. NRQCD-HISQ correlators
The NRQCD b and HISQ u/d, s or c light quark propagators are combined into a meson correlation function in a straightforward way. Since staggered quarks have no spin index, staggered quark propagators must first be converted to 4-component 'naive' propagators so that they can be combined with quark propagators from other formalism such as NRQCD. To do this, the 4x4 'staggering matrix' Ω(x) = 4 µ=0 (γ µ ) xµ that was used to convert the naive quark action into the staggered quark action has to be applied to the staggered quark propagator at each end to 'undo' the transformation [29]. The spin and colour components of the naive propagator and the NRQCD propagator can then be tied up at source and sink with appropriate γ matrices (taking appropriate 2×2 blocks since the NRQCD propagator is 2-component) to form a pseudoscalar or vector meson correlator. We sum over the spatial sites on the sink time-slice to project onto zero spatial momentum in all cases.
One complication is that 'random-wall' sources (i.e. a set of U(1) random numbers over a timeslice) are used for the light quark propagators to improve statistical accuracy in our light meson calculations (see, for example, [28]). When these propagators are tied together the result is equivalent to having a delta function source at each point on the time-slice. As well as the convenience, statistical accuracy is also improved by re-using these propagators in our heavy-light meson calculations. The source for the NRQCD propagators must then use the same random numbers on the same source time-slice and in addition must also include a spin trace over the staggering matrix and appropriate gamma matrix for a pseudoscalar or vector meson [24], i.e. there is a separate NRQCD propagator for each meson that will be made. Combining these NRQCD propagators with the light quark propagators is then equivalent to having a delta function source at each point on the timeslice, as for the light meson case.
A further numerical improvement is to make smeared sources for the NRQCD propagators by convoluting a smearing function with the random-wall source as above. Suitably chosen smearing functions can improve the overlap of the correlator with different states in the spectrum, and this is particularly important for fits to extract radially excited energies [6]. Here we use it to improve the determination of ground-state properties by improving the overlap with the ground-state at early times before the signal/noise ratio has degraded significantly. For each ensemble we then use a local source and 2 smeared sources. The smearing functions take the form exp(−r/a sm ) as a function of radial distance, with two different radial sizes, a sm , on each ensemble as given in Table III. Propagators were calculated, and meson correlators obtained, using 16 time sources on each configuration. The calculation was also repeated with the heavy quark propagating in the opposite time direction. All correlators from the same configuration were binned together to avoid underestimating the statistical errors. When each smearing is used at source and sink this gives a 3 × 3 matrix of correlation functions on each ensemble. In addition we have a 3-vector of correlation functions from using each smearing at the source and a relativistic current correction operator at the sink, to be discussed below.
Meson energies and amplitudes are extracted from the meson correlation functions using a simultaneous multiexponential Bayesian fit [30] as a function of time separation between source and sink to the form Here i, j label the smearing (or current correction operator) included in the correlator and k labels the set of energy levels for states appearing in the correlator. Here we are concentrating on the properties of the groundstate, k = 0. k labels a set of opposite-parity states that appear with an oscillating behaviour in time as a result of using staggered quarks. Energies for groundstates, radially excited states and oscillating states were extracted from these correlation functions in [10]. Here we use similar fits to determine ground-state amplitudes and thereby decay constants. The fits are straightforward and follow the same pattern in all cases. We fit the pseudoscalar and vector correlators simultaneously for each pair i.e. B and B * , B s and B * s and B c and B * c . That enables us to extract a correlated ratio of amplitudes that we need for the ratio of decay constants. We take a prior on the ground-state energy determined from effective mass plots, with a width of 300 MeV. The prior on the lowest oscillating state is taken to be 400 MeV higher than the ground-state with a width of 300 MeV. The prior on the energy splittings in both the oscillating and non-oscillating sectors, E n+1 − E n , is taken as 600(300) MeV and the priors on the amplitudes as 0.1(2.0). The fits include points from t min to t max , close to half the temporal extent of the lattice. t min is taken from 6 − 8 for B and B s fits and t max is taken as 18 on the very coarse lattices, 28 on coarse and 40 on fine. For B c we use time ranges 12 − 24 on very coarse, 8 − 21 on coarse and 10 − 30 on fine. In all cases we have good fits that reach stable ground-state parameters quickly. We take results from fits that use N exp = 4.

D. Determining decay constants
Meson decay constants, f , are hadronic parameters defined from the matrix element of the local current that annihilates the meson (coupling for example to a W boson). For mesons at rest: Here H is one of the pseudoscalar mesons, B q for q = l, s, c. These matrix elements depend on the QCD interactions that keep the quark and antiquark bound inside the meson. They can be calculated directly from the amplitudes obtained from fits to the meson correlators, provided that we can accurately represent the continuum QCD currents, J A0 and J Vi on the lattice. The representation of these currents when combining a lattice NRQCD b-quark with a light quark is discussed most recently in [5,31]. For the temporal axial current whose matrix element gives the pseudoscalar decay constant, we determine matrix elements on the lattice made from light quark fields Ψ q and NRQCD field Ψ Q of: A0 is the leading term in a nonrelativistic expansion of the current operator, and J A0 is simply the operator that corresponds to our local sources for the b quark described in Section II C. Thus the matrix element of J (0) A0 between the vacuum and the ground-state meson is obtained directly from the amplitude of this operator from our fit function, i.e. b loc,0 from eq. (3). By inserting a complete set of states with standard normalisation into the pseudoscalar meson correlation function we have Similarly, by inserting the operator J A0 at the sink for meson correlators made from the three different sources that we use we can determine an amplitude for this operator in the ground-state The way in which J A0 can be combined into an accurate representation of J A0 from full QCD is described in [5]. Here, for most of our results, we will use an expression that is slightly less accurate than that in [5]. We take: Thus, we can combine matrix elements above to obtain: up to sources of uncertainty that will be discussed in the appropriate subsections of Section III. Note that the decay constant appears naturally multiplied by the square root of the meson mass in these expressions. Analogous expressions are used for the vector current case, using amplitudes from the vector meson correlator fits.
Systematic errors are reduced by working with the ratio of vector to pseudoscalar meson decay constants (multiplied by the ratio of the square root of the masses). Hence we define the quantity R q for meson B q , determined from: δz is z Vi − z A0 and will be tabulated, along with the results, in Section III. This expression is accurate up to missing α 2 s pieces of the overall renormalisation factor (i.e. in the term (1 + δz · α s )) and missing additional α s renormalisation factors for the sub-leading current contributions (that would appear multiplying Φ (1) A0 for example). These sources of systematic uncertainty will be estimated in Section III and included in our final error budgets.
We will first calculate R s as the 'calibration' ratio of vector to pseudoscalar decay constants. It is convenient subsequently to calculate ratios of R l and R c to R s . Some systematic errors cancel in these ratios of ratios, allowing us to obtain a more accurate picture of how much the ratio of vector to pseudoscalar heavy-light meson decay constants depends on the light quark mass.

A. b-light correlators
Correlators for B s , B * s , B l and B * l are fitted as described in Section II and results for the ground-state amplitudes of leading, J (0) , and sub-leading, J (1) , currents are tabulated in Tables V and VI respectively. In Table VIII we also tabulate for both B s and B l the ratio of the sum of the amplitudes that make up the NRQCD TABLE V. Amplitudes for J (0) and J (1) for temporal axial and vector currents between the vacuum and the Bs and B * s mesons respectively, extracted from correlator fits and multiplied by √ 2 in accordance with eq. (9). Results are in lattice units and the errors given are statistical/fit errors only. Results for the Bs were previously given in [5]. Results here differ slightly because the fits included both vector and pseudoscalar correlators in a simultaneous fit and also incorporated more correlators that included J (1) amplitudes. (1) TABLE VI. Amplitudes for J (0) and J (1) for temporal axial and vector currents between the vacuum and the B l and B * l mesons respectively, extracted from correlator fits and multiplied by √ 2 in accordance with eq. (9). l denotes a u or d quark, taken here to have the same mass. Results are in lattice units and the errors given are statistical/fit errors only. Results for the B l were previously given in [5]. Results here differ slightly for reasons given in the caption to Table V.  (2) TABLE VII. Coefficients zA 0 and zV i needed for the oneloop renormalisation factor for the pseudoscalar and vector decay constants respectively for the values of m b a used on the different ensembles. z is constructed from results given for the appropriate NRQCD bare masses and massless HISQ quarks in [31] as z = ρ0 −ζ10. In [5] zA 0 is called z0. Column 4 gives δz which is the difference between zV i and zA 0 . Column 5 gives the corresponding values of δz for the case where only the leading-order NRQCD currents (J Bs ) eq. (11). Column 3 gives the equivalent quantity for the B * l /B l mesons. Column 4, R LO s , gives renormalised ratio from eq (12) but including only the leading-order NRQCD currents, J (0) . Finally column 5 gives the renormalised ratio of amplitudes including the current corrections. These numbers are determined from eq. (10) and plotted as the points in Figure 1.  (12)  7 0.9775 (14) 0.9734 (23) 0.8916 (16) 0.9678 (17) vector and temporal axial currents (without any renormalisation factors) defined as: These ratios are determined directly from the fits, including the correlations between the fitted amplitudes for vector and pseudoscalar mesons, and therefore have smaller statistical errors than determining them naively from the results in Tables V and VI. The z factors needed to multiply α s in the one-loop renormalisation for the temporal axial and spatial vector currents are given in Table VII. These are calculated for massless HISQ light quarks and the values of am b in the NRQCD action used on each of the ensembles. The fact that these z coefficients are very small was already noted in [5]. This means that renormalisation factors to the continuum current are close to 1 1 .
From Table V and VI it is immediately clear that the leading order amplitudes, Φ (0) , show a difference between vector and pseudoscalar mesons with the vector result being smaller than the pseudoscalar. This difference is largely down to the 'hyperfine' interaction in the NRQCD Hamiltonian (the term with coefficient c 4 in eq. (2)). The Tables make clear that the impact of this interaction is 1 Note that we do not need an initial nonperturbative step to achieve this, as is used by the Fermilab Lattice/MILC Collaboration [32]. That step is largely required to remove large but generic renormalisation factors associated with the clover action and has been tested nonperturbatively in [33]. to lower the ratio of vector to pseudoscalar decay constants. This effect agrees in sign with that seen in an earlier lattice NRQCD analysis of the impact of different relativistic corrections on heavy-light meson decay constants [34]. It also agrees with early estimates using HQET and QCD sum rules [35].
From the tables it is also clear that the relativistic current correction matrix element, Φ (1) , has opposite sign for the vector and pseudoscalar cases, being positive for the vector and negative for the pseudoscalar. The impact of these corrections is then to raise the ratio of vector to pseudoscalar decay constants. This sign, and the fact that the pseudoscalar J (1) matrix element is approximately three times that of the vector, agrees with HQET expectations [35,36] and earlier lattice NRQCD analyses [34].
Simply dividing the current correction matrix element, Φ (1) , by Φ (0) gives naively a relative contribution to the amplitude from the relativistic current corrections of size -(8-10)% for the pseudoscalar and +3% for the vector. This does not take into account the fact that the addition of the relativistic current correction J (1) which appears at tree-level changes the overall renormalisation of the lattice NRQCD current at O(α s ) because radiative corrections to J (1) can look like J (0) [31]. Therefore to determine more accurately the effect of the relativistic current corrections we have to compare the renormalised result with and without the inclusion of the J (1) current correction. This is done in Table VIII in which we compare the two results for the ratio of f √ M for B * s and B s . The righthand column, denoted R s , is the full result obtained from eq. (10) using the values of δz from Table VII. The values for α s used in that expression are taken in V scheme at the scale 2/a, where a is the lattice spacing on that ensemble, and are given in Table I. The column denoted R LO s gives results using only the leading-order currents and with δz LO values given in Table VII and the same values of α s . Note the difference between δz LO and δz. Both coefficients are small, but they have opposite sign. This then compensates to some extent for the effect of the current corrections and means that, comparing R s and R LO s in Table VIII we see now that the total effect of the current correction terms in the ratio amounts to 7-8%, somewhat less than the naive estimate of 12-15%. There is of course an uncertainty on this estimate coming from missing α 2 s terms in the renormalisation. A similar procedure would be needed to estimate accurately the effect of the hyperfine term on the ratio R q . However, because the hyperfine interaction is embedded in the NRQCD Hamiltonian it is automatically included in the perturbative matching calculation for the NRQCD currents and we do not have the z coefficients without the hyperfine term included.
B. fB * s Figure 1 plots the full results for R s , the ratio of f √ M for the B * s and B s mesons, obtained from eq. (10) and given as column 5 of Table VIII. Statistical errors in R s are small, less than 0.5%, so we see that the value for the ratio is clearly less than 1 and the dependence on the lattice spacing is small, but clear and unambiguous. To derive a physical result we need to fit this dependence, allowing for other systematic uncertainties from lattice QCD.
The key sources of systematic error that need to be allowed for, by inclusion in our fit function, are: • Matching uncertainties -α 2 s . The missing α 2 s coefficient in the overall renormalisation factor for the ratio of amplitudes of the NRQCD currents is potentially the largest source of uncertainty here. We can allow for this by simply taking a fractional error which is α 2 s ≈ 0.1 times a value for this coefficient. However, the value of α 2 s changes with the lattice spacing and the coefficient may also depend on am b , as the known one-loop coefficient δz does, see Table VII. Thus a better estimate is obtained by incorporating a factor to take account of this missing term into the fit. We write the factor as (1 + cα 2 s ) and take c to have the form where c 1 sets the overall allowed size of the coefficient and the δx m terms allow for dependence on am b . δx m = (am b −2.7)/1.5 varies from -0.5 to 0.5 over the range of am b values we use here [6].
• Matching uncertainties -α s Λ/m b . We must also allow for missing α s terms that alter the normalisation of the relativistic current corrections within the NRQCD current and/or include the matrix elements of additional current corrections that only appear first at O(α s Λ/m b ). Such corrections were included in our determination of f Bs and f B in [5] since they are known for the temporal axial current for massless HISQ quarks [31]. They will be discussed further in Subsection III E but here we must include an uncertainty for the fact that they are missing in our ratio. For this we can include an additional term in the factor described above of the form dα s Λ/m b where d has an expansion in powers of δx m of the same form as c above. Here we can take Λ/m b to be 0.08, the size of the relativistic current corrections as determined above.
• Matching uncertainties -(Λ/m b ) 2 . Further current corrections at the next order in the relativistic expansion would appear at (Λ/m b ) 2 . Since we have no information about these we do not include them in the fit but take an additional uncertainty of (0.1) 2 = 1% to account for them.
The improved NRQCD Hamiltonian that we use (eq. (2)) is accurate through O(α s Λ/m b ) in the context of heavy-light power-counting. Thus the hyperfine interaction that contributes to R s is accurate through this order, which is to a higher order than the matching uncertainties discussed above. Errors from the NRQCD Hamiltonian are then smaller than, and are effectively included in, the matching uncertainties already discussed. Likewise missing terms in the NRQCD Hamiltonian are at even higher order, O(1/m 3 b ) [19]. • Discretisation uncertainties. These can come from the gluon action, the HISQ action and the NRQCD action. However, most discretisation uncertainties will cancel between vector and pseudoscalar mesons since the difference between them is a spindependent effect and hence suppressed by Λ/m b . This is clear from Figure 1 which shows very little dependence on a. In all three actions discretisation errors appear as even powers of a. We therefore include a factor (1+(Λ/m b ) j e j (Λa) 2j ) to allow for these uncertainties in the fit. We take a value 0.2 for Λ/m b here to be conservative. For e 1 we allow for dependence on am b coming from the NRQCD action as discussed above for the coefficients c and d.
• Tuning uncertainties -valence quark masses. Our valence masses are tuned very accurately (to 1%) but we allow for effects of mistuning. For the s quark mass these will be negligible since, as we show below, the difference between R s and R l is very small. Mistuning of the b quark mass will affect R s through the hyperfine interaction and the size of the current correction matrix elements, i.e. through a term of the form (Λ/m b )δm b /m b . We therefore allow for a term of this form in the factor that includes discretisation effects above. We determine δm b from the physical values for m b given on each ensemble in [6,10] and these are tabulated in Table III. The largest value of δm b /m b is 1.3% on set 4.
• Tuning uncertainties -sea quark masses. Our results include values on ensembles of gluon field configurations at a variety of values of the u/d quark mass in the sea, varying from 0.2m s down to the physical point. The s and c quark masses in the sea are well-tuned. Dependence on the sea quark masses is very small, as is clear from Figure 1.
We therefore include a simple linear dependence on the sea quark masses, as might be expected from leading-order chiral perturbation theory. This dependence takes the form gδm sea /(10m sea,phys ) where we include u/d and s quarks in m sea and the factor of 10 is a convenient way to introduce the chiral scale of 1 GeV. δm sea = (2m l + m s ) − (2m l,phys + m s,phys ) and is obtained using values for m s,phys given in [6]. We take m l,phys /m s,phys = 1/27.4 [16]. Values for δx sea ≡ δm sea /m sea,phys on each ensemble are given in Table I.
• Uncertainties in the value of the lattice spacing.
Since we are determining a dimensionless ratio of decay constants, uncertainties in the value of the lattice spacing only enter indirectly through the uncertainty in tuning the quark masses. As discussed above the tuning of m b affects the size of the relativistic correction terms that affect the vector/pseudoscalar ratio. We have a 1% uncertainty in our lattice spacing values, largely correlated between the ensembles and so we add an additional overall uncertainty of 0.2 × 0.01 = 0.2% to allow for this. The factor of 0.2 is a conservative estimate for the size of relativistic corrections.
Putting the features above together we arrive at a fit form for R s as a function of a and quark masses as: In dividing by F 2 we follow the convention that we used at O(α s ) in eq. (10) of writing the renormalisation as a multiplicative factor. Thus if F 2 were instead known, rather than fitted, the raw results would be multiplied by this correction factor along with the factor at O(α s ). We use a Bayesian fitting approach [30] to implement the fit function of eq. (13). Priors on all of the coefficients are taken as 0.0(1.0) except for c 1 , which is taken as 0.0(0.2). This allows for an α 2 s coefficient in the overall renormalisation factor that is twice as large as the largest seen at O(α s ) (see δz values in Table VII). The prior on the physical value, R s,phys , is taken as 1.0(0.2).
Applying this fit function to our results gives a χ 2 /dof = 0.13 and a physical result for R s of 0.957(23), when we include the uncertainty from missing higher order current corrections and the lattice spacing. The error budget from the fit is laid out in Table IX. As expected the uncertainty is dominated by that from current matching, although the fit has constrained this uncertainty to be a bit smaller than the naive expectation. The physical value, along with the total error, is plotted as a grey band on Figure 1. R s is the ratio of decay constants multiplied by the square root of the masses. We can convert this to a ratio for the decay constants using the square root of the experimental ratio of the masses of 1.0045(2) [2]. This then gives This is 2σ below 1. To analyse the corresponding ratio, R l , for the B/B * it is convenient to take the ratio to R s . Table VI gives our results for the B and B * amplitudes and Table VIII gives the ratio of the sum of amplitudes for J (0) and J (1) for vector and pseudoscalar. These results include the correlations between the vector and pseudoscalar meson correlators from the simultaneous fit. The results for B * l /B l are very similar, not surprisingly, to those for B s and B * s . The statistical errors are significantly larger, however, as is expected when the light quark mass is reduced [26]. The renormalisation factor (eq. (10)) for R l is the same as that for R s (since mass effects for light quarks are negligible in the matching) and so the renormalisation cancels in the ratio R l /R s . We can therefore simply determine R l /R s from the ratio of the first two columns in Table VIII: Figure 2 shows our results for R l /R s on each ensemble, plotted against the light quark mass in units of the physical s quark mass taken from [6]. The results are very close to the value 1.0, but show a small downward trend as the light quark mass falls towards its physical value. There is no significant dependence on the lattice spacing.
To fit the dependence of R l /R s and extract a physical result, we use much of the same fit function as that given for R s in eq. (13). The two key differences are that the overall renormalisation factor now cancels, so that we can drop the factor cα 2 s from F 2 , and that we now want to include a fitted dependence on m l . We can also use the known similarity of B l and B s to constrain the fit further. For example, we know that decay constants for heavylight and heavy-strange mesons differ by about 20% [5,14]. This is in fact a very strong result, still true even when the light/strange quark is accompanied by a light or strange quark (see, for example, [28] for π, K and η s results). The dα s term in eq. (13) takes account of missing radiative corrections to the sub-leading currents, J (1) . We retain that term here but multiply its coefficient by 0.2 to allow for strange/light differences in the matrix element for J (1) . We reduce the coefficient 0.2 (allowed for the size of Λ/m b ) in front of discretisation errors and m b tuning terms by a further factor of 0.2 for the same reason. Finally, we include an additional term in the fit to allow for dependence on the light valence mass, since we have results for a variety of m l values. For this we include a term in F 1 of the form h(m l /(10m s,phys )). The factor of 10 once again is used to convert m s into the chiral scale of 1 GeV. The prior on h is taken as 0.0(1.0). Since this term is already largely covered by including a term to allow for sea quark mass dependence, it has very little impact.
The fit has a χ 2 /dof of 0.23 and gives a physical result: Since R l /R s measures both SU(3)-breaking and spinbreaking effects in heavy-light meson decay constants we expect a result very close to 1.0, but our value enables us to constrain any difference from 1.0 to a few percent. It also indicates a tendency for the vector to pseudoscalar decay ratio to be smaller for B * /B than for B * s /B s , i.e. for the ratio to increase with increasing mass, m q , in the B q meson. We will return to this in section III D when comparing results for B c mesons. A full error budget for R l /R s is given in Table IX. Combining our result for R l /R s with our earlier result for R s gives R l = 0.945 (26). Combining with the experimental value for the square root of the ratio of the meson masses, 1.0043 [2], we obtain which is more than 2σ below 1.
D. fB c and fB * c B c and B * c meson correlation functions are calculated from NRQCD b and HISQ c propagators in exactly the same way as those described for NRQCD b and HISQ s or l propagators in subsection III A. We do not include the full set of ensembles used for the lighter HISQ quark mass calculations since experience has shown very little sea quark mass dependence for heavy meson correlators that do not include valence light quarks [10,24]. We thus include ensembles at two different values of the sea u/d quark mass for very coarse and coarse sets rather than three. The meson correlation functions are fit simultaneously so that correlations between them can be included in the determination of the ratio of amplitudes needed for the ratio of decay constants.
The results for the matrix elements, Φ, of the leading, J (0) , and subleading, J (1) , pieces of the temporal TABLE X. Amplitudes for J (0) and J (1) for temporal axial and vector currents between the vacuum and the Bc and B * c mesons respectively, extracted from correlator fits and multiplied by √ 2 in accordance with eq. (9). Results are in lattice units and the errors given are statistical/fit errors only. The ground-state energies determined from the fits agree with those given in [10] and we do not repeat them here.   (18)). zA 0 ,c is constructed from results given for the appropriate NRQCD bare masses and massive HISQ quarks with the appropriate values of mca in [31] as zc = η0 − τ10. z1 and z2 are results for massless HISQ quarks [5,31]. The values of αs used with these z coefficients are given in Table I. axial and spatial vector currents are given in Table X. We first discuss combining the results for the temporal axial current into a value for the decay constant of the pseudoscalar B c meson. We will use a formula [5] which is somewhat more accurate than that used in eq. (9): Bc . z 1 α s is an additional radiative correction to the subleading current J (1) . z 2 α s multiplies an additional subleading current which has the same matrix element as J (1) and so does not need to be separately calculated.
The z coefficients now have to be calculated for massive HISQ quarks with a mass in lattice units corresponding to our values for am c on the different ensembles. This has been done for z A0,c and the values are given in Table XI. They differ slightly from those for massless HISQ quarks in Table VII but are still very much less than 1. The z 1 and z 2 coefficients have only been calculated for massless HISQ quarks and these are also given in Table XI. There is then a systematic error in our formula of eq. (18) as a result of using the massless z 1 and z 2 coefficients and we will allow for that in our error budget along with systematic errors from unknown higher order terms in the overall renormalisation factor. The results obtained from applying eq. (18) are plotted in Figure 3 as a function of lattice spacing. We see, as expected, very little change between ensembles with similar lattice spacings but different sea u/d quark masses. To determine a physical value for the decay constant we fit the results to a functional form that includes allowance for systematic errors in the lattice QCD calculation.
The systematic errors have the same sources as those discussed for R s in section III A and we will use the same fit form as that given in eq. (13) and we reproduce that below as eq. (19) with the modifications appropriate here. As in section III A, the major source of uncertainty here comes from missing higher order terms in the matching of the NRQCD-HISQ current to continuum QCD. This is taken account of in eq. (19), as before, by the term F 2 (α s ) which includes an α 2 s term with coefficient c in the overall renormalisation factor and a term with coefficient d that allows for systematic errors in the α s corrections to the J (1) current contribution included in eq. (18) from the fact that z 1 and z 2 are taken for massless HISQ quarks. Given the values we have for z A0,c and the dependence on am c seen in that coefficient, we do not expect coefficients c and d to be large and we take priors on their fit values of 0.0(0.2).
From Fig. 3 we see significant lattice spacing dependence in the results and we must allow both for regular lattice spacing dependence and that coming from the NRQCD action. This dependence is included in factor F 1 . The regular lattice spacing coming the HISQ action can have a scale set by m c in this case and we expect that to dominate. We take m c to be 1 GeV here. The analysis of discretisation errors for c quarks in the HISQ action [7] shows that the dependence comes from terms suppressed by powers of the velocity of the c quark. Since v 2 c ≈ 0.5 in a B c [37] we include a factor of 0.5 in front of the terms allowing for discretisation errors. We must also allow for dependence on the u/d quark mass in the sea, as before, and for mistuning of the b quark mass. For mistuning of the b quark mass we allow a conservative factor of 0.3 based on the variation in decay constants between heavyonium mesons (see Fig. 8).
Our fit function is: We take prior values on all coefficients to be 0.0(1.0) except for the physical value on which we take 1.0(2), c and d, on which we take 0.0(2) and e 1 on which we take 0.0(3) (since it is O(α s )). The error on the plotted values in Figure 3 is dominated by the uncertainty in the value of the lattice spacing (given in Table I). In doing the fit we allow for half of this error to be correlated between ensembles (since it comes from systematic uncertainties from the NRQCD calculation used to fix the lattice spacing [6]) and half to be uncorrelated.
The fit gives a χ 2 /dof of 0.11 and a physical value for f √ M for the B c of 1.087(37) (GeV) 3/2 . The 3.4% uncertainty is split between 1.2% from matching and 3.2% from other sources, dominated by lattice spacing uncertainties and discretisation errors. We have checked that missing out z 1 and z 2 from eq. (18) and allowing a larger prior of 0.0(1.0) on the coefficient d in eq. (19) gives a physical result with almost the same central value and uncertainty.
Our physical result is plotted as a grey band in Figure 3. It agrees very well with our result of 1.070(15) GeV 3/2 [8] based on using the HISQ action for a heavy quark combined with a HISQ c quark and working at a range of heavy quark masses between c and b on lattices with a range of lattice spacings from 0.15 fm down to 0.045 fm. The HISQ-HISQ result has an uncertainty which is a factor of 2 smaller than the NRQCD-HISQ result we give here. This is because the HISQ-HISQ current is absolutely normalised in the calculation of pseudoscalar decay constants and the calculation is done over a wider range of values of the lattice spacing for better control of discretisation errors. Good agreement between the HISQ-HISQ result and NRQCD-HISQ result TABLE XII. Column 2 gives the Coefficient zV i ,c needed for the one-loop renormalisation factor for the vector decay constant B * c . zV i ,c is constructed from results given for the appropriate NRQCD bare masses and massive HISQ quarks with the appropriate values of mca in [31] as zc = η0 −τ10. Column 3 gives δzc, which is the difference between zV i ,c and zA 0 ,c from Table XI. Column 4 gives the difference between δzc and the corresponding value δz for massless HISQ quarks (from Table VII). The values of αs used with these z coefficients are given in Table I. Column 5 gives the unrenormalised ratio of vector to pseudoscalar amplitudes (see text) determined from the simultaneous fit to B * c and Bc meson correlators. Results from columns 4 and 5 are used in the determination of Rc/Rs.

Set
zV i ,c δzc δzc − δz R unren.  (17)  7 -0.013(5) 0.021 (7) 0.058 (8) 0.99766 (22) was already seen for the B s in [5,38] and this further test increases our confidence in our handling of lattice QCD errors. In particular it is an important test of our normalisation of improved NRQCD-HISQ currents that are also in use for semileptonic decay rate calculations underway on these gluon field configurations.
Using the experimental value of the B c meson mass of 6.276(1) GeV [39] we can convert our value for f √ M into a result for the decay constant: f Bc = 0.434 (15)GeV. (20) Again this agrees with our earlier result using HISQ quarks of 0.427(6) GeV [8].
The vector to pseudoscalar decay constant ratio, R c , is obtained in an analogous way to that for the B s and B * s mesons described in Section III A. The formula we use is that given in eq. 10, in which z A0 and z Vi are the coefficients calculated for temporal axial and spatial vector currents respectively using HISQ quark mass values appropriate to c quarks. Values for z A0,c are given in Table XI and values for z Vi,c are given in Table XII. Table XII also gives δz c , the difference between the two, which is needed for the decay constant ratio in eq. 10. Note that we are now neglecting radiative corrections to the current correction J (1) (i.e. the terms with coefficients z 1 and z 2 in eq. (18)) since we do not have these terms for the vector current. Table X gives the results for the amplitudes that we need to construct the decay constants and their ratio. We see that the qualitative features of the results are the same i.e. that the amplitude of the leading order current is smaller for the vector than for the pseudoscalar meson, lowering the vector to pseudoscalar decay constant ratio, whereas the current correction contributions have opposite effect. The impact of the current corrections is a few percent less than in the B s case but varies more strongly with lattice spacing.
In a similar approach to that used for R l in Section III A we will study R c through its ratio with R s . In this case the renormalisation factor does not cancel completely at O(α s ) since δz c is not equal to δz. Instead we have a renormalisation factor for R c /R s which is (1 + [δz c − δz]α s ), i.e. we can write: Values of (δz c − δz) are given in Table XII. We see that these are small and independent of the value of am b within statistical uncertainties. Since the dependence on am b comes from the NRQCD action it is not surprising to find some cancellation between these two cases. The remaining small renormalisation then reflects the fact that the c quark mass is not zero (i.e. m c /m b = 0). Table XII gives in the final column results for the appropriate ratio of sums of amplitudes needed in eq. (21), i.e. R unren.
c . This can be combined with R unren. s from Table V and the small renormalisation applied to form R c /R s . Figure 4 gives results for the ratio of R c to R s from eq. (21) as a function of lattice spacing. We see that the results are independent of lattice spacing and sea quark mass. Importantly the values obtained are all significantly larger than 1.0, showing that the vector to pseudoscalar decay constant ratio is sensitive to the mass of the light quark combined with the b. Half of the difference from 1.0 comes from the raw amplitudes and the other half from the renormalisation factor in eq. (21).
In fitting this ratio as a function of lattice spacing to obtain a physical result we will use the same fit form as that used in subsection III A, eq. (13). The only change in form that we make is to remove am b dependence from the coefficient of unknown α 2 s renormalisation terms in the factor F 2 assuming they follow the same form as discussed above for the O(α s ) term. We take the radiative correction terms for the J (1) currents in F 2 to have the same form as in eq. (13) but allow for cancellation between R c and R s by giving that term coefficient 0.03 rather than 0.08. For the discretisation errors and m b tuning error terms in F 1 we likewise allow for cancellation between R c and R s by giving these terms coefficient 0.1 rather than 0.2.
The fit to our results gives χ 2 /dof of 0.1 and a physical result of R c R s = 1.037(12) where we have allowed a 0.5% uncertainty from missing higher-order relativistic current corrections. This value is 3σ greater than 1.0 giving a clear indication that R q increases as the mass of the quark q increases. This is consistent with what was found (with much lower significance) in Section III A for q = l and s. Our physical result for R c /R s is plotted as the grey band in Figure 4. A full error budget for R c /R s is given in Table IX. Using our earlier value for R s of 0.957(23) we obtain R c = 0.992 (27). We convert R c into a ratio of the decay constants of the B * c and B c mesons by dividing by the square root of the ratio of the masses. For this we use the experimental value for the B c mass of 6.276(1) GeV [39] and our lattice QCD result for the mass difference between B * c and B c of 54(3) MeV [10]. This gives a mass ratio for B * c to B c of 1.0086 (5). We then obtain f B * c f Bc = 0.988 (27).

E. Heavy Quark Mass Dependence
In this subsection we give results for calculations that use NRQCD quarks with masses lighter than that of the b in order to study the heavy-quark mass dependence of decay constants and their ratios and make a link between b and c. Using the HISQ action we have previously mapped out the dependence of pseudoscalar decay constants and quark masses [8,17,33,38] in this region in some detail, and we will be able to compare to these results.
The HISQ action has the smallest discretisation errors of any quark action in current use, since it removes treelevel a 2 errors and has no odd powers of a appearing. It is therefore a very good action for c physics [7,26,41]. Raising the mass from that of c requires fine lattices to keep masses in lattice units below ma = 1, where, naively, it might be expected that discretisation errors would become large. It is possible to reach the b on 'ultrafine' lattices with a lattice spacing as small as a = 0.045 fm. This TABLE XIII. The coefficients c1, c5, c4 used in the NRQCD action (eq. (2)) for values of the heavy quark mass in lattice units given in column 2. c6 is equal to c1 and c2 and c3 are set to 1.0. Amplitudes for J (0) and J (1) for temporal axial and vector currents between the vacuum and heavy pseudoscalar and vector mesons respectively, made from a heavy NRQCD quark and a HISQ s quark and denoted Hs and H * s . The different heavy quark masses in lattice units used on sets 1 and 4 are given in column 2. Results are in lattice units and the errors given are statistical/fit errors only.  Tables V and XIII). The open blue squares are results from a = 0.044 fm lattices using the HISQ formalism for both b and s [38]. The solid error bars on both sets of points include statistics and the (correlated) uncertainty in the value of the lattice spacing. The dotted error bars on the NRQCD points include in addition an estimate of NRQCD systematic errors [5]. The black bursts are the final physical values for the Ds and Bs [5,26].
has given accurate results for m b , f Bs and f Bc [8,17,38] because we can use operators that are absolutely normalised.
For NRQCD the issues are complementary ones. In this case we have systematic control of a non-relativistic effective theory. Discretisation errors are much smaller, having a scale set by internal momenta rather than the quark mass. In this case naive arguments suggest that we need ma > 1 to control coefficients of relativistic correction operators, for high precision. In fact for b quarks on the ensembles we use here, with lattice spacing values ranging from 0.15 fm down to 0.09 fm, values of ma are well above 1 and there is significant headroom to reduce the mass, particularly on the coarser lattices. Since the ratio of c to b quark mass is 4.5 [17], we cannot reach the c quark mass with ma > 1 even on the very coarse lattices. However, it is still of interest to vary the mass and com-  (Tables VIII and XIII). Dotted error bars include an estimate of NRQCD systematic errors. Black bursts indicate the physical result for B * s /Bs mesons from this paper and for D * s /Ds mesons using HISQ c and s quarks from [40].
pare the mass-dependence using NRQCD heavy quarks to that obtained from a completely different perspective, in terms of systematic errors, using HISQ quarks.
We have already shown that using HISQ b quarks and NRQCD b quarks gives results in agreement for the decay constant of the B s [5,38] and the B c ( [8] and subsection III D). Here we will illustrate how well this agreement continues to lighter masses.
We work on one ensemble each from the very coarse (set 1) and coarse (set 4) lattices. We will focus on results using s HISQ quarks where, as we have seen, dependence on the sea u/d mass is negligible. It is most convenient to use the same values of am b as those used before on the finer lattices, since then the determination of the radiatively-improved coefficients in the NRQCD action is simple. In Table XIII we give the coefficients that we use for heavy quark masses am h = 1.91 and 2.66 on very coarse set 1 and for am h = 1.91 and 3.297 on coarse set 4. On very coarse set 1 the lightest am h then corresponds to 1.91/3.297= 0.58 times m b . On coarse set 4 the mass 3.297 is higher than m b (since there am b = 2.66, see Table III), but am = 1.91 corresponds to 0.72 times m b . The coefficients are calculated by combining the one-loop coefficients at the appropriate am b values given in [6,20] with the appropriate α s value (also given in [6]) for that lattice spacing. These coefficients are then used in the NRQCD action (eq. (2)) along with relevant tadpoleimprovement factors given in Table III for that ensemble. We again use a local and two smeared sources for the NRQCD propagators, with smearing radii as given in Table III. We combine the NRQCD propagators with those for the HISQ s quarks on each ensemble. Table XIII gives results for the amplitudes for the leading-order and relativistic correction currents for the heavy-strange pseudoscalar meson (H s ) and vector meson (H * s ). These are obtained from simultaneous fits to the vector and pseudoscalar meson correlators as described in subsection II C.
To determine the pseudoscalar decay constant, f Hs , we use a more accurate formula than the one given in eq. (9). Instead we use the formula accurate through α s Λ/m h given in [5]: A0 . z 1 α s is an additional radiative correction to the subleading current J (1) . z 2 α s multiplies an additional subleading current which has the same matrix element as J (1) and so does not need to be separately calculated. The coefficients z 1 and z 2 are given for the masses we use in Table XIII. The z A0 values are in Table VII and α s values in Table I. Figure 5 shows results for f Hs M Hs as a function of inverse heavy quark mass in units of the physical b quark mass. The results for set 1 are shown as open red circles including the value at the b quark mass (am b = 3.297) from Table V as well as the results for lighter heavy  quark masses from Table XIII. The solid error bar is the dominant error in the raw results coming from the uncertainty in the lattice spacing. The dotted error bar includes an estimate of systematic errors from NRQCD coming from missing α 2 s renormalisation and (Λ/m h ) 2 current corrections. The latter systematic error grows as m h falls. The open blue squares give results from 'ultrafine' (a =0.044fm) n f = 2 + 1 lattices using the HISQ formalism for the heavy quark [38]. These results were part of an analysis of the heavy-strange pseudoscalar meson decay constant that spanned the range from c to b.
The plot shows good consistency between the two sets of results, which use very different formalisms on lattices that differ in lattice spacing by over a factor of 3. The black stars mark the final physical result for the B s [5] and D s [26] decay constants obtained by HPQCD after performing a fit including discretisation uncertainties.
Table XIII also includes results for the vector to pseudoscalar ratio of decay constants, R s = f H * s M H * s /f Hs M Hs . This is defined from eq. (10) up to missing α 2 s and α s Λ/m b matching uncertainties. These are plotted as a function of the inverse heavy quark mass in units of the b quark mass in Figure 6, including also results from Table VIII at the b quark mass. We see, as expected, that the values rise as m b /m h grows towards the c quark mass. The dotted error bars include an estimate of the (correlated) systematic error from missing factors in the matching of the NRQCD current to full QCD. These are estimated by rescaling results from our study here for the B s (subsection III A). The missing terms are: α 2 s terms in the overall renormalisation which are taken to be independent of m h ; α s Λ/m h current corrections which grow linearly with m b /m h and (Λ/m h ) 2 current corrections which grow quadratically.
The black bursts mark the physical result at the B s obtained in subsection III A and the result at the D s obtained using HISQ c and s quarks in [40]. The mass dependence of our NRQCD results is consistent with a value for R s that grows from our result at the B s towards the result we obtained at the D s with a relativistic formalism. The growth of the NRQCD systematic errors and indeed the fact that the c quark in a D s is not very nonrelativistic mean that we cannot accurately extrapolate from results here around the B s to the D s . We can estimate the slope at the B s , however. Our results on the coarse lattices, set 4, give a linear slope for the ratio R s of 0.050(17)m b /m h close to the b, where the uncertainty comes from NRQCD systematic errors in the current matching.

IV. DISCUSSION
Naively we expect heavy mesons with the same valence quark content but with vector or pseudoscalar quantum numbers to be very similar since spin-dependent (hyperfine) intereactions that distinguish between quark and antiquark having spins parallel or anti-parallel are suppressed by the quark mass. Such arguments in respect of the meson masses are straightforward to make even within the quark model. For NRQCD the dominant source of such effects for the vector to pseudoscalar meson mass difference is the term proportional to c 4 in the NRQCD Hamiltonian, eq. (2) [10,22].
For decay constants the arguments are more subtle which is why lattice QCD calculations are important to pin down the results. Viewed from the perspective of a nonrelativistic effective theory, there are three sources for terms that affect the ratio of vector to pseudoscalar decay constants for heavy mesons: one is the hyperfine term in the Hamiltonian as above, the second is the relativistic current correction terms (J (1) in eq. 5) and the third is matching of the current operator to full QCD. The first two give an effect that is proportional to Λ/m h whereas the third gives corrections to 1 that are proportional to α s . The different dependence of the three effects and the possibilities of cancellation between them have  [3] to interpolate beween c and the infinite mass limit. The lowest two sets of results in purple and orange use QCD sum rules [42,43]. The red dashed line marks the value 1.0. The ratio of f √ M for vector and pseudoscalar mesons with at least one heavy quark, plotted against the inverse of the mass for the pseudoscalar meson. Filled blue squares give the results obtained here for the B * s /Bs along with the results from [40] for fD * s /fD s multiplied by the square root of the ratio of the meson masses from experiment [2]. Filled green circles give the results for heavy-heavy mesons using experimental values for the vector decay constant obtained from the leptonic width [2] and results from lattice QCD using HISQ quarks for the pseudoscalar [8,26]. The filled red triangle gives the result from this paper for the B * c /Bc ratio.
given rise to a variety of predictions for the ratio of decay constants for vector and pseudoscalar heavy-light mesons over the years and controversy has surrounded the question of whether the ratio is larger or smaller than 1 at the b quark mass. Our results here show that the ratio is less than 1 (to 2σ) for B * /B and B * s /B s mesons. In this subsection we set this in the context of earlier results.
A baseline that can be used for heavy-light mesons [36] is Heavy Quark Effective Theory (HQET) in which the quark Lagrangian becomes simple, with no spin dependence, in the infinite quark mass limit. The matrix elements of the spatial vector and temporal axial currents between the vacuum and heavy-light mesons become the same in this limit within the effective theory but the renormalisation factors that match the currents to full QCD are not the same. These have been calculated through O(α 2 s ) in [44], giving, to this leading nonrelativistic order and in terms of the M S coupling: with n f is the number of light quarks in the sea and the factor 0.18 takes account of the non-zero mass for the c quark if this is included in the sea. Evaluating the expression in eq. (25) at scale m b =4.18 gives 0.92, well below 1.0. The size of the two-loop coefficient here indicates that the series might better be expressed at a different scale [44] for evaluation. In the Brodsky-Lepage-Mackenzie (BLM) scheme [45], for example, which uses quark loops (i.e. the terms proportional to n f here) to gauge the size of a 'typical momentum' flowing through the diagram, we would obtain a scale for α s of 0.2m b and a series This gives a much better looking series at the cost of a low-scale for α s that means that the series cannot be evaluated for quarks with masses lower than m b . Evaluating eq. (27) for 0.2m b ≈ 1GeV gives a result of 0.85, 15% below 1. Early calculations added sum-rule estimates of Λ/m h hyperfine and current corrections to the one-loop piece of eq. (25) and obtained a variety of results depending on the relative sign of hyperfine and current correction terms. In [35] it was found that the hyperfine and current corrections terms have opposite sign (in agreement with a subsequent lattice NRQCD study [34]) and this gave a vector to pseudoscalar decay constant ratio for b-light mesons of 1.00 (4). The central value would be reduced below 1.0 using the two-loop expression above.
The calculation we give here improves on this approach since it is a fully integrated calculation in lattice QCD.
We use an improved NRQCD action for the b quark accurate (for heavy-light calculations) through O(α s Λ/m h ) from which we can calculate the matrix elements of current operators nonperturbatively. The nonrelativistic current, combining the leading term and first, Λ/m b , relativistic correction is matched to full QCD and the O(α s ) matching correction is found to be very small.
Our results, as described in Subsection III A, show that f B * s /f Bs and f B * /f B are about 5.0(2.5)% below 1, and the ratio for the B * /B is 1.3(1.4)% below that of the B * s /B s . Our values are compared to results from two recent QCD sum-rule analyses [42,43] in Figure 7. Although there is some tension, those results are consistent with each other and with our numbers here. All the results show the same tendency for the ratio for B * /B to be slightly smaller than for B * s /B s , although the difference is not significant in any of the cases.
We also compare to a recent lattice QCD result [3] which used the twisted-mass formalism for both heavy and light quarks on gluon field configurations that included the effect of u/d quarks (only) in the sea. The twisted-mass value is obtained from results calculated for heavy quark masses around the c quark mass and above. An interpolation between those results and the infinite mass limit is performed to reach the b, using the twoloop formula of eq. (25) to rescale results so that 1.0 is obtained in the infinite mass limit. The value quoted for f B * /f B is 1.051 (17) and this disagrees with our value by more than 3 (combined) standard deviations. It is not clear that results using only u/d quarks in the sea will necessarily agree with those, like ours, that include a full complement of sea quarks. This may be a case where the 'quenching' of the s quark produces a visible effect.
It is also interesting to compare results for the ratio of vector to pseudoscalar decay constants between heavyheavy mesons and heavy-light mesons. The decay constant of vector heavyonium mesons can be determined from their experimental decay rate to leptons: The decay constants can also be calculated in lattice QCD [41,46,47] and good agreement with experiment is found. Since heavyonium pseudoscalar mesons do not annihilate to a single particle, there is no direct experimental determination of the decay constant. Again, however, the decay constants can be accurately determined in lattice QCD [8,26]. Figure 8 shows the ratio of vector to pseudoscalar decay constants (multiplied by the square root of the ratio of the masses) for (J/ψ)/η c , Υ/η b , B * s /B s and D * s /D s plotted against the inverse of the corresponding pseudoscalar meson mass. For the J/ψ and Υ decay constants we use the values determined from the experimental annihilation rates [2] and eq. (28). These are 0.407(5) GeV and 0.689(5) GeV respectively. From full lattice QCD the η c decay constant is 0.3947(24) GeV [26] and the η b decay constant is 0.667(6) GeV [8]. The D * s /D s decay constant ratio is taken from [40]. We see that the behaviour for heavyonium and heavy-light mesons is similar but the slope is larger for heavy-light mesons. For heavyonium mesons, similar considerations apply to the decay constant ratio as discussed above for heavylight mesons. A baseline might be considered a simple spin-independent potential model in which the decay constant can be related to the 'wavefunction-at-the-origin'. However there are significant QCD radiative corrections to match ψ(0) to the decay constant in both the vector (see, for example [48]) and pseudoscalar [49] cases, and these need to be included. Going beyond this requires the inclusion of spin-dependent terms in the Hamiltonian and relativistic corrections to the leading-order current. These are taken care of in a lattice QCD calculation, either explicitly when using a nonrelativistic formalism such as NRQCD [47] or implicitly included when using a relativistic formalism such as HISQ [8].
Here we have calculated the decay constant of both the B c and the B * c , using NRQCD b quarks and HISQ c quarks and working through first-order in the QCD matching and relativistic spin-dependent corrections to the NRQCD Hamiltonian and the currents. Our result for the B c decay constant agrees well with that obtained previously using the relativistic HISQ formalism for both b and c quarks [8], adding confidence to our analysis of systematic errors in both the nonrelativistic and relativistic approach. Here we also calculate the ratio of decay constants for the B * c and B c , for the first time from lattice QCD.
B c and B * c decay constants have also been calculated within a potential-model approach, including QCD radiative corrections. See [50] for a discussion. Results are in reasonable agreement with ours, but with a larger uncertainty because the approach has less control of systematic errors.
We find a value for f B * c /f Bc which is larger than that of f B * s /f Bs , indicating that the internal structure of the B c is somewhat different from that of a typical heavy-light meson. Figure 8 shows this clearly. When the decay constant ratio is plotted for the B * c /B c it lies very neatly between the heavy-heavy line and the heavy-light line.

V. CONCLUSIONS
Decay constants, which parameterise the amplitude for a meson to annihilate to a single particle, are as much a part of a meson's 'fingerprint' as its mass. They are often harder to determine, however, and some cannot be accessed directly through an experimental decay rate. The overall picture of meson decay constants gives information about how the internal structure of mesons changes for different quark configurations as a result of QCD interactions. To obtain this picture in sufficient detail, for example even to put the decay constants into an order, requires calculations in full lattice QCD, since only then can we reliably quantify the systematic errors.  The lattice result for f π + is marked with a cross to indicate that it is used to set the scale in some analyses (although not here). The result for the K + is from [28], the B + and Bs from [5], the D + and Ds from [16], the φ from [9], the D * s from [40], the ηc from [26], the J/ψ from [41], the Bc and η b from [8], the Υ and Υ from [47] and the B * , B * s and B * c from this paper.
Here we have expanded range of decay constant calculations from full lattice QCD to include vector heavylight mesons. Our results for the ratio of vector to pseudoscalar decay constants are: f B * s f Bs = 0.953 (23) f B * c f Bc = 0.988 (27).

Thus
• The vector decay constant is smaller than the pseudoscalar decay constant for b-light mesons, at the 2σ level for B * /B and B * s /B s . This is in contrast to results for c-light mesons where the vector has a larger decay constant than the pseudoscalar.
• The ratio of vector to pseudoscalar decay constants shows an ordering so that f B * c /f Bc > f B * s /f Bs > f B * /f B . When correlations between the uncertainties are taken into account using ratios, the first of these relationships has 3σ significance, the second 1σ (see eqs. (22) and (16)).
Using our earlier world's best results for f B (0.186(4) GeV, isospin-averaged), f Bs (0.224(5) GeV) [5] and f Bc (0.427(6) GeV) [8] we derive values for the vector decay constants: f B * = 0.175(6) GeV (30) f B * s = 0.213(7) GeV f B * c = 0.422(13) GeV. Finally, in Figure 9 we give a 'spectrum' plot for the decay constants of 15 gold-plated mesons from lattice QCD, including the new results from this paper. It illustrates the coverage and predictive power of lattice QCD calculations. The decay constants are ordered by value, something that is only possible with sufficiently accurate results. The range of values is much smaller than that for meson masses and the ordering of values is not as obvious because the quark masses do not have the same impact on the decay constants as they do on the meson masses. The plot therefore shows up some interesting features in the ordering, for example that the K and B * mesons have such similar values and that the φ meson appears so far up the list. We see that the decay constants for vectorpseudoscalar pairs are close together everywhere, closer than for the pairings in which an s quark is substituted for a light quark in a meson, for example.
Future work will improve the accuracy of lattice QCD results for the vector-onium states such as the φ (not strictly gold-plated) [33] and the ψ [51], both of which can be determined accurately from experiment. The issues there are mainly from lattice QCD statistical errors.
For b-light meson decay constants the dominant source of uncertainty, as we have seen, is from systematic errors in NRQCD and work is underway to reduce these further.