Precise B, B_s and B_c meson spectroscopy from full lattice QCD

We give the first accurate results for $B$ and $B_s$ meson masses from lattice QCD including the effect of $u$, $d$ and $s$ sea quarks, and we improve an earlier value for the $B_c$ meson mass. By using the Highly Improved Staggered Quark action for $u/d$, $s$ and $c$ quarks and NRQCD for the $b$ quarks, we are able to achieve an accuracy in the masses of around 10 MeV. Our results are: $m_B$ = 5.291(18) GeV, $m_{B_s}$ = 5.363(11) GeV and $m_{B_c}$ = 6.280(10) GeV. Note that all QCD parameters here are tuned from other calculations, so these are parameter free tests of QCD against experiment. We also give scalar, $B_{s0}^*$, and axial vector, $B_{s1}$, meson masses. We find these to be slightly below threshold for decay to $BK$ and $B^*K$ respectively.


I. INTRODUCTION
B meson physics is one of the critical elements of the flavor physics programme. The mass differences between the 'heavy' and 'light' eigenstates of the neutral B and B s mesons are now known experimentally and can be used to precisely constrain the ratio of CKM elements |V td |/|V ts | if the appropriate theory results have been calculated with a matching error. A first lattice QCD calculation of these mixing matrix elements, including the effect of u, d and s sea quarks, was given recently and the critical quantity ξ = f Bs B Bs /f B √ B B was obtained to 3% [1]. Similarly experimental results for the annihilation of charged B mesons to leptons via a W boson can be used to constrain V ub and B semileptonic decays to π or D ( * ) can be used to constrain V ub and V cb if the appropriate theory calculations of decay constants or form factors are known. Again calculations of these in full lattice QCD have been done and, for example, the decay constant of the B meson is obtained to 7% [1].
To improve on these lattice QCD results requires pinning down and eliminating sources of systematic error and testing as stringently as possible that this has been done. Here we provide a calculation of the masses of B mesons along with their decay constants using b quarks in the NRQCD formalism [2] with light quarks in the Highly Improved Staggered Quark (HISQ) formalism [3]. The HISQ formalism has improved discretisation errors compared to the asqtad improved staggered quark formalism used in our previous calculations. Because we use NRQCD which can handle hadrons with either sin- * c.davies@physics.gla.ac.uk gle or multiple b quarks we are able to do an accurate calculation of the B meson masses by linking them to meson masses in bottomonium. This provides a strong test of systematic errors. We are able to handle all of the 4 lightest quarks -u/d, s and c -using the HISQ formalism and are therefore also able to calculate mass differences and decay constant ratios accurately between B, B s and B c mesons.
Section II outlines how the Lattice QCD calculation was done and Section III describes the analysis and results. Section IV compares the results to experiment and to other lattice QCD calculations. Section V gives our conclusion.

II. LATTICE QCD CALCULATION
We are concerned here with mesons with one valence b quark and a lighter valence anti-quark, either c, s, or l. We use the notation l, "light", to refer to either the u or d quark. Everywhere in these calculations m l = m u = m d , but we will correct for the effects of this, along with the effects of missing electromagnetic interactions, when we compare to experiment.
The bottom quark moves sufficiently slowly inside bound states that it is well described by a non-relativistic formulation (NRQCD) on lattices of moderate lattice spacing. The lighter partner quark, with balancing momentum in a meson with zero total momentum, is moving much faster and requires a relativistic formulation. For this we use the HISQ formulation, which offers control of discretization errors to the level where even the c quark can be treated relativistically.
Set β a (fm) au0P m asq l au0P m asq s L/T N conf × Nt 1 6.572 0.1583 (13)   and m asq s . The lattice spacing values in fm are determined using the ηs meson mass and decay constant [9] and given in column 3. Column 7 gives the number of configurations and time sources per configuration that we used for calculating correlators. On set 5 only half the number were used for light quarks.

A. The gluon configurations
We make use of the MILC collaborations's library of 2 + 1-flavor gauge configurations [4]. These have two degenerate flavors of light sea quarks and one flavor of strange sea quark, formulated with the Asqtad action [5][6][7]. The gluon action is Symanzik-improved through O(α s a 2 ) except for terms of O(n f α s a 2 ) where n f is the number of sea quarks. In fact these terms [8] are of similar size to the other α s a 2 terms so in practice the gluon action has α s a 2 discretisation errors. For this work we use five different ensembles at three different lattice spacings, with a ≈ 0.15, 0.11, and 0.09 fm. We refer to these as "very coarse", "coarse" and "fine", respectively. The configurations have large spatial volumes (> (2.4 fm) 3 ). Table I lists the specific ensembles used in this work.
In this work we use values of the lattice spacing, a, on each ensemble determined using the mass and decay constant of the pseudoscalar ss meson, the η s . Although this particle is not seen in the real world because of mixing with uu and dd which can be prevented on the lattice, its properties can be determined from those known from experiment of the π and K mesons, as described in [9]. Table I lists the values obtained in [9] for the ensembles we are using here. The values of a are larger on coarse lattices than those from the more traditional way of setting the lattice spacing using the parameter r 1 , but the results agree, as they should, in the continuum limit [9].
The HISQ action further reduces the residual O(α s a 2 ) discretization errors coming from taste-changing effects found in the Asqtad formulation. It does this with an additional fattening step applied to the gluon field coupling to the quarks [3,10]. The errors are reduced by about a factor of 3, making HISQ therefore a better ac-tion to use for l and s quarks.
We have shown that the HISQ action can even be used for c quarks [3], but in that case an additional step is needed. The 'ordinary' tree-level O(a 2 ) discretization errors coming from the finite difference discretization of the covariant derivative are eliminated in both HISQ and Asqtad formulations using an additional 3-link 'Naik' term. The Naik term corrects errors that would otherwise appear at O(pa) 2 in the quark, and therefore meson, dispersion relation of energy vs momentum. A nonperturbative value for the Naik term coefficient (written as 1 + ) can be derived by studying the dispersion relation of the η c meson and tuning the coefficient until the square of the speed of light in this relation is 1. Here we use the values of the Naik coefficient determined in this way in [10] appropriate to the values of m c a on each ensemble [37]. Tuning the Naik coefficient in this way removes all discretisation errors from the HISQ action at leading order in the square of the velocity of the c quark [3]. At subleading order in v c there will be O(α s a 2 ) errors.
On each configuration of the ensembles in Table I we have generated and stored random-wall source charm, strange and light valence propagators. These were calculated in the following way. On the source time-slice t 0 we generate a U (1) vector of random numbers η(t 0 ) x . Then we invert to get the HISQ propagator g HISQ (x, t 0 ): For the fine ensemble we use a different random-wall at each of four values of t 0 on every configuration. For the coarse and very coarse configurations we use two random-wall sources per configuration. Although the time sources are equally spaced, their position in time varies from configuration to configuration through the ensemble in a pseudo-random manner to further reduce auto-correlations within the ensemble.
We use the charm, strange and light HISQ quark propagators from [10], but note that here we are using a different definition of the lattice spacing and this affects the meson masses in physical units and therefore the tuning of the quark masses. We have included some additional quark mass values on the very coarse lattices to be able to correct for mistuning. We list the HISQ valence parameters used in columns 3 to 7 in Table II. The tuning of the HISQ valence c and s masses to their correct values on each ensemble is important to avoid mistuning effects masquerading as, for example, lattice spacing artefacts or sea quark mass dependence. We tune the c mass to give the correct η c meson mass on each ensemble. Here the correct η c mass has to be adjusted slightly from its experimental value of 2.980 GeV to allow for the fact that our lattice QCD calculation is happening in a world without electromagnetism and without c quarks in the sea and in which we do not allow the η c meson to annihilate to gluons (because we have not included such 'disconnected' pieces in our η c meson correlators). These effects are all small and can be estimated from potential models or perturbation theory. We  TABLE II: Parameters for the valence quarks. am b is the b quark mass in NRQCD, and u0L is the tadpole-improvement parameter used there [11]. We use stability parameter [11] n = 4 in NRQCD everywhere. Columns 4, 6 and 7 give the charm, strange and light bare quark masses for the HISQ action. 1 + is the coefficient of the Naik term in the charm case [3]. On set 1 we give additional values of m b and ms that were used for tuning purposes as described in the text. We also used alternative ci coefficients for the am b = 3.4 case, specifically c1 = c6 = 1.36 and c5 = 1.21, again as described in the text.
take the physical value of m ηc appropriate to our calculations to be 2.985(3) GeV [12], incorporating a shift of 2.6 MeV for electromagnetic effects and 2.4 MeV for annihilation effects. The effect of c quarks in the sea is negligible; this will be discussed in the next section. We tune the s quark mass from the value of the mass of the η s meson. This is a fictitious pseudoscalar meson made of s quarks, which is not allowed in our lattice calculation to mix, by annihilation, with uu and dd mesons. We cannot therefore fix its mass from experiment but must do within a lattice QCD calculation, in which we extrapolate results for π, K and η s meson masses and decay constants simultaneously to the continuum limit and physical point for the l and s quarks. This was done in [9] and the value M ηs = 0.6858(40) GeV was obtained, as the value appropriate to this lattice QCD world.
The light quark valence masses are taken to match approximately the light quark masses used in the sea. The way that this is done is to take the ratio of the light valence HISQ mass to the appropriate HISQ strange mass value to be the same as the ratio of the sea light quark mass to its appropriate s mass value (the u 0P factors cancel in that ratio). This can only be done approximately because the sea 'strange' quark masses are not very well tuned in some cases [4] and the correct value for s is not known very accurately. However, sea quark mass dependence is only a very small effect for everything calculated here, so this is not a big issue.

