New methods for B meson decay constants and form factors from lattice NRQCD

We determine the normalisation of scalar and pseudoscalar current operators made from non-relativistic $b$ quarks and Highly Improved Staggered light quarks in lattice Quantum Chromodynamics (QCD) through $\mathcal{O}(\alpha_s)$ and $\Lambda_{\text{QCD}}/m_b$. We use matrix elements of these operators to extract $B$ meson decay constants and form factors, then compare to those obtained using the standard vector and axial-vector operators. This provides a test of systematic errors in the lattice QCD determination of the $B$ meson decay constants and form factors. We provide a new value for the $B$ and $B_s$ meson decay constants from lattice QCD calculations on ensembles that include $u$, $d$, $s$ and $c$ quarks in the sea and those which have the $u/d$ quark mass going down to its physical value. Our results are $f_B=0.196(6)$ GeV, $f_{B_s}=0.236(7)$ GeV and $f_{B_s}/f_B =1.207(7)$, agreeing well with earlier results using the temporal axial current. By combining with these previous results, we provide updated values of $f_B=0.190(4)$ GeV, $f_{B_s}=0.229(5)$ GeV and $f_{B_s}/f_B = 1.206(5)$.


I. INTRODUCTION
Hadronic weak decay matrix elements containing bquarks that are calculated in lattice Quantum Chromodynamics (QCD) are critical to the flavour physics programme of overdetermining the Cabibbo-Kobayashi-Maskawa (CKM) matrix in order to find signs of new physics. The accuracy of the lattice QCD results often limits the accuracy with which the CKM matrix elements can be determined and with which the associated unitarity tests can be performed [1]. It is therefore important both to improve and to test the accuracy of the lattice QCD results. This includes determining the lattice QCD values using a variety of different formalisms for b quarks and light quarks, in addition to using different methodologies within a given formalism.
It is now becoming possible to study heavy quarks up to the mass of the bottom quark using relativistic formalisms [2,3], but this is relatively expensive numerically. Consequently, to date, the most extensive studies of heavy quarks in lattice QCD have been done with nonrelativistic formalisms, such as NRQCD [4] or the Fermilab formalism [5] and its variants [6]. Relativistic formalisms have the advantage of simple continuumlike current operators that couple to the W boson that can be chosen to be absolutely normalised, for example through the existence of a partially conserved axial current (PCAC) relation [7]. The main issue with these * chughes@fnal.gov † christine.davies@glasgow.ac.uk ‡ URL: http://www.physics.gla.ac.uk/HPQCD formalisms is then controlling discretisation errors [8]. In nonrelativistic formalisms the numerical calculation itself is more tractable, along with the control of discretisation errors, but the current operators have a nonrelativistic expansion and must have their normalisation matched to that of the appropriate continuum operator. The expansion and the normalisation are the main sources of systematic uncertainty in these lattice QCD results. The comparison of lattice QCD values derived using nonrelativistic and relativistic formalisms provides a test of systematic uncertainties (see, for example, [9] and [3]). However it is also important to provide tests of systematic uncertainties within a given formalism using different methods. Here we provide such a test of the NRQCD approach by normalising new sets of current operators that have not been used in this formalism before, then comparing results for the decay constants and form factors obtained to the previous determinations.
The archetypal heavy meson weak decay process is annihilation of a B meson to τ ν. The hadronic parameter which controls the rate of this process is the B meson decay constant, proportional to the matrix element to create a B meson from the vacuum with the temporal axial current containing a heavy quark field and a light antiquark field. The most precise calculation to date of the B meson decay constant, f B , uses improved lattice NRQCD and Highly Improved Staggered light quarks on gluon field configurations that include u, d, s and c quarks in the sea with multiple values of the lattice spacing and a u/d quark mass going down to the physical point [9]. That calculation used lattice QCD perturbation theory [10] to normalise the temporal axial current operator through O(α s ), O(α s Λ QCD /m b ) and O(α s aΛ QCD ) and obtained a final uncertainty of 2%, including uncertainties from current operator matching and missing higher order current operators.
Decay constants can also be defined in continuum QCD from pseudoscalar current operators using the PCAC relation. This is typically the method of choice for lattice QCD calculations using relativistic formalisms where a lattice PCAC relation allows the pseudoscalar current to be absolutely normalised. This enables the D and D s decay constants to be obtained with 0.5% uncertainties using the HISQ formalism [7,11,12]. Here we normalise the NRQCD-light pseudoscalar current through O(α s ), O(α s Λ QCD /m b ) and O(α s aΛ QCD ) and obtain a value for f B with similar uncertainty to that determined from the temporal axial current, providing a test of the systematic errors.
B meson exclusive semileptonic processes are important for the determination of CKM matrix elements through the matching of experimental decay rates to theoretical expectations as a function of momentum transfer.
Here the hadronic parameters that encapsulate the information needed on QCD effects are the form factors, calculable in lattice QCD. For the case in which both initial and final mesons are pseudoscalars (e. g. B → π ν) there are two form factors, a vector form factor and a scalar form factor. It is the vector form factor that gives the decay rate in the light lepton mass limit, but both form factors appear in the lattice QCD determination of the matrix elements of the vector current. The form factors can be separated by comparing spatial and temporal vector current matrix elements but additional information can also be obtained by determining the scalar form factor directly from the scalar current. Indeed this method has been used for the accurate determination of D and K meson semileptonic form factors in lattice QCD using the HISQ formalism [13][14][15]. Here we compare results using NRQCD-light scalar currents to those obtained using vector currents for B → π ν. We discuss how this method will be used in improved 'second generation' B meson semileptonic form factor calculations now underway.
The paper is organised as follows: in section II we derive the normalisation of the NRQCD-light scalar and pseudoscalar current operators; in section III we combine this with the lattice calculation of the matrix elements of different components of the current to give results for decay constants and form factors; section IV gives our conclusions, including planned future work using these results.

