Decay constants of B-mesons from non-perturbative HQET with two light dynamical quarks

We present a computation of B-meson decay constants from lattice QCD simulations within the framework of Heavy Quark Effective Theory for the b-quark. The next-to-leading order corrections in the HQET expansion are included non-perturbatively. Based on Nf=2 gauge field ensembles, covering three lattice spacings a (0.08-0.05)fm and pion masses down to 190MeV, a variational method for extracting hadronic matrix elements is used to keep systematic errors under control. In addition we perform a careful autocorrelation analysis in the extrapolation to the continuum and to the physical pion mass limits. Our final results read fB=186(13)MeV, fBs=224(14)MeV and fBs/fB=1.203(65). A comparison with other results in the literature does not reveal a dependence on the number of dynamical quarks, and effects from truncating HQET appear to be negligible.


Introduction
In the ongoing quest for new effects in high-energy particle physics, flavour physics provides information complementary to that from the direct searches performed at ATLAS and CMS. Indeed, low-energy processes and rare events can be sensitive probes of New Physics, in particular when they are mediated by virtual loops, in which non-Standard Model particles can circulate, or when they involve new couplings occurring at tree-level. However, any analysis of experimental data in the quark sector depends on theoretical inputs, such as hadron decay constants, that encode the long-distance dynamics of QCD, which cannot be reliably estimated in perturbation theory.
In this regard, B-physics is an emblematic case. For example, it is crucial to understand the origin of the current discrepancy in the Cabibbo-Kobayashi-Maskawa matrix element V ub measured through the exclusive processes B → τ ν [1,2] and B → π ν [3,4], where the latter makes use of the B → π form factors computed on the lattice. Is it due to an experimental problem, or due to new physics such as the presence of a new, righthanded, tree-level coupling to a charged Higgs boson in the B leptonic decay [5], or due to a severe underestimate of the uncertainty on the decay constant f B governing that decay. For comparison, the recent measurements of B(B s → µ + µ − ) at LHC [6,7] are in excellent agreement with the Standard Model prediction [8,9], where the latter depends on the decay constant f Bs , whose estimate is dominated by lattice results.
The methods that have been used to estimate f B and f Bs include applications of quark models, as discussed in [10][11][12] and references therein, and QCD sum rules in the analysis of two-point B-meson correlators [13][14][15][16]. Several strategies have been proposed to determine f B and f Bs from first principles using lattice field theory, including the extrapolation of simulation results obtained in the region between the charm quark mass m c and a mass ∼ 3m c to the physical b-quark mass m b [17,18], simulations of relativistic b-quarks using an action tuned so as to minimize discretization errors [19,20], and the use of Non-Relativistic QCD [21,22]. We here use Heavy Quark Effective Theory (HQET) [23][24][25][26], regularized on the lattice with the parameters of HQET determined by a non-perturbative matching to QCD [27][28][29]. The virtue of this approach is that perturbative errors are absent and the continuum limit exists. The matching at order O(1/m h ) has been performed in the N f = 2 theory [30]; this is the first step required, for example, in order to compute the b-quark mass [31], which we use in this letter to extract f B , f Bs and f Bs /f B from our simulations. In Section 2, we review the methods employed, before presenting the results in Section 3. Section 4 contains our conclusions.

HQET on the lattice
Heavy Quark Effective Theory regularized on the lattice is a well-defined approach to Bphysics. It is based on an expansion in powers of 1/m h of QCD correlation functions around the limit 1/m h → 0. The continuum limit can be taken order by order in the expansion, since it only requires correlation functions computed in the static theory, which is nonperturbatively renormalizable. Applying the strategy previously discussed in [27,29,32] and employed to measure f Bs in the quenched approximation [33], the HQET action and and We use the labels h to denote the heavy (static) quark field appearing in the HQET Lagrangian, and q = u/d, s for the light and strange quark channels, respectively. The normalization is such that the classical values of the coefficients are ω kin = ω spin = −c (i) A = 1/(2m h ). A bare mass m bare has to be added to the energy levels computed with this Lagrangian in order to obtain the QCD ones. At the classical level it is m h , but in the quantized theory, it has to further compensate a power divergence. The heavy quark field ψ h obeys 1+γ 0 2 ψ h = ψ h , and all spatial derivatives are symmetrized: In QCD the decay constant f Bq is given by 0|qγ 0 γ 5 b|B q (p = 0) = f Bq m Bq , with the relativistic convention B q (p )|B q (p) = 2E q (p)δ(p − p ). We are thus interested in extracting matrix elements from correlation functions defined at zero spatial momentum; the operator A (2) 0 therefore does not enter into our present computations at all. In order to assure the renormalizability of HQET at next-to-leading order, the O(1/m h ) terms in (2.1) are treated in the usual way as operator insertions into static correlation functions, where the suffix "stat" reminds us that expectation values are computed in the static theory. In our strategy we can treat HQET non-perturbatively to leading order (static) or next-to-leading order, including terms of O(1/m h ). The corresponding sets of parameters, ω stat ≡ m stat bare , Z stat A or ω ≡ m bare , ω kin , ω spin , Z HQET A , c (1) A absorb power and logarithmic divergences of the effective theory regularized on the lattice. For technical reasons they have been determined in [30] for a series of heavy quark masses parameterized in terms of the renormalization group invariant (RGI) heavy quark mass z ≡ M L 1 . 1 After our recent determination of the RGI b-quark mass z b ≡ M b L 1 [31], we now choose a quadratic polynomial to interpolate ω(z) -computed at z = 11, 13, 15 -to the physical point z b = 13.25. Similarly, we interpolate ω stat (z) to z stat b = 13.24 at the static order. As expected from [30], all individual interpolations of the ω (stat) (z) parameters to z (stat) b are smooth and do not deviate much from the closest point at z = 13. In the following we will refer to the HQET parameters at the physical b-quark mass only. For completeness they are collected in Table 1 for the three lattice spacings a(β) and two static discretizations (HYP1, HYP2) in use.