C. NRQCD b quarks
One can estimate the velocity of the quarks in bottomonium mesons by comparing radial excitation energies to masses. This shows that in the Υ the b is very non-relativistic with v 2 b ≈ 0.1 (in units of c 2 ). By comparison, in charmonium v 2 c ≈ 0.3. In a mixed system with a lighter quark the b quark is even slower than in bottomonium. Consider that in a B c meson the reduced quark mass is roughly 1.5 times that of the bb system and 4.5 times that of a cc system. For a constant mean kinetic energy [13] across all three systems, we then expect v 2 b ≈ 0.04 and v 2 c ≈ 0.35 in the B c . For a b quark combined with an even lighter quark, v 2 b will get even smaller. Assuming a mean momentum of Λ QCD for the b quark inside a B or B s would give v 2 b ≈ 0.01. An NRQCD approach is then well justified for the b quark in B c , B s and B l mesons when combined with a relativistic approach for the lighter quarks.
As we have used a random-wall source for the HISQ propagators, it is critical that we initialize the NRQCD b propagators b with the same random-wall function η(t 0 ) x as we used for the HISQ propagators. This is non-trivial because the HISQ propagators have one only Dirac component, whereas NRQCD propagators have explicit spin. The spin information for the HISQ propagators is tied up with the positions of source and sink, however the source site information is lost for a random wall source once the propagator has been calculated. If we had access to a HISQ propagator from all source points to all sink points we could reconstruct the full 4×4 spin structure by multiplying by Ω(x)Ω † (y) where Ω is the staggering operator, given as a product of Dirac gamma matrices as: Here we cannot apply the source Ω to the HISQ propagators, but instead we can apply it as a source for the NRQCD propagators. This then effectively 'undoes' the staggering transformation and gives a naive quark source that can be combined with a b quark source, but it is done after the staggered propagators have been made. The bquark source then has to have 4 spin components, rather than the usual two for NRQCD because we cannot separate the upper and lower components of the naive quark source after the fact. Our b quark source is constructed from Ω(x) multiplying the appropriate random noise at each site, η(t 0 ) x , updating the standard method of combining staggered quarks with other formalisms that have explicit spin components derived in [14]. A further issue is that we must project onto pseudoscalar or vector heavy-light mesons directly at the source of the b quark propagator rather than combining appropriate spin components at the end. We then have to calculate separate b quark propagators for each of the pseudoscalar and vector mesons.
Finally, to enhance our ability to isolate the groundstate behavior, we smear the b propagator source with a smearing function, S, of various functional forms and differing radii, r i . Combining all of these factors we therefore, on timeslice t 0 , initialize the NRQCD propagator as: where Γ is an element of the Dirac algebra chosen to project out either a pseudoscalar or a vector heavy-light meson state.
On subsequent timeslices we evolve the NRQCD propagator recursively in the standard way [15]. Note that here, because upper and lower components do not mix in NRQCD, the upper and lower halves of our b-quark source are simply evolved separately with the same evolution equation. Our b-quark propagators then have only 2 spin components at the sink end. The evolution equation is: with and E andB are improved versions of the naive lattice chromo-electric and chromo-magnetic fields, E and B.
All the gauge fields appearing are tadpole-improved by dividing by a tadpole factor, u 0L , derived from the mean link in Landau gauge. The equations given above represent the standard NRQCD action, used in many previous lattice QCD calculations (for example [15] and [11]), correct through v 4 b . The largest source of remaining systematic error from this action is from radiative corrections to the coefficients of the v 4 b terms required to match full QCD through O(α s v 4 b ) (the v 2 b term is tuned nonperturbatively in fixing the b quark mass as described below). Here we generally use the tree-level values of c i = 1 for the constants as we have done before. However, we have also done some calculations on set 1 using values of c 1 , c 5 and c 6 that include radiative corrections, to be able to gauge the size of the systematic error from these terms. How the radiative corrections are calculated and further tests of them will be described elsewhere [16,17]. The values we use on set 1 are those appropriate to a value of m b a of 3.4 and α s in the V scheme appropriate to the lattice spacing of that ensemble. These are c 1 = c 6 = 1.36 and c 5 = 1.21 [16].
We list the NRQCD valence b masses used in column 2 of Table II along with the u 0L parameters. We have used two different masses on set 1, again to test for systematics from quark mass tuning. Since NRQCD quarks propagate in one direction in time only we improve statistics by generating propagators both forwards in time (for T /2 time units) and backwards in time from each light quark source.
The b quark mass is tuned by determining Υ and η b meson masses [17]. Because the zero of energy has been shifted in NRQCD we cannot determine meson masses directly from their energy at zero momentum as is done in a relativistic formulation. Instead we must calculate the 'kinetic mass' of a meson, M , which appears in the relationship between E(p) and p 2 : We are able to do this very precisely using random wall sources [9]. We determine the kinetic mass of both the Υ and the η b mesons and use the spin-average of the two, i.e.
to tune the b quark mass. The reason that we use the spin-average is to avoid systematic errors from the terms in the NRQCD action that give rise to spin-splittings. These terms are only included at leading (v 4 b ) order in the action above (equation 6), and we know that there are sizeable discretisation errors in the hyperfine splitting between the Υ and η b as a result [11]. These errors also have a small effect on the kinetic masses [17] and we take the spin-average to remove them.
The experimental value for M bb that we tune to must be adjusted, as described above in the η c case, to be the value appropriate to our lattice QCD world which is missing some elements of the real world. We estimate the effect of the electromagnetic attraction of the b and b to be 1.6 MeV from a potential model. The absence of electromagnetism shifts both M Υ and M η b upwards. The shift from not allowing the η b to annihilate to gluons we take to be the same as for the η c at 2.4 MeV. In addition we must consider the absence of c quarks in the sea. This has only a very small effect on charmonium states [12] but a somewhat larger effect on bottomonium, because these states are more sensitive to exchange of higher momentum gluons that can generate a cc pair. The effect can be estimated perturbatively from that of a massive quark loop in the gluon propagator that gives rise to the heavy quark potential [12]. It gives rise to a correction to the potential which is proportional to a delta function at the origin: This correction is very similar to the hyperfine potential in bottomonium, except that it is not spin-dependent and the mass that appears in the denominator is m c and not m b . The hyperfine potential induces the mass difference between M Υ and M η b of 69 MeV [18]. The potential above has a coefficient 14 times smaller, where this factor is given by [80π/(3α s )] × (m c /m b ) 2 and we take m b /m c = 4.51 from [19]. We therefore expect a shift of around 5 MeV to both Υ and η b . Again the absence of c quarks in the sea shifts the mass values upwards. The experimental η b mass is now 9.391(3) GeV [20] and the Υ, 9.460 GeV, giving a spin-average of 9.443(1) GeV. Applying the shifts above with 50% errors, we find the appropriate value for M bb for us to tune to is 9.450(4)(1) GeV, where the first error is from applying the shifts and the second is the experimental error.