II. NORMALISATION OF LATTICE NRQCD CURRENT OPERATORS
Here we discuss the normalisation of the lattice NRQCD-HISQ current operators when the light quark is taken to be massless and follow the methodology laid out in [16,17], along with most of the notation. We will start with a discussion of the temporal axial cur-rent and show the modifications that need to be made to those results to yield the normalisation of the pseudoscalar current. Results for the temporal vector/scalar case are then identical because of the chiral symmetry of the HISQ action.
The matrix element of the appropriate temporal axial current defined in continuum QCD between the vacuum and pseudoscalar meson, H, at rest yields the meson decay constant, f H , via the relation where M H is the meson mass. The continuum QCD current operator can be systematically expanded in terms of lattice NRQCD-HISQ current operators (whose matrix elements can be determined in a lattice QCD calculation) as where increasing j corresponds to operators that are higher order in a relativistic expansion. The C j are dimensionless coefficients that compensate for the different ultraviolet behaviour between the continuum and lattice regularisations of QCD and hence they can be calculated in perturbation theory as a power series in the strong coupling constant, α s . The coefficients of powers of α s will depend on the bare heavy quark mass in lattice units, am b , which is the parameter appearing in the lattice NRQCD action (we use b rather than the generic label h since this is almost always the b quark). Here we work through O(α s ) and include the three operators (j = 0, 1, 2) that allow us to match the current through O(α s Λ QCD /m b ) and O(α s aΛ QCD ). The determination of the C j is done most conveniently by choosing to match matrix elements of the left-and right-hand sides of eq. (2) for a heavy quark to light quark scattering process induced by the current. The procedure then [16,17] is to: • calculate the amplitude for such a process through O(α s ) in continuum QCD; • expand this amplitude through first order in powers of 1/M where M is the heavy quark pole mass; • choose lattice NRQCD-HISQ operators that reproduce the terms in this expansion and calculate the one-loop mixing matrix of these operators in lattice QCD perturbation theory using the same infrared regulation procedure as used in the continuum. Infrared divergences must cancel between the continuum and lattice calculations in the end, since the two only differ in ultraviolet physics. Note also that the mixing matrix should be calculated at a pole mass that matches that of the continuum calculation; • invert this mixing matrix to determine the (finite) C j coefficients that will give the correct linear combination of lattice currents to produce the same one-loop scattering amplitude as in continuum QCD.
The continuum calculation for the temporal axial current q(x)γ 5 γ 0 h(x) was done in [16] in the M S scheme using Feynman gauge, on-shell mass and wavefunction renormalisation and a gluon mass (λ) to regulate infrared divergences. q(x) is the light quark field and h(x) the heavy quark field satisfying the Dirac equation and the γ matrices are the standard ones in Euclidean space-time.
The key diagram to be calculated in continuum QCD is shown in Figure 1, where the double line represents an incoming heavy quark of momentum p, the single line an outgoing massless quark of momentum p and the cross represents the current. The self-energy diagram must also be evaluated to determine the wave-function renormalisation. The result for the temporal axial current amplitude is given through 1/M as a combination of five matrix elements of Dirac spinors multiplied by factors of p 0 , p 0 and p · p in [16]. By using the Dirac equation for the light quark, and expanding the heavy quark energy and Dirac spinor to 1/M , this is reduced to The coefficients are  (15)) needed to determine the renormalisation of the lattice NRQCD-HISQ pseudoscalar/scalar current from that of the temporal axial vector/temporal vector current for massless HISQ quarks. Column 5 gives the one-loop NRQCD mass renormalisation coefficient. These results were calculated and presented as the linear combination relevant for eq. (18) in [10] using the standard v 4 -accurate NRQCD action (with stability parameter n = 4) and the individual values are given here. We also include new results for a lighter b quark mass, am b = 1.22, suitable for the MILC 'superfine' (0.06 fm) lattices. The constituent matrix elements are where u Q is a two-component spinor related to the Dirac spinor u h (p) (to O(1/M 2 )) by and where u Q satisfies γ 0 u Q (p) = u Q (p). We now carry out the continuum calculation to the same order for a pseudoscalar current P = q(x)γ 5 h(x) and obtain where u q and u h are Dirac spinors and Here µ is the scale parameter from dimensional regularisation and λ is the gluon mass. We expand the heavy  [9] from numbers in [10] and are reproduced here. The results for z A0 2 have changed slightly for the heaviest masses because of an improved calculation of ζ02.
with η (j) P defined in an analogous way to eq. (4) and Note that B (0) P and B (1) P have the same value ((a 1 − 1)/α s ) coming from the first term on the right-hand side of eq. (8). A check on these results comes from applying the continuum PCAC relation. This shows that the leading order term B (0) should differ between P and A 0 by an amount that is the one-loop conversion factor between the pole and M S quark mass at scale µ. Using γ 0 u Q = u Q the relationship between the operator matrix elements for the pseudoscalar and temporal axial current cases are with a sign change for the leading relativistic correction, j = 1. Exactly the same relations are obtained for the scalar case with respect to the temporal vector calculation.
.297 -0.0008 (22)  The current operators needed in the lattice NRQCD calculation are readily identified from the Ω (j) by replacing spinors with fields and converting momentum factors to derivatives. This gives, for the temporal axial-vector case, where m b is the bare lattice NRQCD quark mass and Q(x) is the two-component NRQCD field (i.e. a fourcomponent field with zero in the lower two-components). The analogous expressions for the pseudoscalar case mirror eq. (12). The next step is to calculate the mixing matrix for the lattice operators in lattice QCD perturbation theory through one-loop. For the A 0 case with Z A0,ij written as [16,17] (15) Z q is the coefficient of the one-loop term in the wavefunction renormalisation for massless lattice quarks, here in the HISQ formalism. Similarly, Z b is the coefficient of the one-loop term in the lattice NRQCD wavefunction renormalisation and Z m b the coefficient of the one-loop mass renormalisation between the bare NRQCD quark mass and the pole mass [18]. This latter factor appears for j = 1, 2 because of the explicit mass factor in the operator and our choice to use the bare NRQCD mass in the The different contributions that make up the renormalisation coefficients z A 0 0 and z P 0 as a function of the bare heavy quark mass in the lattice NRQCD Hamiltonian, as given in eqs. (21) and (24). The infrared-finite pieces are plotted for infrared-divergent contributions and Zq is not included since it does not vary with am b .
NRQCD operators relevant for the lattice calculation. ζ ij are the coefficients of the one-loop terms obtained from the renormalisation of the vertex diagram with J (j) at the vertex. Note that Z b , Z m b and ζ ij are all functions of am b and must be evaluated at the value of am b being used in the lattice QCD calculation.
Peeling off the external states we can then write, using A 0 as an example, which determines the C j coefficients of eq. (2). To O(α s ) ij (17) so that, substituting in the results for the η (i) from eq. (4), we have [16] Z q and Z b have logarithmic infrared divergences with aλ as do ζ A0 00 and ζ A0 11 . These cancel against the logarithmic divergences in B  A0 (see eq. (5)) so that C 0,A0 and C 1,A0 are finite. Similarly the linear infrared divergence of B (2) A0 is cancelled by a matching divergence in ζ A0 02 . Note that the explicit factors of aM remaining in the B (j) A0 will now be replaced by am b , which is the same as aM to this order in α s . Combinations of the ζ A0 ij that allow the C j,A0 to be determined are given for the FIG. 3: The z0 and z1 factors for the O(αs) matching of the temporal axial and pseudoscalar NRQCD-HISQ currents to continuum QCD (eqs. (19) and (23)) plotted against the bare lattice b-quark mass.
NRQCD-HISQ case in [10]. The calculation is done for the standard v 4 -accurate NRQCD action [19,20], but the results are also correct for the α s v 4 -improved NRQCD that we will use here [21], because the impact of the α s v 4 improvement terms will only appear in the matching at α 2 s .
Here we recast the expansion of the QCD currents into a more natural combination of lattice QCD currents as The values for z A0 0 , z A0 1 and z A0 2 were given in [9] and are reproduced here in Table II. The values are the same for the temporal vector current from the chiral symmetry of the HISQ action.
To perform the equivalent calculation for the pseudoscalar current we note that the Ω A0 . Hence we do not need to perform a new calculation in lattice QCD perturbation theory. We simply need to reconstruct the mixing matrix for the pseudoscalar case from that of the temporal axial vector. We can then write and find Note the overall minus sign for C 1,P as well as the fact that all of the ζ ij factors with either i or j equal to 1 now come in with opposite sign. These factors are all finite, so the C j,P are still manifestly infrared finite. Table I gives results for the finite ζ factors, ζ A0 10 , ζ A0 01 and ζ A0 12 , as well as Z m b that allow us to determine the C j,P from the C j,A0 for a variety of values of the heavy quark mass in lattice units, am b . These correspond to the values of b quark masses used in our lattice NRQCD calculations that will be discussed in Section III.
For the pseudoscalar current case, we multiply both sides of eq. (20) by the heavy quark mass in the MS scheme at the scale µ and then, on the right-hand side, convert these into lattice NRQCD bare quark masses using the relation Values for Z m b at a variety of am b values for lattice NRQCD are given in Table I [10]. Then we have We find The values of appropriate ζ ij and Z m b given in Table I enable us to determine the z P j values from the z A0 j . The z P j values are given in Table III. The values are the same for the scalar current due to the chiral symmetry of the HISQ action. Fig. 2 shows the different contributions to z A0 0 and z P 0 as a function of am b . This includes the finite pieces of each of the terms in eq. (21) that vary with am b (thereby excluding Z q ). The different contributions are all of moderate size and show mild dependence on am b over the range of am b values that we use. In Fig. 3 we plot z 0 and z 1 for the temporal axial vector and pseudoscalar cases. The magnitudes of z A0 0 and z A0 1 are both very small and both have very little dependence on am b , a fact previously remarked on in [9]. z P 0 and z P 1 have larger magnitude and somewhat more dependence on am b . However both are still smaller than 1 across the range of am b values we use. Note that using J (0) alone in NRQCD to approximate either the temporal axial vector or pseudoscalar currents gives a larger renormalisation factor at O(α s ). This is because the coefficient that represents the 'mixing-down' of J (1) into J (0) , ζ 10 , reduces the size of the one-loop renormalisation of the combined current in both cases. That J (0) + J (1) is much closer to the continuum current than J (0) will be demonstrated in an order-by-order comparison of results in Section III.
The coefficients z A0 2 and z P 2 have stronger am b dependence dominated by that in the mixing coefficient ζ 02 . This grows linearly with am b at large values of am b so that, as am b → ∞, the contribution of J (2) becomes an α s aΛ QCD correction term. In Fig. 4 we show ζ 02 /am b up to large values of am b (much above those that we use in practice) where this behaviour becomes clear. At the am b values that we use the α s aΛ QCD and α s Λ QCD /m b behaviour is intertwined. Values of z A0 2 and z P 2 with the linear am b term removed are shown in Fig. 3.
Lattice QCD results can be combined with eqs. (19) and (23) to determine the hadronic matrix elements of the temporal axial vector/vector and pseudoscalar/scalar currents up to systematic uncertainties coming from missing higher order radiative and relativistic corrections (which will differ between the currents). In the Section III we will compare results for hadronic decay constants obtained using temporal axial or pseudoscalar currents and form factors from temporal vector and scalar currents. The extent to which they agree is a test of our systematic uncertainties.