Isolating the ground state
In extracting hadronic quantities, it is crucial to have good control over the contributions of excited states to the correlators. Since we are interested in the lowest-lying state in a given channel, we attempt to suppress excited-state contaminations by considering appropriate linear combinations of correlation functions. Specifically, we form the matrices which depend on a basis of operators O i , i = 1, . . . , N . The key step is to find a solution to the generalized eigenvalue problem (GEVP) in the static approximation only depend on the static generalized eigenvalues λ stat n (t, t 0 ), the vectors v stat n (t, t 0 ), and the O(1/m h ) correlators C kin/spin (t) [35], in analogy with perturbation theory in quantum mechanics.
We define the operator basis [36] relativistic quark field. The gauge links in the covariant Laplacian ∆ have been triply APE smeared [37,38] in the spatial directions.
The parameters κ G and R k have been chosen such that they correspond to approximately the same sequence of physical radii at each value of the lattice spacing, see [31] for details. We solve the GEVP for the matrix of correlators in the static limit, eq. (2.11) for N = 3. The resulting eigenvalues and eigenvectors, together with the matrices in eq. (2.10) and the correlators are used to build optimal interpolating fields such that the matrix elements 0|A stat 0 |B (n) and their O(1/m h ) corrections can be extracted, from the correlation functions above, up to contaminations from excited states which are . Notice that for the ground state the energy difference in the exponential correction is of the form E N +1 − E 1 rather than E 2 − E 1 , with N the rank of the correlator matrices in eq. (2.10). This asymptotic convergence holds for t 0 ≥ t/2, as discussed in detail in [35], to which we refer for any unexplained notation. In particular, the symbols , with x=kin/spin, are used in the following to indicate the static and O(1/m h ) contributions to the matrix elements of the axial current as in [33], where one can find expressions for the expected time dependence of the different terms, which read  (15) 215(7) 213 (8) 232 (8) 231 (9) 1.077 (28) 1.086 (25)  F6 0.0562 (9) 203 (8) 201 (8) 228 (7) 228 (7) 1.120 (48) 1.138 (39)  F7 0.0449 (7) 201 (6) 200 (6) 222 (6) 223 (7) 1.103 (26) 1.119 (24)  G8 0.0260 (5) 190 (8) 190 (8) - 0.0940 (24) 222(16) 221(15) - 0.0662 (10) 205(14) 205(15) 229(15) 231(15) 25 as determined in [31]. The last two rows summarize our results of a combined chiral and continuum extrapolation using either the LO or the NLO fit ansatz (3.1) for each individual observable.
We provide in Figure 1 an illustration of typical plateaux for the heavy-light and heavystrange mesons matrix elements p stat and p spin . Those plateaux are chosen following the procedure detailed in [33,39]. The criteria use the results of the GEVP analysis to ensure that in the plateau region the systematic errors due to excited-state contributions are less than a given fraction (typically one third) of the statistical errors. As a consistency check, we have also employed a global fit of the form of eqs. (2.14) to our data. The values of p n obtained from the fit were consistent with the plateau values, albeit with smaller statistical errors. Our errors may therefore be seen as a conservative estimate.