D. HISQ-NRQCD two-point functions
At the sink end of the propagators we must tie a twocomponent NRQCD propagator with the one-component HISQ propagator. We first convert the HISQ propagator back into a multi-spin object by multiplying by Ω(x) at each sink site: remembering that the Ω † factor that would normally be at the source end has now been included in the NRQCD propagators.
Finally we generate the B meson correlator by combining the HISQ propagator and NRQCD propagators at the sink timeslice t, again with smearing functions and appropriate Dirac structure: Note that this is a sum over 4 spin components at the source and 2 at the sink, so Γ is either a 2 × 2 unit matrix for the pseudoscalar meson or a Pauli spin matrix for a vector meson.

E. Extracting physical masses
We do simultaneous constrained fits [21] to the form where i and j index the 3 × 3 matrix of source and sink smearing functions. The second term fits the oscillating "parity-partner" states that appear in most staggered meson correlators. Our correlation functions cover the range of t − t 0 values from zero to T /2, although we do not fit all the way to t − t 0 = 0. Instead we start at t min = 2 − 4 for B s and B l fits and 6 − 8 for B c fits to reduce the effect of excited states. We constrain the parameters of the fit with prior values and widths, which are fed into the augmented χ 2 function that the fit minimizes. These priors represent very general information about mass splittings and amplitudes. The prior value for the ground state mass is simply taken from an effective mass plot, with the prior width taken to be a factor of 1.5 from this value. The mass parameters for the higher mass states enter the fit as the logarithm of the mass difference with the state immediately below, so that the mass difference is positive and higher mass states remain, by definition, higher mass. We take a prior value on these mass differences, both for normal states and oscillating states, to be ≈ 600 MeV, converted to lattice units for each fit. The prior width on the mass differences is taken to be a factor of 2. The prior on the mass of the lowest oscillating state is taken to be ≈ 400 MeV above the ground state, with a prior width of a factor of 1.5. We use 0.1 ±1.0 for the prior value and width for all amplitudes. Here the prior value of 0.1 is simply to provide a non-zero starting point for the fit. The width of 1.0 can be uniform across correlators of different smearing functions because we normalise them all so that x S(x) 2 = 1 across a timeslice.
In this way, we are able to obtain high-confidence fits which are stable, both in the central value for the ground state mass and amplitude and their errors, with respect to varying the number of exponential functions included, N exp . We take our results from fits with N exp = 5, since all our fits are stable by this point.
Where possible (sets 2-4) we simultaneously fit multiple light valence channels. That is, by fitting simultaneously C Bs and C B l , we can eliminate correlated errors from estimations of mass differences such as M Bs − M B l .
One important issue with B meson correlators is their exponentially falling signal/noise ratio, which means that the statistical accuracy that can be obtained on masses and decay constants is not as high as that of lighter mesons, for example D/D s [12]. The variance of the B s correlator, for example, contains bbss propagators and can rearrange them into an η b and an η s . Thus the noise (square root of the variance) falls exponentially with a lowest energy (E η b + m ηs )/2 at large times while the signal falls with a lowest energy E Bs . This means that the signal/noise degrades exponentially with a physical energy which is the mass difference between the B s and (M η b + M ηs )/2 (330 MeV). This is illustrated in Figure 1 where we explicitly compare the effective mass of the B s correlator and the effective mass of its statistical error, and show that the 'mass in the noise' is as expected. This physical mass difference cannot be altered but if we use smearing functions, as we have done here, it is possible to extract the ground state B s mass from early t − t 0 values, where the noise is less of an issue. [38] As discussed earlier, the zero of energy is changed in the NRQCD formulation so the energy parameters E k and E k include an energy shift for which we must correct before comparing to experimental values. We can do this by comparing the B meson state of interest (containing 1 b quark) to a reference state, which can also be calculated with NRQCD b quarks on the lattice and whose mass is known measured experimentally. That is, where E ref is calculated on the same lattice ensemble, and M ref comes from experiment (adjusted if necessary for the absence of electromagnetism etc. from our calculation). n is the number of b quarks in the reference state. The reference state can also be a linear combination of states, such as the spin-average of bottomonium states that we will use below.
To minimize the contribution of the 0.8% uncertainty on a −1 to the overall uncertainty in M B , it is important to choose a reference state that makes the quantity (E B − E ref /n) as small as possible. We will sometimes do this below by subtracting the masses of additional reference states, for example ones made purely of c quarks for the B c case.