A. Lattice configurations and simulation parameters
The gluon field configurations used here are listed in Table IV. They are 'second-generation' MILC configurations [22,23] using a gluon action fully corrected through α s a 2 [24] and HISQ quarks [8] with u, d, s and c (n f = 2 + 1 + 1) flavors in the sea. They include multiple values of the lattice spacing and multiple values of the u/d (taken to be degenerate) sea quark mass varying from one fifth of the s quark mass down to the physical value. On these gluon field configurations the B and B s decay constants were calculated in [9] us-  [22,23] used here with their (HISQ) sea quark masses, m l (=(mu + m d )/2), ms and mc in lattice units. β = 10/g 2 is the QCD gauge coupling and the lattice spacing, a, is determined using the Υ(2S − 1S) splitting [9,21]. The lattice size is L 3 s × Lt. Each ensemble contains around 1000 configurations and we take 16 time sources per configuration to increase statistics. ing a radiatively improved (through α s v 4 b ) NRQCD action for the b quark [21,25], the HISQ action for the lighter quark and an NRQCD-HISQ temporal axial current matched to continuum QCD following the process described in Section II. Here we will compare results using the pseudoscalar current matched to the same level of accuracy. In a similar way, the systematic uncertainties in the semileptonic form factor for B → π obtained from (the traditional method of) using a vector current can be tested by employing a scalar current. Since we are largely re-using results from earlier papers [9,30] we do not repeat technical details, for example on the NRQCD Hamiltonian, but refer the reader to those papers for more detail.

