$B$ Meson Decay Constants $f_{B_c}$, $f_{B_s}$ and $f_B$ from QCD Sum Rules

Finite energy QCD sum rules with Legendre polynomial integration kernels are used to determine the heavy meson decay constant $f_{B_c}$, and revisit $f_B$ and $f_{B_s}$. Results exhibit excellent stability in a wide range of values of the integration radius in the complex squared energy plane, and of the order of the Legendre polynomial. Results are $f_{B_c} = 528 \pm 19$ MeV, $f_B = 186 \pm 14$ MeV, and $f_{B_s} = 222 \pm 12$ MeV.


Introduction
The decay constant of a pseudoscalar meson B q consisting of a heavy bquark and a lighter q-quark, with q = u, d, s, c, is defined through the matrix element of the pseudoscalar current where |Ω is the physical vacuum, M Bq is the mass of the B q meson and f Bq the corresponding leptonic decay constant. These decay constants are of great phenomenological interest since they enter as input in non-leptonic B decays, in the hadronic matrix elements of B −B mixing, and in the extraction of CKM matrix elements |V cb |, |V ub | from the leptonic decay widths of B mesons. The so-called hadronic B parameter, which is directly related to the deviation from the vacuum saturation hypothesis, also requires knowledge of these leptonic decay constants, which are of central interest to the ongoing experiments carried out at B factories. Since their determination involves non-perturbative QCD effects, one has to rely on essentially two approaches, QCD sum rules (QCDSR) and lattice QCD (LQCD) simulations. Since the pioneering work of Shifman, Vainshtein and Zakharov [1] the sum rule method has been used successfully to calculate various low energy parameters in QCD [2]. Specific sum rules are based on Borel transforms, Hilbert transforms, positive moments or inverse moments. Sum rule calculations of the decay constants have been performed since the 1980's, with results in the range f B = 160 − 230 MeV, f Bs /f B = 1.1 − 1.4, and f Bc = 160 − 360 MeV [3]- [11]. A heavy quark effective theory calculation [12] gives f B = 206 ± 20 MeV, while historical LQCD determinations fall in the wide range f B = 161 − 218 MeV, and f Bs /f B = 1.11 − 1.16 [13]. Recent QCDSR calculations have narrowed the range of results to f B = 189 − 207 MeV, and f Bs = 216 − 242 MeV [14]- [15], while LQCD calculations, claiming high accuracy, give f B = 186 − 205 MeV, and f Bs = 224 − 244 MeV [16]. For B c a compilation of many LQCD determinations [17] gives values in the wide range f Bc = 380 − 680 MeV. The experimental situation remains unclear. For instance, using the recent Belle result for f B [18], and |V ub | = (4.01 ± 0.56) × 10 −3 [19], one obtains However, a recent value of the branching fraction of B + → τ + ν from BaBar [20] gives f B = 295 MeV (f B = 221 MeV), depending on their value of V ub from exclusive (inclusive) charmless semileptonic B-decays. In a previous calculation of f B and f Bs by some of the present authors [9], a method was used based on finite energy QCDSR (FESR) which equates positive moments of data with QCD theory. On the theoretical side, a large momentum expansion in powers of m 2 b /s was taken from reference [21], where m b is the mass of the bottom quark and s the square of the center-of-mass energy. The perturbative expansion was considered up to second order in the strong coupling constant, and up to seventh order in m 2 b /s. On the phenomenological side, a combination of positive moments involving Legendre polynomials was used to optimize the experimental data, enhancing the lowest lying B q meson. The contribution of the unknown continuum data was shown to be negligible, after a suitable choice of Legendre polynomials [9]. The present paper is devoted to determining f Bc , and updating results for f B and f Bs . We use more recent information on the QCD side of the sum rule, and employ a new criterion for optimizing the stability of the result. For the QCD correlator we use an O((m 2 b /s) 7 ) expansion of the pseudoscalar two-point function up to O(α 2 s ) [21]. For f Bc and f Bs we supplement this with O(m 4 q ) light quark mass corrections up to O(α s ) [8]. To account for non-perturbative corrections we include terms up to dimension six in the operator product expansion (OPE). Rather importantly, the correlator is expressed in terms of the running quark mass, rather than the pole mass [22], as it is well known that this improves the convergence of the perturbative series [8].