A. Bs meson mass
To determine M Bs we follow the strategy described earlier, using the spin-average of bb states as a reference, and calculate: From this we can reconstruct M Bs using Here E bb is the spin-average of Υ and η b energies at zero momentum calculated with the same NRQCD action and on the same configurations as used for calculating the B s meson energies. M bb is used to tune the b quark mass, as discussed earlier, and M bb,phys is the physical value taken from experiment, but adjusted (to 9.450 GeV) for the lattice QCD world (missing electromagnetism, η b annihilation and charm-in-the-sea). To compare our results for M Bs,latt to experiment we have to add corrections to put back in missing electromagnetism and charm-in-thesea effects. These corrections are negligible, however, as we will discuss below. Table III lists all our fitted values needed for determination of the B s meson mass. Note that the error on the fitted B s meson energy is larger than any of the errors on the other fitted energies. This is because of the signal/noise problem in the B s correlator discussed earlier. The other correlators used here do not have that problem and the fits give much more precise results for ground state masses. Details of these other fits are given elsewhere [9].
In Figures 2 and 3 we show how ∆ Bs varies with the square of the η s mass and the spin-averaged bb mass from our results on set 1. These results allow us to correct for, and estimate the errors from, mistuning b and s masses. The lines are simple linear fits in M bb and M 2 ηs . The slope of ∆ Bs against M 2 ηs is 0.19, in good agreement with that expected from the experimental data comparing B s and B. The slope against M bb is very small because the b quark mass effects naively cancel in ∆ Bs . However some residual dependence remains and gives a slope of 0.017, somewhat smaller than the experimental result of 0.033 obtained over a much larger mass range from comparing B s and D s . Table IV gives the values of ∆ Bs , adjusted for mistuning by using the slopes given above and the mismatch of M bb and M 2 ηs on each ensemble compared to the physical values (9.450 GeV and (0.6858 GeV) 2 respectively). We take a 50% error on any shift applied for mistuning.  For each set we list the valence b quark mass and its associated kinetic mass and energy for the spin-average of Υ and η b states. We also give the valence s quark mass and its associated ηs meson mass. These values are also given in [9]. Where we have used multiple b and s masses on set 1, we give the ηs and E bb values only once to avoid confusion. In column 7, we give the fitted energy of the Bs meson (i.e. E0 from fits to the form given in equation 12). In column 8 we give the hyperfine splitting, ∆ hyp This column is largely from [22] but includes some additional values on set 1 that we use for studying systematic errors. Note also that the value for set 1 on line 1 is different from that in [22] although consistent with it. Here we use a value from a fit to the Bs and B * s correlators alone, rather than from a full fit including B l and B * l , to be in keeping with the other Bs values given on set 1. Columns 9 and 10 give the values of mass differences between scalar and pseudoscalar and between axial vector and vector respectively, discussed in section III E. Errors from mistuning are smaller than the statistical errors except on sets 1 and 2. Table IV also gives the error from the uncertainty in the lattice spacing from Table I. The error is smaller by a factor of 2 than the naive result of multiplying ∆ Bs by the percentage error in a −1 . The reason is that changing the lattice spacing requires the masses to be retuned and this affects ∆ in the opposite direction. The resulting errors on the tuned values for ∆ Bs are typically less than 1%, around 4 MeV. Within these uncertainties we are not able to distinguish any dependence on sea quark masses or the lattice spacing. Sea quark mass effects are expected to be very small, because the B s is a gold-plated particle and has no valence light quarks. The lattice spacing dependence depends on the quantity Set ∆B s (GeV) δx l δxs 1 0.6392 (16)(28)(25)  chosen to fix the lattice spacing. Earlier reporting of these results [23], using the variable r 1 to set the scale, did show visible lattice spacing dependence. Here it appears, perhaps not surprisingly, as if ∆ Bs has the same discretisation errors as the η s used to the fix the scale.
The tuned results from Table IV are used to reconstruct M Bs,latt (using equation 15) and this is plotted in Figure 4 against the square of the lattice spacing. In order to quote a physical value that can be compared to experiment we need to fit our results as a function of lattice spacing and sea quark mass so that systematic errors from such dependence can be fed into our final error. The sea quark mass dependence we take to be a simple polynomial form in the variables δx s and δx l , defined by δx q = (m q,sea − m q,sea,phys )/m s,sea,phys . These variables were used in [12] but must be adjusted here consistently for the change in definition of the lattice spacing and the new values are given in Table IV. Any sea quark mass dependence identified in our fit can be extrapolated to the physical point where δx l = δx s = 0, and our errors allow for dependence not resolved by our fit.
The lattice spacing dependence is a trickier issue in NRQCD because we cannot extrapolate naively to a = 0. What we need to do is to fit the lattice spacing dependence and assess, using information from the fit, how much of the dependence is physical and how much unphysical, and allow for both in the final error. Physical dependence on the lattice spacing will arise from discretisation errors in the gluon and sea quark actions, and in the light valence quark (HISQ) action. We expect this dependence to be O(α s a 2 ) at leading order as discussed earlier.
The NRQCD action also has discretisation errors. These are corrected at tree level by the terms with coefficients c 5 and c 6 in equation 6. Beyond tree level c 5 and c 6 have an expansion in powers of α s , required for NRQCD to match QCD at that order, whose coefficients depend on am b . This dependence will typically be mild for large am b but become singular as am b → 0. This has been explicitly checked for the c 5 and c 6 coefficients for a slightly different action in [24] and results have also been derived for this action [16,17]. They show almost no am b dependence for am b > 1. In general, however, the coefficients of discretisation corrections can be am b dependent in NRQCD and therefore our discretisation errors can be am b dependent. We need to allow for a mild nonsingular dependence (i.e. appropriate to the values of am b that we are using) in our fits, so that the systematic error from this can appear in our final results. Since any smooth function can be expanded over a limited range using a polynomial, we simply allow for linear and quadratic terms in the variable δx m = (m b a − 2.7)/1.45. The factors 2.7 and 1.45 are chosen so that δx m changes from -0.5 on the fine lattices to +0.5 on the very coarse lattices as m b a covers the range that we have used.
We therefore fit our results for ∆ Bs to the following form: We take the prior on ∆ Bs,phys to be 0.6(2). The priors on the sea quark mass dependence, b l and b s , are taken to be 0.00 (7). Sea quark mass effects are suppressed by a factor of 3 over valence mass effects and here valence effects correspond to a slope in quark mass of less than 0.2. The priors for the b parameters corresponding to the quadratic sea mass dependence are then set to (0.2) 2 /3 i.e. 0.000 (13). We take the scale of the physical a-dependence to be the scale of Λ = 400 MeV, since we expect it to be set by typical internal meson momenta in QCD. The coefficients of the quadratic a dependence, c 1 , c l and c s , should be O(α s ) so we take priors of 0.0(5). For c 2 , and the am b dependence of the discretisation errors, c jb and c jbb , we take a very conservative prior of 0(1). The result from the fit is ∆ Bs,phys = 0.638(6) GeV. The fit sees no dependence on lattice spacing, am b or sea quark masses but our final error allows for all of these. The fit, and its error, is robust under changes in the number of fit parameters, for example, including or not including the a 4 terms in equation 16. It is also robust under changes the prior widths. For example, doubling the prior width on the lattice spacing or sea quark mass dependence changes the final result by less than 1 MeV.
In the error budget in Table V we separate the 6 MeV error into component parts coming from the errors on the original data points (statistics, tuning and uncertainty in the lattice spacing) and the errors coming from the lattice spacing and sea quark mass dependence of the results, using the method described in [12]. The error on the original data points dominates.
To reconstruct M Bs,latt we must add 9.450/2 GeV to ∆ Bs,phys as in equation 15. This gives M Bs in the lattice world with no electromagnetism or c quarks in the sea. The latter effect should be negligible for the B s since it  Table V. The black star is the experimental result [20], offset from a 2 = 0 for clarity. is a much larger particle than the Υ or η b and therefore much less sensitive to the gluon exchange that could create a cc pair. The effect of electromagnetism can be estimated following [12]. There we gave a phenomenological formula for electromagnetic and m u /m d mass difference effects in heavy-light mesons: M sim is the mass of the meson in the absence of electromagnetism and with m u = m d . From experimental charged and neutral B and D meson masses we determined A ≈ 4 MeV and B ≈ 3 MeV. For the B s then this formula gives a shift between M (Q, q) and M sim (Q, q) as a result of electromagnetism of -0.1 MeV, a very small effect. We make no correction for this, but add an error for it to our error budget. Additional systematic errors that must be added in to the error budget are those from relativistic corrections that are not included in our NRQCD action. These errors affect results at all lattice spacings equally and so cannot be estimated from our results as we have done for the discretisation errors. We must consider the effect of relativistic corrections on both the B s mass and on the Υ and η b masses because they both appear in ∆ Bs . In fact we expect relativistic corrections to have a bigger effect in bottomonium than on the B s . Our NRQCD action is cor- Our NRQCD action already includes high order terms in Λ/m b at tree level and so there are no significant tree level errors for the B s . The leading error is at α s v b in missing radiative corrections to c 4 , the coefficient that multiplies the σ · B term giving rise to the hyperfine splitting. As we will discuss in section III D we have plenty of evidence that errors coming from this term are small, at most 10% of the hyperfine splitting itself. Since this error would vanish for the spin-average of the masses of the B s and the B * s , which is 3/4 of the hyperfine splitting above the B s mass, we take the error in the B s mass to be 3/4 of the error in the hyperfine splitting, 3.5 MeV.
Errors from the uncertainties in M ηs and M bb that we use to tune the s and b quark masses can be estimated from the slopes in Figures 2 and 3. The 4 MeV uncertainty in M ηs feeds into a 1 MeV uncertainty in ∆ Bs and therefore M Bs . The uncertainty in ∆ Bs from the 5 MeV uncertainty (simply adding the statistical and systematic errors) in M bb is very small -less than 0.1 MeV because of the cancellation of the b quark mass inside ∆ Bs . However, the uncertainty reappears when we reconstruct M Bs by adding M bb /2 to ∆ Bs . This then gives a sizeable 2.5 MeV error in M Bs .
Finite volume errors are expected from chiral perturbation theory to be negligible for the masses of mesons containing heavy quarks on volumes exceeding (2.4fm) 3 , that we are using here.
The full error budget is given in Table V. The systematic errors discussed above, added in quadrature, give 9 MeV, dominating the 6 MeV errors coming from the sta-tistical errors of the data and its lattice spacing and sea quark mass dependence. The final result is then M Bs = 5.363(6)(9) GeV. Figure 4 shows a dark shaded band for the first error and a lighter shaded band to encompass the full error, adding 6 MeV and 9 MeV in quadrature to give 11 MeV. To reduce the full error will require improvements to the NRQCD action, currently underway. The experimental result for the B s mass is 5.3663(6) GeV.