B. B and Bs meson decay constants
Using the PCAC relation of continuum QCD we can determine the B meson decay constant, f B , from the matrix element of the temporal axial current between the vacuum and a B meson (at rest) as (25) or from the product of the pseudoscalar density and the quark masses as Here M B is the B meson mass. In [9] the temporal axial current relationship of eq. (25) was used. A 0 was constructed from the leading and next-to-leading NRQCD-HISQ currents in a non-relativistic expansion and was matched to continuum QCD according to eq. (19). This involves writing A 0 in terms of the lattice currents, J A0,lat , J A0,lat and J A0,lat . In [9] the matrix elements of each current between the vacuum and a B meson are determined in lattice QCD, so that the matrix element of A 0 in eq. (25) can be obtained to the specified level of accuracy (the matrix elements for J (1) A0,lat and J (2) A0,lat are the same for a meson at rest).
In eq. (23) we give an expansion to the same order for the combination of quark mass and pseudoscalar density, m b P , in terms of the same NRQCD-HISQ currents multiplied by the bare NRQCD quark mass. Because the matrix elements for each of the lattice NRQCD-HISQ currents are given in [9] we can reconstruct the matrix element of m b P required on the left-hand side of eq. (26) and so determine the decay constant in a different way. This decay constant should agree with that determined from the temporal axial current up to the uncertainties quoted. These are dominated by systematic errors from missing higher order matching terms and relativistic current corrections [9]. Figure 5 shows how well this process works order-byorder as relativistic current corrections and α s matching terms are added in. The plot shows the ratio of the decay constant obtained using the temporal axial current to that using the pseudoscalar density. We use the results from [9], which gives matrix elements for each contribution to the current on each of the ensembles in Table IV. For each ensemble the bare NRQCD quark mass am b is tuned to that of the b-quark using the spin-average of Υ and η b masses and the u/d quark mass, am l , is given the value used for the light quark mass in the sea. The s quark mass is tuned using a fictitious ss pseudoscalar meson whose properties are well-determined in lattice QCD [26]. Values for the π, K and η s meson masses made from these light quarks are given in [26,27]. We use α s in the V-scheme at a scale 2/a in the operator matching as in [9].
In the determination of f B from the pseudoscalar density in eq. (26) there is a factor of (1 + m l /m b ) on the left-hand side. We neglect this for the u/d quark because, at the physical point, m l /m b = 1/(52.55 × 27.4) [28]. This is negligible compared to the other uncertainties. For the additional factor of m B on the right-hand side, which must be removed to take a ratio of the two different decay constants, we use the average of the charged and neutral experimental B meson masses [1]. It has already been demonstrated that the lattice QCD result for the B meson mass, using NRQCD for the b quark, agrees  ) and (23). The dashed lines show the relative uncertainty on fB values quoted in [9].
with experiment at the physical value of the u/d quark mass [27].
In the upper plot of Figure 5 (for the B meson) the lowest (zeroth) order result includes only the J (0) A0,lat current at tree-level, whose matrix element cancels in the ratio, and so the result is simply m b /M B . Not surprisingly, substantial differences are seen between the results of O(20%), being the size of the binding energy of a B meson. A significant improvement is seen on including α s radiative corrections to the normalisation of J  Figure 5 is close to the solid line at 1.0, which indicates the same result is obtained for the decay constant from both the temporal axial current and pseudoscalar density. More importantly the differences from the value 1.0 of this ratio lie within the dashed lines that correspond to the 2.2% relative uncertainty quoted for f B from the temporal axial current in [9]. This uncertainty was dominated by an estimate of the uncertainties expected from missing higher order relativistic current corrections and α 2 s matching errors. Since the results from the temporal axial current and pseudoscalar density will have different uncertainties from the missing α 2 s matching, the result of Figure 5 is a demonstration that the estimates of these uncertainties are realistic.
The lower plot of Figure 5 repeats this exercise for the decay constant of the B s meson, again using results from [9]. In this case, using eq. (26), we do not neglect the s quark mass on the left-hand side. Instead we write and take m s /m b = 1.0/52.55(55) [28]. We use the experimental value of the B s meson mass to determine a ratio of the decay constants from temporal axial current and pseudoscalar density. The lower plot of Figure 5 has identical features to that of upper plot. This is not surprising because the relative effect of the matrix elements of the relativistic current corrections is the same for u/d and s quarks as can be seen from the results in [9]. We can complete the analysis by fitting the decay constant results for f B and f Bs , obtained from the pseudoscalar density, as a function of u/d quark mass and lattice spacing, to extract physical results for comparison to the final answers obtained using the temporal axial current. In the temporal axial current case, because of a normalisation factor from the meson states, the hadronic quantity naturally obtained from the lattice QCD calcu- A0,lat and J (1) A0,lat are tabulated in [9]. The equivalent hadronic quantity corresponding to matrix elements of m b P is f B (M B ) 3/2 . The values we need for m b P matrix elements are obtained by combining the appropriate Φ values for J (0) and J (1) as in eq. (23), including multiplcation by am b and then by an additional power of 1/a to convert to GeV units. In this case, the additional multiplication by 1/a will slightly increase the uncertainty in the values we are fitting because of the uncertainty in the determination of the lattice spacing.
We  Table IV, obtained from the pseudoscalar current and plotted against the light quark mass in units of the strange quark mass (given as M 2 π /M 2 ηs ). The grey bands show the results of the fit described in the text. Errors on the data points include statistics/fitting only; the grey band includes the full error from the fit to lattice spacing and quark mass effects along with the perturbative matching uncertainty. obtained from the pseudoscalar current. The data-points are as shown in Figure 6 and the grey band is the result of the fit described in the text, including uncertainties from lattice spacing and quark mass effects along with uncertainties from higher-order relatvistic corrections to the current. ure 6, as well as plotting our fits to the light quark mass and lattice spacing dependence that enables us to extract a physical result. Following [9] we use a fit form Here d 1 and d 2 allow for discretisation effects and are also δm b dependent (suppressed for clarity). b 1 allows for dependence on the light quark mass, including the chiral logarithm, l(M 2 π ). For the B s case the chiral logarithm term is not present and b 1,l → b 1,s . e 1 allows for α 2 s corrections from only matching to one-loop in perturbation theory, while e 2,3 allow for the fact that the higher order matching coefficients can in principle have am b dependence. These priors are given identical values as in [9], except for e 1 = 0.0(3) as the pseudoscalar matching coefficients are slightly larger than their temporal axial-vector counterparts. Extrapolating to the physical point in the absence of electromagnetism, i.e. M π = M π0 , where m l = (m u + m d )/2, we find f B (M B ) 3/2 = 2.37(7) GeV Our complete error budget is given in Table V, with a breakdown that follows [9]. Errors arising from statistics, the lattice spacing, operator matching and chiral parameters are estimated directly from the fit. The remaining source of systematic error in the decay constants comes from missing higher order relativistic corrections to the current. As discussed in [9], for the heavy-light system under consideration the higher order relativistic corrections will be of the size (Λ QCD /am b ) 2 0.01, which we take for this component of the error.
We can convert the above results into values of the decay constants using the PDG masses [29] for M B l = (M B0 + M B ± )/2 = 5.27963(15) GeV and M Bs = 5.36689 (19). Our final results for the decay constants obtained from the pseudoscalar current are f B = 0.196(6) GeV, f Bs = 0.236(7) GeV and f Bs /f B = 1.207 (7).
As the same matrix elements are used as input when determining the decay constant from temporal axialvector and pseudoscalar currents, we must include correlations when performing an average of the results obtained here and in [9]. As the different values for f B (s) are slightly outside the 1σ error, we scale the error of the weighted averaged value by FIG. 8: The ratio of the scalar form factor at zero recoil for the B to π decay obtained from using the temporal vector current to that from using the scalar current. Results are from lattice QCD calculations of the matrix element of a non-relativistic expansion of the appropriate current operator between a π meson and a B meson. Results from each of the ensembles of Table IV are shown plotted against the square of the lattice spacing. Symbols are as in Figure 5. The dashed line gives the relative error on f0(q 2 max ) quoted in [30].

