HQET at order 1/m: III. Decay constants in the quenched approximation

We report on the computation of the $B_s$ meson decay constant in Heavy Quark Effective Theory on the lattice. The next to leading order corrections in the HQET expansion are included non-perturbatively. We estimate higher order contributions to be very small. The results are extrapolated to the continuum limit, the main systematic error affecting the computation is therefore the quenched approximation used here. The Generalized Eigenvalue Problem and the use of all-to-all propagators are important technical ingredients of our approach that allow to keep statistical and systematic errors under control. We also report on the decay constant $f_{B'_s}$ of the first radially excited state in the $B_s$ sector, computed in the static limit.


Introduction
Flavour physics is becoming a precision field. B-physics measurements may produce stringent tests of the Standard Model (SM) and consequently reveal possible effects coming from New Physics. They are complementary to direct searches and they provide constraints on the flavour structure of any possible extension of the Standard Model. At the moment the significance of such tests is limited by the uncertainties on the theoretical side [1]. A typical example is the process B s → µ + µ − . The SM prediction for the branching ratio is O(10 −9 ) [2][3][4] and the best experimental upper bound (from D0) is 4.2 × 10 −8 @ 90% CL [5]. The decay is very sensitive to an extended Higgs sector and may be strongly enhanced in various extensions of the Standard Model (e.g. the supersymmetric model discussed in [6]). LHCb has a potential to measure a branching ratio as small as 9 × 10 −9 at 3 σ with 0.1 fb −1 of data [7]. The hadronic parameter entering the SM prediction is the B s meson decay constant f Bs , which is known from the lattice with an uncertainty of about 15% [8,9].
More precise lattice computations are needed to make progress, however heavy quarks on the lattice are difficult due to O((am b ) n ) discretization errors, where a is the lattice spacing. A description of heavy-light systems which is suitable for lattice QCD simulations is given by Heavy Quark Effective Theory (HQET) [10,11] with nonperturbatively determined parameters [12].
In this paper we report on a quenched computation of f Bs performed entirely in HQET including 1/m b corrections non-perturbatively. The plan of the paper is the following. In section 2 we restate the strategy that we have used and already explained in [13], with particular emphasis on the use of the GEVP variational method [14]. In section 3 we give the numerical values of f Bs and f B s obtained at the 3 lattice spacings that we have considered and discuss the extrapolation to the continuum limit. We briefly conclude in section 4.

Non-perturbative HQET
We aim at computing the decay constant f Bs , defined in QCD as with the normalization of states B s (p)|B s (p ) = 2E(p)δ 3 (p−p ), from matrix elements defined in HQET. To this end we need to match the HQET Lagrangian and the currents to their QCD counterparts. To order 1/m b , the HQET Lagrangian reads where ψ h satisfies 1+γ 0 2 ψ h = ψ h , and ω kin and ω spin are matching parameters whose tree-level values are ω kin = ω spin = 1/(2m b ), and δm is a counter-term that absorbs the power-divergences of the static quark self energy.
Again to order 1/m b , the time-component of the QCD axial current A QCD where all derivatives are symmetrized The renormalization constant Z HQET A depends on the ratio m s /m b . This is a small effect, which is further reduced by a factor of the coupling constant α(m b ). We will ignore this dependence and use the value of Z HQET A determined with a massless light quark [12]. Note in addition that the operator A (2) 0 does not contribute to correlation functions and matrix elements at zero spatial momentum, such as those we are interested in here.
At the static order the Lagrangian is automatically O(a) improved, therefore the current and its on-shell matrix elements are O(a) improved if one sets c As long as we restrict our studies to the decay constants only, to fully specify HQET the parameters δm, ω kin , ω spin , Z HQET A , and c (1) A must be determined by matching the effective theory to QCD. Using the Schrödinger functional, our collaboration has performed a fully non-perturbative determination of the parameters of HQET [12]. Here we employ the same discretization of QCD and HQET and in particular use the determined values for the parameters of the effective theory.

The Generalized Eigenvalue Problem
We follow here the application of the GEVP [16,17] described in [14]. For the sake of completeness we recall the basic ingredients of the method. The matrix of Euclidean space correlation functions between the zero-momentum projection a 3 provides the basis for the Generalized Eigenvalue Problem (GEVP) An effective creation operator for the n th state can be defined bŷ Defining the vector of correlators of a composite field P (which does not have to be among theÕ i ) C P,i (t) = P (t)Õ * i (0) , i = 1, . . . , N , (2.14) the effective matrix elements approximate the matrix elements of the corresponding operatorP as The definition of R n in eq. (2.12) is slightly different from the one in [14] and has the advantage of being defined at all (and not only even) values of t, thus giving better statistical precision for the final result while preserving the same control (2.17) over the contamination from excited states as proven in [14]. After expanding the correlators to first order in ω ∼ 1/m b we consider the GEVP in perturbation theory in 1/m b and find .
Thus, in order to obtain the effective matrix elements, the GEVP has to be solved for the static correlation functions only With these definitions, and by organizing the 1/m b expansion in the way we discussed in the previous section, the decay constant of a pseudoscalar B s meson (n = 1) or of radial excitations (n > 1) computed in the static approximation and in HQET (i.e. including terms of order 1/m b ), respectively, read where p stat n , p kin n , p spin n and p A (1) n are the plateau values of the corresponding effective matrix elements (see [14] where however p A (1) n is called p δA n ). For the improvement term proportional to b stat A am q we use the 1-loop estimates of the coefficient b stat A from [18].
In the formulae am q is the bare subtracted strange quark mass 1 2 1 κs − 1 κc , with κ c the critical value of the hopping parameter defined through the vanishing of the quark mass derived from the axial Ward identity.
In order to consistently truncate the expansion at order 1/m b , it is convenient to take the logarithm of (2.23) and expand the logarithms (rather than expanding directly the product of the factors from the correlation function times its renormalization constant). Parameters of the simulations: inverse coupling β, approximate scale parameter r 0 in lattice units [26], spacetime volume, hopping parameter corresponding to the strange quark mass [27], critical hopping parameter [21], and numbers of low-lying eigenmodes and stochastic noises used.
3 Numerical results

Simulation parameters
We are now ready to present the result of our numerical simulations to extract f Bs . The parameters of the simulations are given in Table 1. Each ensemble contains 100 quenched configurations. The heavy quark is described by the HYP1 and HYP2 static actions [19][20][21] while the valence strange quark is described by the non-perturbatively O(a)-improved Wilson action [22,23]. Our lattices are L 3 × T with L ≈ 1.5 fm, T = 2L, and periodic boundary conditions are applied in all directions. We use all-to-all propagators based on the Dublin method [24], but with even-odd preconditioning and N L approximate (instead of exact) low modes; for details of our method the reader is referred to [25]. No low modes have been computed for β = 6.4956 because the numerical cost would have been too high with respect to the gain in statistical precision; instead, we have improved the statistics by using N η = 4 stochastic noises, twice the number of noise sources used at the other lattice spacings.

Bare matrix elements
In at each of our three lattice spacings for both the HYP1 and the HYP2 static quark action.
The interpolating fields are constructed using quark bilinears built from the static quark field ψ h (x) and different levels of Gaussian smearing [28] for the light quark field with APE smeared links [29,30] in the Laplacian with exactly the same parameters as in [25].   For these bilinears, we compute the following correlators: where the O(1/m b ) fields and A 0 have been defined in eq. (2.4) and eq. (2.6). We have followed the procedure explained in detail in [25] to choose the time ranges over which we fit the various plateaux. Some examples of the plateaux found are shown in Figure 1; it can be seen that without some knowledge of the analytical form of the leading corrections it would often be difficult to tell whether a reliable plateau has been found.
We first fit the matrix elements to the expected form (where x ∈ {kin, spin}) using the energy levels extracted by the procedure described in [25] as input parameters. Then, in a second step, we form plateau averages starting from t 0 = t 0,min at each value of N and ∆t = t − t 0 , and take as our final estimate that plateau for which the sum σ tot = σ stat + σ sys of the statistical error σ stat of the plateau average and the maximum systematic error σ sys = π(t, t 0,min ) becomes minimal, subject to the constraint that σ sys < 1 3 σ stat . We impose the latter constraint in order to ensure that the total error is dominated by statistical errors. The extracted bare matrix elements are given in table 2, quoting not only the final plateau average, but also the result of the fit, which generally agrees rather well with the final result.

Continuum limit
From the bare matrix elements and the parameters of HQET, determined in [12] with the HYP1 and HYP2 static actions, we form dimensionless quantities in which the divergences cancel exactly to order O(1/m b ), and use them to compute f x For a comparison to the static limit and previous work, we also consider where Z stat A,RGI is the renormalization factor of the Renormalization Group Invariant static-light axial current, as defined in [31]. In contrast to Z stat A,RGI , the HQET parameter Z stat A in eq. (3.5) has been determined by a non-perturbative matching at finite mass. The correspondence is (3.9) in terms of the conversion function C PS introduced in [31] and now known up to threeloops [32][33][34]. For Z stat A,RGI we use the non-perturbative value from [21]. Since both of the static actions used are discretizations of the same continuum theory, we perform a combined continuum limit by fitting a function of the form (k = 1, 2 for HYP1, HYP2 actions) incurred from using the one-loop value of c stat A , we compute the continuum limit also for c stat A = 0 using a quadratic extrapolation and compare the result to the continuum limit obtained using the oneloop value. As can be seen from table 3, the influence of c stat A is negligible at this level.
Statistical errors are computed by a jackknife analysis that also includes the correlation among HQET parameters. We find that the results obtained from taking the continuum limit ofΦ x n and using it to compute Φ x n , and from taking the continuum limit of Φ x n directly agree well within the errors. In Figures 2 and 3, we show the continuum extrapolations of some relevant quantities. It is easily seen on those plots that the combination of HYP1 and HYP2 results is legitimate because they point to the same continuum limit within the errors. Numerically we finally get:    with the first three continuum-limit results one has to consider the combination C PS × Φ RGI 1 = 2.20(4) (1-loop c stat A ). Here the conversion function C PS [31,32] has been computed with the three-loop anomalous dimension from [33,34].  Figure 4) and the non-perturbative result in [35] for the decay constant around the charm quark mass (blue triangles in Figure 4). The static limit of this quantity is r 3/2 0 Φ RGI 1 , which we have also non-perturbatively computed here (purple square in Figure 4). As explained above we rely on perturbation theory only for the evaluation of the conversion function C PS (M/Λ MS ) relating the RGI matrix elements in static HQET with their counterpart in QCD defined at a given heavy quark mass [31,32]. Thus, dividing by C PS (M b /Λ MS ) compensates for the well-known logarithmic scaling of the decay constant with the heavy-quark mass [36,37]. One can see that our result is falling rather well on the straight line expected from heavy quark scaling, indicating that the neglected O(1/m n b )| n≥2 corrections are small. We note, however, that this comparison and conclusion rely on the perturbative evaluation of C PS , and that the associated α s (m) 3 errors are very difficult to estimate.
Our use of the GEVP method also allows us to extract some information on the matrix element for n = 2, i.e. of the first excited state of the B s system, for which we quadratically extrapolated to the continuum limit. The unimproved version of this quantity (i.e. p stat 2 /p stat 1 ), quadratically extrapolated to the continuum, gives the same result of 1.24 (7). We have obtained the same qualitative result as [38] concerning this ratio: it is noticeably larger than 1, in good qualitative agreement with predictions from quark models that become Lorentz covariant in the heavy quark limit [39] and relativistic quasi-potential quark models [40,41], while other models predict a value less than 1 for this quantity [42].