B-meson decay constants at different orders in HQET
The quantities of interest are obtained by combining the lattice parameters of HQET, computed non-perturbatively, and the bare matrix elements evaluated in the static theory. All divergences of the effective quantum field theory are thus properly removed and the continuum limit can be safely taken at a fixed order in the 1/m h expansion.
To obtain f B and f Bs including O(1/m h ) terms in HQET, we compute  and their ratio f stat Bs /f stat B , using static HQET parameters at the physical point ω stat (z = z stat b ), with z stat b = 13.24 as determined in [31]. The last two rows summarize our results for a combined chiral and continuum extrapolation using either the LO or NLO fit ansatz (3.1) for each individual observable.
current. 2 The corresponding expression in static HQET is given by where c stat A is another O(a) improvement coefficient that is needed at the static order. Both b stat A and c stat A have been computed perturbatively in [40]. If treated as an independent observable, the ratio f Bs /f B is easily obtained through eq. (2.15) at next-to-leading order, viz.
f Bs /f B = exp{φ s − φ} m Bs /m B , (2.17) or in the very same way through eq. (2.16) at the static order. Note that in the ratio the leading dependence on the scale setting procedure, which explicitly enters via the lattice spacing appearing as a 3/2 , cancels. Furthermore, terms in eq. (2.15) or eq. (2.16), which do not carry an explicit label i drop out and the term multiplying b stat A becomes independent of the critical hopping parameter κ crit (g 2 0 ). Concerning the parameters and overall statistics of the large-volume simulations used in the present analysis we refer the reader to Table 1 of [31]. The light quark is treated in a unitary setup with m π covering a range from 190 MeV to 440 MeV while the bare strange quark mass is tuned on each CLS ensemble using the kaon decay constant [34]. The lattice spacings are a/fm ∈ {0.048, 0.065, 0.075} for β ∈ {5.5, 5.3, 5.2}, corresponding to the CLS ensemble ids N-O, E-G and A-B respectively. In Table 2 Table 3.

Error analysis and propagation
We follow [41,42] Table 2 and 3), as well as functions of additional input Y , such as the HQET parameters ω i . Also the results of fits to the data are considered as functions of the original data, where the weights in the fits (we always use only the diagonal errors as weights) are precomputed and then kept fixed, i.e., a dependence of f on the weights is not considered.
The error σ f of such a function is then The block-diagonal covariance matrix C Y of the additional input is known: a block [34] for the axial current renormalization factors at the three different β (entering the lattice spacing determination and f π ) and a block [30] for the ω i . The contributions from the individual ensembles e are , The term proportional to τ e exp accounts for the difficult-to-estimate contribution of the tails to the autocorrelation function Γ e f [42]. For τ e exp we insert our previously estimated values  (see e.g. Table 1 of [31]), and W is chosen as the point where Γ e f comes close to zero within about (1-2 ×) its estimated statistical error. The required derivatives ∂f ∂pα are computed numerically [41].
We note that there are many hidden correlations which are all taken into account, e.g., the lattice spacing at one β depends on information from other β through the combined chiral extrapolation in [34]. A straightforward implementation of eq. (2.19) would be cumbersome and numerically expensive. We compute it iteratively instead [43,44].
As explicit example we choose the ensemble with the highest statistics available, e = N6. Figure 2 shows the numerical estimate of the normalized autocorrelation function ρ(t) in terms of the simulation time in molecular dynamic units (MDU). After summing up the autocorrelation function explicitly within a window where it is still rather well determined, the sum up to infinity is determined by modelling it with a single exponential exp(−t/τ exp ) plotted as "tail". On the left hand side of the figure the observable chosen is the lattice spacing (see [34]); the relevant contribution is from the kaon decay constant. On the right hand side we have chosen φ, i.e., essentially the B-meson decay constant in lattice units (see eq. (2.15)). While on the left, where only light-quark physics enters and measurements were taken more frequently, the tail is seen quite well at t < 100 MDU, on the right the autocorrelation function appears to drop to a lower value at short time and in fact is not significant at t ≈ 30 MDU. Still, our somewhat conservative procedure estimates a ≈ 35% contribution to τ int on the left and ≈ 82% on the right.