C. Scalar form factor for B → π decay
Another process where accurate determination of the matrix elements of heavy-light currents is required is for the weak semileptonic decays of B mesons to light mesons. The archetypal process here is B → π ν. The hadronic parameters needed to determine the rate for this decay are known as form factors and they are now functions of q 2 , the squared 4-momentum transfer between the initial and final meson. For B → π ν there are two form factors which we will denote f + and f 0 , but the experimental rate is only sensitive to f + if the final state lepton is light. The form factors are related to the matrix elements of vector and scalar currents by and There is a kinematic constraint that f + (0) = f 0 (0). At the zero recoil, maximum q 2 , point we can compare the matrix elements of the temporal vector and scalar currents directly. At that point, where q 0 = M B − M π , and FIG. 9: The ratio of RBπ values (see text for definition) obtained using a combination of temporal vector and temporal axial vector currents to that from a combination of scalar and pseudoscalar currents. Results from each of the ensembles of Table IV are shown plotted against the square of the lattice spacing. Symbols are as in Figure 5. The dashed line gives the relative error on RBπ from using the temporal vector current for f0 and temporal axial current for fB quoted in [30].
For currents made of NRQCD b quarks combined with HISQ light quarks, the chiral symmetry of the HISQ action guarantees that the nonrelativistic expansion of the temporal vector current has the same form as for the temporal axial vector current given in eq. (19). Likewise m b S has the same expansion as m b P given in eq. (23). As for the case of the decay constant discussed in Section III B, the currents that appear in each order of the nonrelativistic expansion are the same for S and V 0 . Thus we can construct the matrix element of m b S given the lattice matrix elements of the different current contributions for V 0 . These are given for the zero recoil situation in [30] for the same 2+1+1 gluon field configurations as used for the decay constants in Section III B, and listed in Table IV. Figure 8 shows the ratio of the scalar form factor at zero recoil determined from eqs. (31) and (32) using successively more accurate representations of the NRQCD-HISQ temporal vector and scalar currents from eqs. (19) and (23). The matrix elements for the individual lattice current pieces are calculated in [30], noting that the matrix element of J (2) V 0 is equal to that of J (1) V 0 at zero recoil. In the additional mass factors on the left-hand side of eq. (32) we ignore m l compared to m b as it is less than a 0.5% effect across our range of light quark mass values. For M B we average the charged and neutral B meson masses, as in Section III B and for M π we use the values appropriate to these ensembles given in [9,27]. Figure 8 shows a very similar picture to that of Figure 5 with the ratio of the two results becoming closer to 1.0 as non-relativistic current corrections are included and radiative corrections to them added in. With the most accurate matching that we have, the ratio of the results for the scalar form factor differs from one by less than the relative uncertainty on f 0 (q 2 max ) of 3% quoted in [30]. In [30] the ratio of f 0 (q 2 max ) to the decay constant ratio f B /f π was calculated to see if this ratio became one in the massless π meson limit, as expected from soft pion theorems [31][32][33][34]. This was indeed found to be the case, resolving a long-standing issue in the literature. The quantity calculated in [30] was which becomes f 0 (q 2 max )f π /f B as M π → 0. The piece of this ratio that involves b quarks is f 0 (q 2 max )/f B and this was determined for NRQCD b quarks from the ratio of the matrix element between π and B of the temporal vector current divided by the the matrix element between the vacuum and B of the temporal axial vector current. Because of the chiral symmetry of HISQ quarks the matching of these two NRQCD-HISQ currents to their continuum counterparts is the same. Then the overall renormalisation factor (1 + α s z 0 + . . .) (see eq. (19)) cancels between them and the renormalisation uncertainty from missing α 2 s and higher order terms is much reduced. The ratio R Bπ can then be determined to high accuracy.
Here we can also calculate R Bπ , using the scalar current for f 0 (q 2 max ) and the pseudoscalar current for f B . Again the overall renormalisation factor between the two will cancel (see eq. (23)). Figure 9 shows the ratio of R Bπ calculated in these two different ways as, once again, successively more accurate representations of the b-light currents are used. Now, because of cancellation of the overall renormalisation factors, there is no difference between the zeroth order result and the O(α s ) result. Once α s Λ QCD /m b corrections are included the ratio of R Bπ values is very close to one and well within the uncertainty of 2% on R Bπ quoted in [30].