B. Bc mass
For the B c meson mass we could use exactly the same procedure as for the B s . However, there is a better method, in which we subtract in addition the mass of a charmonium reference state, the η c , to reduce the energy difference calculated on the lattice to a very small value: We call this the "heavy-heavy" (hh) subtraction method.
Here M ηc is the value of the η c mass calculated on the lattice and M ηc,phys is its value from experiment appropriately adjusted for the lattice QCD world, as described earlier.
For the HISQ quarks that we use for c the energy obtained from fits to charmonium correlators at zero momentum is the charmonium mass, so there is no issue with the zero of energy. We simply use the additional charmonium subtraction here to reduce errors from the uncertainty in the lattice spacing. We will compare results of this to a second method: that we call the "heavy-strange" (hs) subtraction method. Here we are using the B s meson to cancel the NRQCD shift of the zero of energy in the B c . Again the subtraction of the D s meson mass, calculated with HISQ c and s quarks is simply to reduce lattice spacing errors from the mass difference. We first discuss results from the hh method. The B c energies and M ηc masses are given in Table VI and the bb energies and masses, already used in the determination of the B s mass, are given in Table III. As before, we have to tune quark masses on each ensemble to their correct value. We show in Figure 5 how the splitting depends on M bb from our results on set 1 at two values of am b and two values of am c . We see that the slope is very small, 0.014, because very little b quark mass dependence is left after cancellation in this mass difference.   An estimate can be derived for the expected slope by comparing results for the b quark mass set to the value of the c quark mass (when ∆ Bc,hh becomes -3/8 times the charmonium hyperfine splitting). This gives a slope of 0.016 over a much wider range. Figure 5 also shows the value of the mass difference for the case where we use an NRQCD action with c 1 , c 5 and c 6 set to the values including O(α s ) radiative corrections appropriate for set 1. We see that this makes negligible difference. Figure 6 shows the results as a function of M ηc . The slope here is very small but in the opposite direction to that for the dependence on M bb . The value of the slope is -0.004. Based on the arguments above we would expect a slope of opposite sign but about 60% that of the b-quark mass dependence. However, as stated above, this estimate is made over a much larger range than the Figure. We can use the slope against M bb and against M ηc to correct for the slight mistunings of the b quark and the c quark that we have on some ensembles. Even though the shifts from M bb and M ηc dependence are very small they are not negligible. This is because ∆ Bc,hh itself is very small and also because it is very precise, since all of the energies involved have tiny statistical errors. We take a 50% error on b quark mistuning but a 200% error from c quark mistuning to allow for the fact that we may be underestimating the slope with c quark mass because of discretisation errors in the HISQ action for c on the very coarse lattices. Table VII gives tuned values for ∆ Bc,hh on each ensemble, along with three errors; that from statistics, from tuning and from the uncertainty in the lattice spacing. These latter two errors dominate. As before, variation in the value of the lattice spacing means that masses must be retuned. Here this has the effect of producing a net change equal to the naive lattice spacing error. Figure 7 plots these results against the square of the lattice spacing, after reconstructing the B c mass by adding back in (M bb,phys + M ηc,phys )/2 = 6.2175 GeV.
Clear lattice spacing dependence is visible in Figure 7 but there is no sign of sea quark mass dependence. As for the case of ∆ Bs we fit the results for ∆ Bc,hh as a function of lattice spacing and sea quark mass to extract a physical result. The fit form is essentially the same as for ∆ Bs . However, because ∆ Bc,hh is such a small quantity it cannot set the scale for the discretisation and sea quark mass effects. So, instead of allowing a function of a, δx l and δx s to multiply ∆ Bc,hh we add such a function with a multiplicative factor of 0.4 GeV, representing a typical scale for QCD binding energies. We include more terms for discretisation errors than in the B s case and set their scale by m c , rather than Λ, because in this case discretisation errors will come largely from the HISQ action for We take the prior on ∆ Bc,hh,phys to be 0.05 (5). The priors on the sea quark mass dependence, b l and b s , are taken to be 0.00 (7), and on the parameters corresponding to the quadratic sea mass dependence 0.000(13), as before. We take the scale of the physical a-dependence to be the scale of m c ≈ 1 GeV, as discussed above. The coefficients of the quadratic a dependence, c 1 , c l and c s , should be O(α s ) so we take priors of 0.0(5). For c 2 , and the am b dependence of the discretisation errors, c jb and c jbb , we take a very conservative prior of 0(1).
The result for ∆ Bc,hh,phys is 0.0616(42) GeV, giving a B c mass of 6.279(4) GeV. The fit is robust under changes of the prior values. For example, we tried the following changes: • taking the multiplier of a-and m sea -dependence to be 0.8 instead of 0.4; • taking the priors for a-dependence to have width 2 instead of 1;    Table III as are the ηs meson masses for each s quark mass and the corresponding Bs meson energies. In column 4 we give the ηc meson mass corresponding to each value of amc, and in the final column we give the corresponding Ds meson mass. Note that the values for Mη c and MD s are different from those reported in [12] because here we are using a nonperturbatively determined Naik coefficient as discussed in the text. In column 5, we give the fitted energy of the Bc meson (i.e. E0 from fits to the form given in equation 12). In column 6 we give the hyperfine splitting, ∆ hyp c = E(B * c ) − E(Bc), discussed in section III D. This column is largely from [22] but includes some additional values on set 1 that we use for studying systematic errors.
Set ∆ Bc,hh (GeV) -∆ Bc,hs (GeV) 1 0.0980 (4)(12)(8) 1.034(2)(16)(4)   VII: Results for ∆B c , the mass difference between the Bc meson and a particular reference mass, on different ensembles after tuning to the correct valence b, c and, where appropriate, s quark masses. The 3 errors listed are statistical, tuning and from the uncertainty in the lattice spacing. Column 2 gives results from the hh method and column 3 from the hs method, as described in the text.
• taking the priors on sea quark mass dependence to be 0.5 rather than 0.2.
None of these changed the result by more than 1 MeV. Our result is for a world without electromagnetism or charm quarks in the sea. The effects of electromagnetism on the B c can be estimated from a potential model in the same way as we have done for bottomonium and charmonium. The quark and antiquark in the B c have the same sign of electric charge, however, and so the effect now is repulsive rather than attractive. We estimate that the effect of switching on electromagnetism is to shift our B c mass upwards by 2 MeV. The effects of c quarks in the sea can be estimated following the discussion in subsection II C as approximately 1/60 of the hyperfine splitting in the B c system, or 1 MeV. This effect is attractive and so will counteract the effect of electromagnetism. We take the net shift in the B c mass as 1(1) MeV, moving our result to 6.280(4) GeV. This is the value given by the dark shaded band in Figure 7.
The complete error budget is given in Table V . We must estimate the effect of these terms on both M bb and M Bc . v 2 b is about half the size in the B c compared to bottomonium. This means that there is some cancellation of the α s v 4 b errors in ∆ Bc,hh , since we estimated α s v 4 b errors previously as α s v 2 b (500 MeV). Independent confirmation that these terms have a small net effect comes from the calculations that we have done here for the case where c 1 , c 5 and c 6 take the values that include the O(α s ) radiative corrections, see Table VI and Figures 5 and 6. There will be little cancellation of the v 6 b errors, however, since v 6 b is much smaller in the B c than in bottomonium. We take the systematic error from spin-independent terms then to be 4 MeV for α s v 4 b (i.e. half that for B s ) and 5 MeV for v 6 b , added in quadrature to give 6 MeV. The leading spindependent error is from missing radiative corrections to the σ.B term. This affects the B c only because of the spin-averaging of M bb . As for the B s we take 3/4 of a possible 10% error in the hyperfine splitting, estimated in [22] at 53 MeV, i.e. 4 MeV. Errors from the uncertainty in M bb and M ηc do not affect ∆ Bc,hh significantly but M Bc inherits an error of half their uncertainty when it is reconstructed from ∆ Bc and those masses. These two uncertainties are partly correlated, because they both contain estimates of electromagnetic and annihilation effects in the two very similar charmonium and bottomonium systems. The systematic errors from electromagnetism and c quarks in the sea for the B c are also correlated with the errors for these effects from charmonium and bottomonium. These errors are marked with a star in Table V and Figure 9 for the equivalent as a function of M bb and Figure 10 for the equivalent as a function of Mη c . out the different components and add them linearly with appropriate signs before squaring and accumulating into the total error. The error from these three components is then increased by their correlation from 3 MeV to 4.6 MeV. We estimate finite volume errors to be negligible for the B c . Our total error, from adding statistical and systematic errors in quadrature, is 9.5 MeV, giving a mass for the B c meson from the hh method of 6.280(10) GeV. The total error is plotted as the more lightly shaded band in Figure 7.
The hs method has different systematic errors from the hh method and so provides a good cross-check. The raw results needed for this method are given in Table VI and in Figures 8, 9 and 10 we show results for the quantity −∆ Bc,hs (because ∆ Bc,hs is negative) defined by: In the figures −∆ Bc,hs results from set 1 are plotted against the different quark masses involved in the calculation, with M 2 ηs , M bb and M ηc acting as proxies for the s, b and c quark masses respectively (we have two different values for the masses of each quark). We see that there is fairly strong dependence on the s quark mass but very little on the b quark mass or the c quark mass. The slope against M 2 ηs is 0.41, which agrees well with that expected from experiment if we substitute light quarks for s quarks in the formula for ∆ Bc,hs above. The slope against M bb is 0.005 and against M ηc is 0.07. These are only in very rough agreement with the expectations    Figure 8 for the equivalent as a function of M 2 ηs and Figure 9 for the equivalent as a function of M bb . As before we can use these results to compensate for mistuning of the quark masses. Again we take a 50% error on tuning shifts for b and s but a 200% errors on those for c to allow for discretisation errors in our estimates of those effects. The tuned results are given in Table VII. All of the statistical, tuning and lattice spacing errors are larger than those of the hh method. The tuning error dominates on the very coarse and coarse lattices, but on the fine ensemble it is comparable to the other errors. The lattice spacing error is reduced by a factor of two over the naive error by the retuning required when the lattice spacing changes. Lattice spacing dependence is small but visible in these results; no sea quark mass dependence is evident.
As for the B s and the B c hh method we fit the results for ∆ Bc,hs to a functional form allowing for lattice spacing dependence, including NRQCD effects, and sea quark mass dependence. The functional form is the same as that used for the hh method, except that now we can include these dependences as a multiplicative factor since ∆ Bc,hs is not unusually small. We use: ∆ Bc,hs (a, δx l , δx s ) = ∆ Bc,hs,phys 1 + (23) We take the same prior values as for the hh method except for the prior for ∆ Bc,hs,phys which we take to be -1.0 (2).
We obtain the result ∆ Bc,hs,phys = -1.054(17) GeV. This changes by less than 1 MeV on doubling the prior width for the lattice spacing dependence or the sea quark mass dependence. To reconstruct the B c mass from this we need to add the appropriate values for the B s and D s masses in a world without electromagnetism or c quarks in the sea. As discussed earlier, electromagnetism has negligible effect on the B s mass. The D s mass is lower by 1.3(7) MeV, however, in a world without electromagnetism, from the phenomenological formula in equation 17. This gives a total for the appropriate sum of M Ds + M Bs of 7.334 GeV. Figure 11 shows our tuned results for m Bc as a function of lattice spacing. The dark shaded band is the result from the fit just described including the error obtained from it. We have shifted the B c mass obtained upwards by 1(1) MeV, to a central value of 6.281 GeV, as described earlier to allow for electromagnetic and charm-in-the sea effects that are not included in our calculation. The lighter shaded band gives the total error, from the error budget of Table V, the systematic error components of which we will now discuss. The first three entries in the Table are the split of the 17 MeV error obtained from the fit among its different components.
Errors from missing relativistic corrections to the NRQCD action affect the B c energy and the reference B s energy. The leading missing spin-independent terms in the NRQCD action are O(α s v 4 b ) and O(v 6 b ). We earlier estimated the shift from O(α s v 4 b ) terms on M bb at 15 MeV. For the B c we expect a systematic error of about half this value, so we take 7.5 MeV, since v 2 b is roughly half as big. For B s α s v 4 b terms have very little effect and neither meson will be sensitive to v 6 b terms. Spin-dependent NRQCD errors come chiefly from missing radiative corrections to the σ · B term, but there will be cancellation here between the B c and the B s since both mesons will respond in the same way to a change in c 4 . We therefore take a systematic error of 1 MeV rounding up the difference between the 4 MeV systematic previously allowed for this for the B c and the 3.5 MeV systematic for the B s .
Systematic errors from uncertainties in the physical values of M ηs , M ηc and M bb which affect the quark mass tuning can be estimated from the dependence of ∆ Bc,hs on these quantities discussed earlier. The error from the uncertainty in M ηs is sizeable at 2.3 MeV; the others are very small. We must also allow for systematic errors from uncertainties from electromagnetism and charm-inthe-sea for the reference masses of the D s (0.7 MeV is half the shift applied in that case) and the B s (negligible) as well as for the B c itself (1 MeV as discussed above).
This gives a total error of 19 MeV, dominated by the statistical and tuning errors of the raw data. Our final result for the B c mass from the hs method is then 6.281 (19) GeV. This is plotted as the more lightly shaded band in Figure 11. The agreement between the hs and hh methods is very good, although their systematic and statistical errors are very different, with the hh method being significantly more accurate. The agreement is in fact not surprising when we consider that the B s mass determined in section III A agreed well with experiment. The B c hs method replaces M bb with M Bs and M ηc with M Ds so if the B s and D s masses are known to agree with experiment given masses tuned from M bb and M ηc then the B c from hs and hh will agree. However, the fact that they were derived completely independently is a good consistency check of the method and of our error estimates.