Continuum and chiral limit extrapolations
We use formulae from Heavy Meson Chiral Perturbation Theory (HMχPT) when applicable [45][46][47]: As in [31] we have parameterized the chiral behaviour through the variable y = m 2 π 8π 2 f 2 π , with y exp representing the physical value at f exp π = 130.4 MeV and m exp π = 134.98 MeV. Since we employ two static discretizations, we also need to account for different cutoff effects, parameterized by D δ , with δ = 1, 2 corresponding to HYPδ. Due to the universality of the continuum limit, the other coefficients do not depend upon δ. After having fixed the HQET parameters at the physical b-quark mass, see Table 1, we treat the B-meson masses as constants, which are taken as m B = 5279.5 MeV and m Bs = 5366.3 MeV from the PDG [48]. Together with all data points that enter the joint chiral and continuum  Table 2.
extrapolation, we list our results from different, independent fit ansätze in Table 2 and 3. For consistency, we decided to treat f Bs as the dependent observable, to be derived from our final results The first, statistical, error as obtained from the NLO HMχPT fit ansatz also includes the discrepancy to the static result, the uncertainty from the HQET parameters and the lattice spacings. We add a second, systematic, error to account for the uncertainty in the chiral extrapolation. It is given by the difference between the quoted value and its counterpart obtained by employing the LO fit ansatz for the chiral extrapolation. While we show only the NLO extrapolation of f B in the left panel of Figure 3, we also add the continuum extrapolated value from the LO fit ansatz. With all correlations taken into account, our estimate for the B s -meson decay constant becomes In the right panel of Figure 3, we contrast this result (filled triangle) with an extrapolation of our f Bs lattice data as if treated as an independent observable, c.f.  In Table 5 we split the statistical error of our observables among different sources. Obviously, the errors from the HQET parameters ω and renormalization factor Z A which enters through the scale setting, largely cancel in the ratio. Although one looses precision, in general, due to the increased variance in HQET observables compared to observables in the light quark sector (such as m π or f π ), one is in the fortunate position that the former couple less to the slow modes of the Monte Carlo chain, and therefore their integrated autocorrelation times are smaller than for "light" quantities.

A quick look at phenomenology
The Flavour Lattice Averaging Group (FLAG) [49] has made a selection of lattice results for f B , f Bs and f Bs /f B with N f = 2, 2 + 1 and 2 + 1 + 1 dynamical quarks [17][18][19][20][21][22]. Only one determination entered the two-flavour average and has been updated [18] since  (22). As a phenomenological application, we can insert our results for f B and f Bs into the formulae describing the branching ratios of B → τ ν τ and B s → µ + µ − transitions: Here Y ≡ Y (x tW , x Ht , α s ) takes into account various electroweak and QCD corrections, parameterized by x tW = m 2 t /m 2 W and x Ht = M 2 H /m 2 t with M H being the Higgs boson mass. Using as inputs the experimental value B(B → τ ν τ ) exp = 1.05(25) × 10 −4 quoted by the PDG [2,48,[50][51][52] and our estimate of f B , we get |V ub | = 4.15 (29) where the errors come from f B and the branching ratio, respectively. The value is roughly 1.5 σ above the exclusive determination from B → π ν. Moreover, using the recent combination of experimental measurements at LHC, namely B(B s → µ + µ − ) = (2.9 ± 0.7) × 10 −9 [6,7,53], together with our determination of f Bs , and all input parameters of (3.5) set as in [8], we obtain The number is in good agreement with the extraction from global fits, which is mostly constrained by B 0 s − B 0 s mixing.

Conclusions
In this paper we have reported on our lattice measurement of the decay constants f B and f Bs performed with two dynamical flavours of O(a) improved Wilson fermions. The b-quark is treated in HQET, with the matching to QCD performed non-perturbatively. This makes the computation entirely non-perturbative, with no reference to continuum renormalized perturbation theory at any point. After an extrapolation to the chiral and continuum limit, we obtain Though it is important to check the dependence of these results on the number of dynamical flavours, and therefore to repeat the computation with a dynamical strange quark, it may still be interesting to compute the ratios f B * /f B and f B * 0 /f B on the N f = 2 ensembles. The first one is often used to check the reliability of sum rules in the B-sector [54]. A lattice measurement at O(1/m b ) requires the matching coefficients that are being computed by the ALPHA Collaboration to extract B → π ν form factors [55]. The second ratio, already in the static limit, can be used to gain some insight into the precision of phenomenological applications of HMχPT, in particular concerning the relevance of the contributions from the J P = {0 + , 1 + } doublet states in chiral loops [56].
The method of the present paper to compute B-meson decay constants has been used previously in the framework of quenched QCD to estimate f Bs without inclusion of virtual quark loops [33]. There, the scale r 0 defined via the static quark potential [57] was employed to express the decay constant in physical units, corresponding to f N f =0 Bs = 216(5) MeV for r 0 = 0.5 fm and f N f =0 Bs = 252(7) MeV for r 0 = 0.45 fm. Given the rather reliable evidence that the true r 0 in physical units lies in between these values (see [58] for a review of the current status), our final result in eq. (4.1) is compatible with the quenched one at the present level of precision. Hence, no significant N f -dependence can be stated.
An interesting piece of information is also contained in the technical Table 5. It shows that the uncertainties in the non-perturbatively determined HQET parameters contribute only at the level of 8% in the static limit and 14% when 1/m b terms are included. Moreover, we find the O(1/m b ) corrections to be very small, 2.5%. This, together with the fact that the computation of the ω i can be much improved with today's machines, gives us confidence that errors can be significantly reduced in the future computation with 2+1 dynamical flavours.