Preliminaries
To study the decay constants of B q mesons, with q ∈ {u, d, s, c}, we consider the pseudoscalar two-point correlator where j 5 (x) is the divergence of the axial-vector current and m b and m q are the masses of the bottom quark b and the lighter quark q, respectively. We then analytically continue Π(q 2 ) over the complex squared energy s-plane and invoke Cauchy's theorem valid for all holomorphic functions P (s), and all closed curves Γ which do not encircle a singularity of Π(s). The correlator has singularities only on the positive real axis for |s| > s thr , the physical threshold. Choosing Γ to correspond to a circle of radius |s| = s 0 , along both sides of the cut on the real axis, and using the Schwarz reflection principle, one finds This leads to a relation between QCD parameters and experimental observables, after invoking quark-hadron duality, i.e. assuming that Π(s) in the contour integral is given by QCD if s 0 is large enough, i.e.
The use of non-trivial integration kernels P (s) in FESR was pioneered in [23] in order to account for potential quark-hadron duality violations, e.g. in the Weinberg sum rules. Since then they have been frequently and successfully used in a variety of QCDSR applications. In particular, Legendre polynomial kernels subject to global constraints have been employed in extractions of the chiral condensates from τ -decay data [24], and on the chiral corrections to the Gell-Mann-Oakes-Renner relations [25]. The current most precise determinations of the charm-and bottom-quark masses [26]- [27] are also based on FESR involving these kernels. The general purpose of these kernels is to tune the FESR so as to emphasize or quench energy regions where the information is well or poorly known, respectively. In this way systematic uncertainties can be considerably reduced.

Phenomenological Contribution and P (s)
As usual we parametrize the phenomenological correlator with a single pole for the B q meson. For the unknown hadronic continuum we define the physical threshold where B * is a vector meson and P q is the lightest pseudoscalar meson with q = u, d, s, c quantum numbers, namely, P q = π, K, D, respectively. The spectral density can then be written as where ρ h (s) is given by the sum over all hadronic intermediate states with the quantum numbers of B q . To minimize the contribution of the unknown hadronic continuum we shall make a judicious choice of P (s). We choose a polynomial P n (s) = a 0 + a 1 s + a 2 s 2 + . . . + a n s n , and determine the coefficients, a i , subject to the global constraint To fix the last coefficient we use an arbitrary overall normalization condition. The functions P n (s) are then Legendre polynomials defined in the interval [s cont , s 0 ], where we fix s cont ∼ s phys by demanding maximum duality in the sense explained below. The introduction of this polynomial kernel in the sum rule minimizes the continuum contribution to the phenomenological side. In fact, to the extent that Im Π cont (s) can be approximated by an (n − 1)-th degree polynomial, these conditions lead to an exact cancellation of the continuum contribution on the left hand side of Eq. (9). At the same time, the role of the B q pole will be enhanced. Increasing n increases the OPE truncation error, but this is compensated by increasing the integration radius |s 0 |. It will turn out that the latter is of reasonable magnitude on the relevant scale, i.e., m 2 B . If, on the other hand, the radius of integration is chosen too high, the polynomial fit will deteriorate. It is hoped, and actually confirmed, that there is a wide intermediate region of stability. To find this region of stability we begin by choosing a value of n = 3, 4, 5, .. and then increasing s 0 until stability is reached in each case. In order to make the duality region more pronounced we allow the value of the threshold s cont to vary around the physical threshold s phys , as in [28], and in such a way that the dependence of f Bq on s 0 shows maximum stability. This means requiring that the first and second derivatives vanish at the same point. Reassuringly, we find that the value of s cont which maximizes this stability is very close to s phys , thus validating this procedure. The contribution of the threshold region is then estimated separately. We have reasons to believe that if the resulting f Bq is independent of n and s 0 in some finite region, then the unknown theoretical and hadronic contributions effectively cancel there. This follows from general properties of Legendre polynomials [29], the details of which will appear in a separate work. We are then left with the key relation