IV. CONCLUSIONS
Here we have given the matching calculation that enables matrix elements of heavy-light scalar and pseudoscalar currents accurate through O(α s Λ QCD /m b ) to be determined in lattice QCD using NRQCD b quarks and HISQ light quarks. This expands the range of methods we can apply to B physics using NRQCD and the tests we can do of our systematic error budget.
In Section III B we determined the B and B s meson decay constants using the pseudoscalar current and then compared to the previously determined values obtained from using the temporal axial current [9]. Since there is no PCAC relation connecting these two currents in lattice NRQCD, they do not have to give the same answer. The way in which the nonrelativistic approximation to the continuum current is built up (albeit from the same ingredients) and the way in which it is renormalised are both different in the two cases. The systematic uncertainties from missing higher order terms will then also be different. We see in Figure 5 the results coming closer together as we include higher orders in the nonrelativistic expansion and matching on both sides. The final results, at the best accuracy that we can currently achieve, agree to within the expected remaining systematic uncertainty. This uncertainty is dominated by unknown α 2 s terms in the overall renormalisation (multiplying the leading current). The agreement is confirmation that the error budget is a reasonable one. Our results for the B and B s decay constants from the pseudoscalar current are: f Bs = 0.236 (7) GeV f Bs /f B = 1.207 (7) .
Our new results for the B and B s decay constants (and their ratio) can be considered as independent values from those obtained using the temporal axial current in [9]. They use the same raw lattice data in the form of matrix elements for the current components and so their statistical uncertainties are correlated. Their systematic uncertainties are not the same, however, since they come largely from unknown, and different, relativistic and α 2 s matching corrections. We can then perform a weighted average of the results arising from the two methodologies, including the statistical correlations, to obtain: These results have very similar uncertainties to those in [9] but do contain more information. Figure 10 gives a summary of lattice QCD results for f B , f Bs and their ratio. It includes results from a variety of light and heavy quark formalisms for calculations that have included at least 2 flavours of quarks in the sea. The most realistic version of QCD corresponds to the results in the top box, including the values we give here, where u, d, s and c sea quarks are incorporated. Our results have the additional advantage of including physical values for the u/d sea quarks, taking m u = m d . The results for f B plotted in Figure 10 correspond to a B meson made with a light quark with the u/d average mass. The grey bands show average values from [1], where there is also discussion of the effects of isospin-breaking. The main message from Figure 10 is that of good agreement between the different lattice QCD results, which is another good test of systematic uncertainties.
A similarly encouraging picture was given of the comparison of temporal vector and scalar current results for the scalar form factor for B → π decay in Section III C. The calculations compared were done at zero recoil where the π is at rest in the rest frame of the B. Here we discuss briefly the potential uses of our new method to determine the vector form factor f + away from the zero recoil point, where connection to experimental decay rates can be made for the determination of Cabibbo-Kobayashi-Maskawa matrix element |V ub |.
Power-counting in powers of the inverse heavy quark mass must be modified for current matrix elements away from the zero recoil point as the momentum of the light meson in the final state increases. Sub-leading currents that include a spatial derivative on the light quark field will have matrix elements that grow as |p |/m b and these can become relatively large compared to the leading order current if p > Λ QCD . The issue is discussed for the NRQCD-asqtad case in [42] where ratios of the subleading matrix elements to the leading matrix elements between the B meson and a heavy pion of the spatial vector current are shown. The matrix elements for currents denoted J (2) k and J (4) k grow as a proportion of the leading order (J (0) k ) matrix element as the pion momentum is increased. This is also true, although not shown there, for the matrix elements of the temporal vector current J (2) 0 . These three currents are analogous to the current J (2) A0,lat (eq. (13)) considered here, in having a derivative on the light quark field; for the spatial vector there are two such currents with different γ matrix structures. Note that the matrix elements for the sub-leading currents that contain a derivative on the heavy-quark field show much more benign behaviour, as might be expected.
The sub-leading currents with derivatives on the light quark field do not appear at tree-level in the expansion of the continuum heavy-light current (see eq. (19) and so are suppressed by powers of α s . In the NRQCD-asqtad case the α s coefficients of these sub-leading currents were calculated and turned out to be small for the largest contribution, from J (4) k [42]. For the NRQCD-HISQ case these coefficients are only known for the temporal vector current (see Table II [10]). As we move away from zero-recoil in B → π decay, systematic uncertainties from these sub-leading currents will grow if they are not included in our nonrelativistic expansion of the continuum current. It is therefore important to work with NRQCD-HISQ currents that do include the subleading currents with derivatives on the light quark field so that accuracy can be maintained as far from zero recoil as possible.
Here we have provided a way of doing this by using the temporal vector and scalar currents (and eqs. (29) and (30)), as used for example in calculations with purely HISQ quarks [13,14]. As we have shown, both of these continuum currents can be written as a nonrelativistic expansion in NRQCD-HISQ currents that includes terms that will become O(α s |p |/m b ) away from the zero recoil point. Uncertainties are then O(α 2 s |p |/m b ) and O(α s |p | 2 /m 2 b ). This approach will be used in NRQCD-HISQ work on the second-generation 2+1+1 HISQ configurations, extending [30] away from zero-recoil.
It should also be noted that the expansion for the pseudoscalar heavy-light current will allow more form factors to be separated out in the analysis of B meson decays to light vectors [43,44]. Processes such as B s → φ + − and B → K * + − provide key opportunities for stringent tests of the Standard Model [45,46] and will need increasingly accurate lattice QCD results for comparison. authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U. S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publica-tion, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.