C. B mass
The correlators for the B meson are noisier than those for the B s , as will be clear from the discussion of the signal/noise earlier. This means that the B meson mass is the least well determined of all our masses. The best way then to pin down the B mass is to consider the mass difference between the B s and the B. NRQCD systematic errors should entirely cancel in such a difference. However, because the difference is a small number we have a fairly sizeable statistical error even when we fit both correlators together, as described earlier, and extract E Bs − E B directly from the fit. Table VIII gives values for the energy difference between B s and B extracted from our fits on each ensemble. Statistical errors are 10-15% of the splitting. However, this amounts to less than 10 MeV in terms of the absolute mass, so still provides a good test against experiment for M B .
We plot the results for M Bs −M B (= a −1 (aE Bs −aE B )) against M 2 ηs − M 2 π which is a useful physical proxy for m s − m l in Figure 12. We expect this mass difference to be largely linear in m s −m l and our results are consistent with this. Given the statistical errors, we fit a relatively simple form to this difference: Here the functions f j are simple ones that respect the fact that M Bs − M B vanishes by definition when the light quark mass is equal to the strange quark mass. So π log M 2 π . We allow these terms to have lattice spacing dependence with a scale set by Λ = 400 MeV. We also allow sea quark mass dependence in the terms multiplying f 1 . The coefficients a j are given priors of 0.0(5) (we expect a slope of 0.2 if the dependence on M 2 ηs −M 2 π were purely linear). For the c j2 a 4 -dependence coefficients we take 0.0(1.0) and for c j1 we take 0.0(5) since a 2 terms should be suppressed by an additional power of α s . For the sea quark mass dependent coefficients b l and b s we take priors of 0.00(7) as discussed earlier.
The physical value for M Bs − M B is then the value at M ηs = 0.6858 GeV and M π 0 = 0.135 GeV, for a = δx l = δx s = 0. We obtain 0.073 (14) GeV. This value is to be compared to the experimental mass difference between the B s meson and the average of the charged and neutral B mesons which we denote B l (thus averaging the u and d quark masses). However it has to be adjusted for electromagnetic effects not included in our lattice QCD calculation. Following the discussion in section III A we see that electromagnetic effects in the B s and B d mesons are very small. For the B u however, because it is charged, the shift is more substantial at 2.2 MeV [39]. Adding in electromagnetism then shifts our M Bs − M B l splitting down by 1 MeV. The result 72 (14) MeV is shown as the shaded green band in Figure 12. It is in reasonable agreement with the experimental result of 87 MeV [20].  Results for energies and masses needed for the determination of the mass of the B meson. Column 6 gives the energy difference between Bs and B l mesons, for different valence b, s and l quark masses given in columns 2, 3 and 4. Column 5 gives the corresponding π meson mass, taken from [9]. Column 7 gives the hyperfine splitting for the B l meson, from [22], Our final result for M B l = (M B ± + M B 0 )/2 is 5.363 -0.072 = 5.291(11)(14) GeV. The first error comes from the mass of the B s and is discussed in subsection III A, the second comes from the mass difference between the B s and B l . We do not expect any significant additional systematic errors from NRQCD, beyond those that the B l inherits from the B s in this method, because those errors should cancel in M Bs − M B l . The error budget for M B l is then as given for M Bs in Table V  By projecting out the vector states at the source and sink we can measure the correlator for the B * , B * s , and B * c . As they come from exactly the same configurations and valence HISQ propagators as the corresponding pseudoscalar states, they are highly correlated. In this case we can do simultaneous fits to both the pseudoscalar and vector meson propagators, and extract a value for the This hyperfine splitting is generated by the σ · B term in the NRQCD action, equation 6. This is O(v 4 b ) in the relativistic power counting for heavyonium and O(Λ/M b ) in heavy-light power counting. In our action we take the coefficient of this term, c 4 to be 1, but it will have radiative corrections when matched through O(α s ) with full QCD. We are also missing higher dimension operators that correct for discretisation errors and add relativistic corrections. For B systems, which are relatively large with very slow-moving b quarks, we do not expect these latter effects to be as important as the issue of the determination of c 4 beyond tree level. The heavy-light hyperfine splitting generated by the σ · B term is proportional to c 4 and so uncertainty in c 4 leads directly to an O(α s ) i.e 25%, uncertainty in the splitting which decreases only slowly on finer lattices. Thus to determine this splitting accurately we need a determination of c 4 .
Since we use exactly the same NRQCD action for all our calculations, however, we can effectively determine c 4 by comparing one set of heavy-light hyperfine splittings to experiment and then predicting the others. Equivalently we can take ratios of hyperfine splittings in which the normalisation factor, c 4 , cancels. This is what we did in [22]. By using the B * s − B s mass difference, which is 46.1(1.5) MeV from experiment [20], we showed that this splitting does not depend on the mass of the lighter quark even for as heavy a quark as the charm quark, and we were able to predict a B * c − B c splitting of 53 (7) MeV. We will not discuss that analysis further here, but we give the table of results of the hyperfine splittings for completeness in Tables III, VI, VIII. They include some additional values over those in [22] for the purposes of further testing systematic errors. Figure 13 shows such a test in a plot of the mass difference between B * s and B s as a function of M bb for two different s quark masses (the first 4 rows of entries in Table III). Dependence on the b quark mass is visible, but dependence on the s quark mass is very small. Results are also shown for B * c − B c for the same b quark masses and they show a parallel slope. In addition we show a result for the the NRQCD action with c i coefficients different from 1 and no change is seen.
Here we are interested in analysing the systematic error in the B, B s and B c meson masses from any uncertainty in c 4 . We do this by comparing our B * s − B s splitting to experiment and interpreting any mismatch as a signal for c 4 = 1. We use the B * s − B s because this is the most accurately determined splitting from our analysis that is also known experimentally. This method can be used as a nonperturbative determination of c 4 , and we used this previously to bound the errors on our prediction of the hyperfine splitting in bottomonium based on earlier B s and B hyperfine splitting results [11]. Figure 14 shows our results as a function of lattice spacing. We have adjusted them for mistunings of the b quark mass according to the results in Figure 13 but the corresponding shifts are small, and less than the statistical errors in all cases. The results show little sign of any lattice spacing dependence or sea quark mass dependence and we see that a value of c 4 of ≈ 1.1 would give agreement with experiment for all the values. We therefore estimate that the correct value of c 4 for this NRQCD is 1.1(1) and that we make a systematic error of about 10% in the heavy-light hyperfine splitting by using c 4 = 1. This produces a systematic error on the B c and B s meson masses discussed in earlier subsections and included in Table V.
Note that the behaviour of the hyperfine splitting in bottomonium is quite different from that of the B s , being strongly dependent on the lattice spacing [11]. However, it is the same operators in the NRQCD action, with the same coefficients, that control the fine structure in both systems. The matrix elements of the operators can behave quite differently, and bottomonium is expected to be a lot more sensitive to the lattice spacing than the B s . This means that the B s is a good system from which to determine c 4 because it is really only sensitive to that coefficient.