Conclusion
In this paper we have reported on the computation of the B s meson decay constant by using lattice simulations in quenched HQET. Including 1/m b corrections introduces power divergences ∼ 1/(am b ) which have to be subtracted non-perturbatively. These non-perturbative subtractions have here been carried out successfully for the first time in lattice gauge theory computations. The necessary couplings of the effective theory had been determined non-perturbatively by matching it to QCD [12].
Our strategy had already been developed earlier [13] but its implementation revealed relatively large statistical errors in the matrix elements of the 1/m b operators (not due to the computation of the non-perturbative parameters of the theory). This shortcoming has now been cured by exploiting (i) a method based on solving a GEVP to reduce the systematic errors on bare matrix elements coming from the contribution of excited states to correlation functions and (ii) all-to-all propagators to improve the statistical precision. For example, at the finest lattice resolution considered, we have obtained a result for the bare static decay constant (HYP2 action), which is three times more precise than the result in [35] at β = 6.45 where ten times more configurations were analyzed in the Schrödinger Functional setup.
We used three lattice spacings to extrapolate to the continuum limit. With r 0 = 0. corrections are very small. In addition, we have shown that the GEVP method is useful for studying phenomenologically interesting quantities involving radial excitations of mesonic states, such as the ratio f B s /f Bs . In this respect we confirm a recent lattice calculation [38] finding that f B s /f Bs > 1 at least in the static approximation.
We intend to apply the approach described in this paper to the computation of f Bs and f B on dynamical N f = 2 configurations in the near future. Note that the problem posed for charm physics on the lattice by the rapid slowing-down of the topological modes of the gauge fields with decreasing lattice spacing [43][44][45] is less relevant in this case, since in HQET we can afford to work with coarser lattices.