Perturbative Contribution
We now turn to the perturbative part of the QCD correlator, Π pQCD (s). To calculate f Bs and f Bc we need to take into account the effect of the lighter quark. We separate the correlator into its expression in the m q = 0 limit, and add to it the lighter quark mass corrections For Π pQCD mq=0 (s) we use the high-energy expansion given in [21]. This is an and is given to where L z = −(ln(−z))/2 and m q is the pole mass of the light quark (note that we keep the light quark mass in the pre-factor arising from the divergence of the axial-vector current). It is well known that the convergence of the QCD correlator is improved if it is expressed in terms of the running quark masses [8,9]. We thus write the pole masses in terms of the running masses, m b (µ) and m q (µ), expand the final expression as a series in α s (µ), and truncate at O(α 2 s ). As the perturbative expansion is not known to all orders in α s (µ) the final result will have some dependence on the renormalization scale, µ. To improve the convergence we take µ = m b , which re-sums the ln( m 2 b µ 2 ) terms appearing in the running quark mass [9]. We later investigate the effect of varying µ, retaining the ln( truncation has a very good convergence over the range of s used here, and it introduces a negligible error. To calculate at µ = m b , and to investigate the µ-dependence, we use the running coupling α s (µ) to four-loop order [30,31], and running quark masses to four-loops [32]. We performed the contour integral in two complementary ways, i.e. analytically and numerically. Analytically, we followed the method outlined in [9]. Numerically, we used standard computational methods to perform the contour integral. The final results agree to a few MeV, and the error estimates are entirely compatible.

Light Quark Mass Corrections
The light quark mass corrections were originally obtained to O(α s ) in [33,34]. These corrections to the full correlator are somewhat unwieldy, but the imaginary part to O(m 4 ) is presented in a compact form in [8]. For example, the O(α 0 s ) term is where x = m 2 b /s. The O(α s ) term can be found in the appendix of [8]. In the case of f Bc , where the lighter (charm) quark mass correction is almost 10 %, the O(m 4 ) terms contribute less than 1 %, so we are justified in making this approximation. Of the 10 % correction, around 4 % comes from the O(α s ) term, indicating that the convergence is not ideal. It is, however, good enough for the present level of precision, and in any case it is somewhat included in the determination of the error from varying µ in the range 3 -6 GeV. In the case of f Bs the light quark mass correction is around 5 %. For f B we take m u = m d = 0. To evaluate the integral of P n (s)Π pQCD mq (s) we again use the Cauchy integral theorem to rewrite it as the integral over the imaginary part along the real axis, now starting from s = (m b + m q ) 2 , the start of the cut, and use Eq. (7). This is purely a mathematical device, unrelated to the discussion leading to Eq. (9), so the result is exact and no error is introduced due to the poor convergence of the QCD expansion at low s.

Non-perturbative Contribution
The non-perturbative contributions to the correlator can be parametrized through Wilson's operator product expansion. In our calculation we include the effects of quark and gluon condensates up to dimension six [33]- [34] Note that we include the effect of the light quark mass only in the prefactor. The non-perturbative part contributes around 10 % to f B , and it is dominated by the lowest dimensional condensate qq . Its value is usually given at 2 GeV, and determined e.g. from the Gell-Mann-Oakes-Renner relation, with the result [25] qq = (−267 ± 5 MeV) 3 . To run this condensate to the renormalization scale µ we use the fact that m q (µ) qq is renormalization group invariant. For the strange quark condensate, and the mixed condensates we use where R sq = 0.6 ± 0.1 is from [36], and m 2 0 = 0.8 ± 0.2 GeV 2 from [37]. The gluon condensate has been determined from data on the hadronic decays of the τ -lepton [35], i.e. αsG 2 π = 0.07 ± 0.02 GeV 4 . The terms involving αsG 2 π and qσGq both contribute around 0.1 %, while the term involving qq 2 contributes around 0.001 %, so we are justified in ignoring higher dimensional operators. In f Bs the non-perturbative part plays a smaller role, while for f Bc it is negligible. The O(α s ) correction to the qq condensate term contributes around 1 % to f B , so this approximation is accurate enough. As with the perturbative part, we rewrite the pole mass in terms of the running mass [22], expand in powers of α s and truncate the series at O(α 2 s ). The resulting integral is easily evaluated using Cauchy's residue theorem.

Results
Given the complete QCD expansion, and the correct hadronic continuum, the final results would be independent of the two parameters n and s 0 . In practice, though, there is a finite region in this parameter space where results are stable. The wider this region, the more accurate the results. As explained earlier, we allow s cont to vary so that the first and second derivatives of f Bq with respect to s 0 vanish at some point. We then extract our prediction for each n from this inflection point. In To understand the effect of the variation of s cont we estimated the contribution of the continuum near the physical threshold in the determination of f B , where some phenomenological information is available. We considered the phase space contribution of the Bππ and B * π intermediate states near threshold, with form factors taken from chiral perturbation theory [38]- [39] and LQCD [40], respectively. In both cases the contributions are negligible compared with the dominant pole contribution of the B meson. We also investigated the contribution of a possible excited B state such as the B 1 (5721) predicted in [41] from heavy quark effective theories. The contribution of such a resonance will increase f B by 4 − 5 MeV, which must be considered as part of the systematic error. We expect a similar situation for f Bs and f Bc . The final result should be independent of the renormalization scale, µ. However, due to the various truncations of expansions in α s (µ) some µ dependence is introduced, which is roughly related to the convergence of the asymptotic expansion. To estimate the error we either remove or double the O(α 2 s ) correction. This is shown in the first line of table II.  [42], very close to the most recent and accurate value m b (m b ) = 4.171 ± 0.009 GeV [27], [43]. In table II we show the change in the decay constants when m b (M S) is changed by ± 1σ. Despite the relatively large uncertainty in the s-quark mass [39], [43], m s (2 GeV) = 94 ± 9 MeV, it only has a small effect on f Bs , reflecting the relative smallness of this correction. The uncertainty in the larger c-quark massm c (m c ) = 1.278 ± 0.009 GeV [43], [44] has also a small impact on f Bc .
Although the non-perturbative part contributes only around 10 % to f B , and less to f Bs , the errors of the condensates are not all negligible. The uncertainty in the quark condensate, qq (2 GeV), affects f B and f Bs by around 2 MeV and 1 MeV, respectively (the effect on f Bc is negligible). The relatively large uncertainty in α s G 2 /π has a small effect, while the errors in the higher dimensional condensates have a negligible effect, reflecting the fact that the quark condensate dominates the non-perturbative contribution. The uncertainty in m 2 0 has a non-negligible effect only for f B . There is some disagreement in the literature on the value of R sq . A recent LQCD result [45] claims R sq = 1.08 ± 0.16, which is in conflict with almost all determinations leading to R sq < 1 (for an exception see [46]). This value translates into f Bs = 230 MeV. The uncertainty in the strong coupling constant, α s (M Z ), has a small effect on the decay constants. The masses of the lowest lying pseudoscalar resonances are now very well known [42], and the only nonnegligible error arises from the uncertainty in M Bc . Adding these errors in quadrature we find f Bc = 528 ± 9 ± 17 MeV , where the first error comes from the inputs and the second is systematic, arising from the dependence on n, the O(α 2 s ) truncation, and the light quark mass expansion. The ratio of decay constants f Bs f B (which is unity in the chiral limit) can be determined quite accurately since many of the errors are correlated. We find f Bs f B = 1.19 ± 0.02 ± 0.04 , where again the first error comes from the inputs and the second is systematic. The effect of including the light quark mass terms in the perturbative expansion, Π (0)pert mq (s), is similar in size to the O(α 2 s ) correction for f Bs but is considerably larger for f Bc . The mass corrections increase f Bs by 11 MeV, whereas f Bc is increased by 45 MeV.

Conclusion
In this paper we discussed a new QCDSR determination of the heavy meson decay constant f Bc , and an update of f B and f Bs . We used Legendre polynomial integration kernels, and proposed a new technique allowing the continuum threshold, s cont , to vary. This variation leads to an effective suppression of the unknown hadronic continuum, and the unknown terms in the QCD correlator. For the latter we used an O(α 2 s ), and O(( m 2 b /s) 7 ) expansion of the pseudoscalar two-point function, supplemented by O(α s ), and O(m 4 q ) lighter quark mass corrections. Non-perturbative terms were incorporated through the OPE up to dimension six. For all three decay constants we found excellent stability over a wide range of s 0 and good convergence in n, the order of the Legendre polynomials.
The result for f Bc agrees with some LQCD, and a few determinations in other frameworks reviewed in [17], but it is around 20 % higher than that obtained in [47] from LQCD. However, it has been argued [48] that LQCD results for this constant may underestimate its value by some 23%. It must be pointed out that results for f Bq are very sensitive to the mass of the q-quark. For instance, in the case of f Bs the impact of a strange-quark mass m s 100 MeV results in a 20% increase in f Bs relative to f B . For f Bc , with m c 1.3 GeV, this rough argument would imply an increase over f B of some 200 %. New and more accurate determinations should resolve this issue. It should be stressed in closing that for f B and f Bs our results are perfectly consistent within errors with LQCD and recent experimental results.