E. Scalar and Axial vector meson masses
When generating the NRQCD propagators we choose Dirac structures Γ to explicitly project out pseudoscalar and vector mesons. Parity partners of both of these contribute to their correlators, as shown in equation 12, and must be included in the fit. The parity partner state of the pseudoscalar is a scalar meson and the vector meson has as its parity partner a axial-vector state. So, by carefully fitting the correlators of the 0 − and 1 − states we get also the spectra of the 0 + and 1 + states for free.
In fact our fit results return directly the mass difference between the ground state in the oscillating channel and the ground state in the non-oscillating channel i.e the 0 + − 0 − and 1 + − 1 − mass differences. We report these results in Table III for the B s (i.e. for the B * s0 and B s1 mesons). For B l and B c our errors are too large on some fits to give a full picture across all ensembles.
The results for the scalar-pseudoscalar mass difference are shown in Figure 15. There is no signal for any systematic dependence on the b or s quark mass, or on the lattice spacing. In deriving a physical result we allow for both physical and unphysical lattice spacing dependence as described for the B s mass in subsection III A, as well as sea quark mass dependence. We use the same fit form as for the B s mass, given in equation 12. The priors are taken to be the same except that we take the prior on the physical value of the scalar-pseudoscalar mass splitting to be 0.4 (2). We also allow for more sea quark mass dependence than in that case, because the scalar meson is not gold-plated (this will be discussed further below). We therefore do not take the factor three suppression of sea quark mass effects in this case, so the prior on the sea quark mass dependent terms is simply 0.0(2) for the linear terms and 0.00(4) on the quadratic terms. The physical result we obtain is 0.385 (16) GeV and this is plotted as the shaded band on Figure 15.
Exactly the same procedure is followed for the axial vector -vector splitting. The results are plotted in Figure 16. From the same fit as that described above we obtain the physical result for the mass difference between the axial vector and vector of 0.391 (15) GeV, plotted on the Figure as a shaded band.
Since we have calculated the mass differences between the scalar and axial vector B s mesons and the corresponding pseudoscalar and vector B s mesons we expect only very small systematic errors coming from NRQCD. Because the b quark is very nonrelativistic in these systems, as discussed in subsection III A the errors from missing higher order relativistic corrections in NRQCD are very small. They will be reduced further here by cancellation in the mass difference. The main source of systematic error from NRQCD will come from radiative corrections to spin-dependent terms in the NRQCD action. In subsection III D we showed that these systematic errors are not large, at least for the σ · B term. There the errors amounted to 10% of the hyperfine splitting, around 5 MeV, split between the vector and pseudoscalar states. Assuming a similar error for other spin-dependent terms which would affect p-wave states, we take a systematic error of 5 MeV from NRQCD in the mass differences.
A potentially much larger source of systematic error is the fact that the scalar and axial vector mesons have strong decay modes, i.e. they are not 'gold-plated'. This will be discussed further in the next section. The strongest decay mode, if kinematically allowed, will be to BK (for the B * s0 ) or B * K (for the B s1 ). If the masses are such that the mesons are below threshold for this decay mode, there will still in principle be coupling between the meson and this virtual decay channel which can shift the meson mass. There is in addition a Zweig-suppressed decay mode to B s π/B * s π which will be kinematically possible.
On the lattice the coupling between single and multiparticle states is distorted by the fact that u/d quark masses are heavier than their physical values and the volume of the lattice is relatively small. The fact that u/d masses are unrealistic means that decay thresholds are higher than in the real world. In principle sensitivity to decay thresholds would be seen in the results as sea quark mass dependence, but that may not become visible until much closer to real world u/d mass values. The finite volume of the lattice restricts the decay momenta that real or virtual multiparticle states can have. A lattice analysis on multiple volumes allows single and multiparticle states to be separated. In practice [31] it seems that bilinear operators of the kind that we have used here have very small overlap with multiparticle states. So, although in principle there may be a multiparticle state (such as B s π) at a lower mass value than the B s0 it is very hard to pick it out of a lattice QCD calculation without explicitly using multiparticle operators, which we have not done.
A simple model to analyse the effect of multiparticle states is one in which point-like meson states are coupled together via a perturbation which is a simple point-like vertex. We can then calculate the shift on the single particle energy from this coupling by integrating over the momenta of the decay products in the initial particle rest frame. For example, for B * s0 coupling to BK: Λ represents an ultraviolet cut-off required for this model to make sense. Λ ≈ 500 MeV and g 2 ≈ (0.5/Λ). If our calculation is correct that the B * s0 is close to, but below, threshold then we can treat the B and K as nonrelativistic and, dropping the B kinetic term, where ∆M = M B + M K − M B * s0 (unperturbed values). The ∆M -dependent piece of the mass shift is then given by: Numerically this gives a shift downwards of a few tens of MeV for ∆M values of a few tens of MeV. From this we conclude that a reasonable systematic error for the absence of coupling to strong decay channels is 25 MeV (which we take to be a symmetric error). This then gives the following mass differences: where the first error is statistics/fitting, the second is the NRQCD systematic error and the third is the error from not including coupling to strong decay channels.
IV. DISCUSSION Figure 4 shows that our result for the mass of the B s meson agrees well with experiment with total errors of 11 MeV (0.2%). The errors are dominated by statistical errors and systematic errors from NRQCD, both of which are being improved in work underway.
As discussed earlier, because we fix the b quark and s quark masses from other mesons, the B s mass determination is completely free from any parameter tuning. An alternative for the b quark mass, adopted by some other lattice QCD calculations is to fix the b quark mass from the B s meson mass itself. However, it is still possible then to determine ∆ Bs = M Bs − M bb /2, as a test of the b quark systematic errors. The only other full lattice QCD calculation of this quantity is from the Fermilab Lattice/MILC collaboration using the Fermilab formalism for the b quark [25]. They determine in fact the quantity ∆ Bs = M Bs − M bb /2 where B s is the spin average of the B s and the B * s masses. This quantity has reduced systematic errors from the spin-dependent terms in the action, in the same way that the use of M bb reduces these systematic errors for the bottomonium system. The Fermilab Lattice/MILC collaboration obtain the value 1359 ± 304 +31 −0 MeV for 2∆ Bs with a partial error budget [25]. We can also determine ∆ Bs in exactly the same way as we determined ∆ Bs . We obtain 0.671(7) GeV for the physical result from our calculation. This becomes 0.675(11) GeV when corrected for electromagnetic, annihilation and charm-in-the-sea effects in bb and with a full error budget (essentially the same as in Table V but with a reduced error for NRQCD systematics in the B s ). The experimental result is 0.6817(11) GeV [20]. Figure 17 shows the results for ∆ Bs and ∆ Bs from this paper and from the Fermilab Lattice/MILC collaboration compared to experiment. Both results agree with experiment but we are able to provide a 2% test of these mass differences, which is a nontrivial test of QCD.
Back in 2004 we predicted the mass of the B c ahead of the CDF experimental discovery in a collaboration with the Fermilab Lattice collaboration [26]. We used NRQCD for the b quarks, as here, but the Fermilab formalism for the c quarks, and the asqtad formalism for the  s quarks in the hs method. As a result, we had larger statistical and systematic errors than we have here, particularly for the hs method. Figure 18 shows the comparison between our old results and the new ones given here, as well as the current experimental value. The improvements in lattice QCD calculations since 2004, including the development of the HISQ action for c and s, give us a substantial improvement in errors and consistency between the hh and hs methods today.
Our result for the mass difference between the B s and the B l meson is the first full lattice QCD calculation of this quantity. As discussed in subsection III C our result is in agreement with experiment, but with substantial statistical errors. These will be improved in further work which is under way.
In subsection III E we gave results for the masses of the 0 + and 1 + B s mesons. The 0 + has not been seen experimentally. A 1 + state has been seen but may not be the one whose mass we have calculated. Figure 19 shows how our results fit into the current experimental picture of 'p-wave' charm-light and bottomlight mesons. Charmonium and bottomonium p-wave mesons are also shown for comparison [20]. For heavylight mesons the p-wave states are expected from heavy quark symmetry [27] to appear in two doublets, classified according to the J of the light quark which can be either 1/2 or 3/2, when L=1 is combined to s l = 1/2. The j l = 1/2 doublet then separates into a 0 + and 1 + meson at non-infinite heavy quark mass, when the heavy-quark spin is coupled in. The j l = 3/2 doublet is likewise made up of 1 + and 2 + states. This is in contrast to the heavyonium case where there is a triplet of 0 + , 1 + and 2 + states with total quark-antiquark spin of 1, and a single 1 + state with total spin 1. The existing experimental results are shown as solid points in Figure 19 divided appropriately according to the picture above. For D and D s mesons both doublets have been seen; for B and B s mesons only the j l = 3/2 doublet has been seen (assuming that the 1 + state seen is associated with that doublet). For charmonium all 4 states of the χ c triplet and the h c are known; for bottomonium the h b has not been seen. The experimental masses are given relative to the spin-average of the s-wave states, a pseudoscalar and a vector in every case. That removes the overall mass scale of each system from the plot and shows, as is well-known but still somewhat surprising, that the orbital excitation energies of heavy degrees of freedom in heavyonium are very similar to those of light degrees of freedom in a heavy-light system.
Since mass splittings between the S = 1 states in heavyonium and between the members of the j l = 1/2 or j l = 3/2 doublets in the heavy-light case are caused by heavy quark spin effects proportional to the inverse of the heavy quark mass we expect to see larger splittings in the c case than in the b case. This is borne out in the experimental data for charmonium and bottomonium and in the comparison of D and B results for the j l = 3/2 doublet (although the disagreement between B and B s might indicate that the doublet assignment for the B 1 + in the Figure is wrong). The splitting between j l = 1/2 and 3/2 doublets is a light quark effect that does not vanish as m Q → ∞. However the splitting will vary with m Q slightly because of Λ/m Q terms in the effective heavy quark action (NRQCD) away from that limit. A variation in the splitting of order 100 MeV out of 500 MeV is then reasonable between D and B.
Our results are entered on Figure 19 as open circles in the B s column. Since we have 0 + and 1 + states we have placed them as the j l = 1/2 doublet. However, it should be stressed that we do not know that that assignment for the 1 + is correct. In any case the 1 + states from the two doublets can mix and we have not allowed for that.
Given the discussion above, our results fit fairly naturally into the picture described. As a j l = 1/2 doublet they sit below the known j l = 3/2 doublet. They sit closer to the j l = 3/2 doublet than for the D s case, but this can be a Λ/m Q effect as discussed above. The splitting between the two states in the doublet is about one third that in the D s case, consistent with this splitting being a 1/m Q effect.
In the D s case the discovery of the lower 0 + /1 + doublet [28][29][30] caused much surprise because the states were low compared to model calculations. The states had been expected to be above threshold for strong decay to DK and D * K respectively and therefore broad (unlike the upper doublet which has to decay in a d-wave). Instead they are below threshold and so decay to the Zweig-suppressed D s π and D * s π channels and are narrow. We have marked the DK, D * K, BK and B * K thresholds on Figure 19. Like the D * s0 our 0 + B s state is also below, but very close to, its Zweig-allowed decay threshold. A similar situation holds for the 1 + state. This might indicate that these states would be narrow. However, there will also be effects from coupling to the decay channel that are not included in our calculation. On our lattices the light sea quark masses are heavier than in the real world and hence the K mesons containing a valence s quark and a sea light quark would be too heavy to allow the B * s0 to decay to BK. B s π decay is allowed but with a very restricted phase space compared to the real world. Coupling to these channels would in principle show up as sea quark mass dependence but would need a bigger range of sea quark masses than we have used. In subsection III E we allowed a 25 MeV systematic error for these coupled-channel effects, noting that the coupling to BK decay will tend to push the mass down.
Our results are the only ones in full lattice QCD to date and with realistic b quarks. There have, however, been several recent lattice QCD results for the case of u and d sea quarks only and taking b quarks in the static limit [31][32][33][34]. The most complete is that of the ETMC collaboration [34]. They give a mass difference between the B s j l = 1/2 doublet and the spin-average of B s and B * s of 413(12) MeV, with an estimated additional possible systematic error of 20 MeV, including coupling to multiparticle states. Using experimental results from charmed mesons to estimate 1/m b corrections to the static limit, they conclude, as we have done, that the scalar B s state is close to the BK threshold.

V. CONCLUSIONS
We have given the first accurate result for the B s meson mass from lattice QCD including the effect of u, d and s sea quarks, and with a full error budget. We have improved significantly on an earlier value for the B c meson mass, achieving smaller errors and better consistency between two different methods. The determination of both of these masses provides a strong test of our lattice QCD approach to b physics, because they test the consistency of heavyonium and heavy-strange or heavycharm physics from the same heavy quark action. All of the QCD parameters used here are tuned from other calculations so our results are parameter free tests of QCD against experiment.
The mass of the B meson, specifically the difference between B s and B meson masses, depends on light quark physics since heavy quark effects cancel. Our result agrees with experiment but needs higher statistical precision for a good test.
We also discuss scalar and axial vector meson masses for the B s . Our results indicate masses below, but close to, threshold for Zweig-allowed decay modes. From our current calculation, however, it is not possible to include effects of coupling to either allowed or suppressed decay channels, so significant shifts to our results from these effects are possible.
Further improvement to these results will come with improved statistical accuracy in calculations now underway. This will lead also to improved determination of decay constants and other B meson matrix elements. Confidence in those calculations and the error analysis associated with them is strongly bolstered by this analysis of the associated meson masses.