A proposal for B-physics on current lattices

A method to extract B-physics parameters (b-quark mass and fB, fBs decay constants) from currently available lattice data is presented and tested. The approach is based on the idea of constructing appropriate ratios of heavy-light meson masses and decay constants, respectively, possessing a precisely known static limit, and evaluating them at various pairs of heavy quark masses around the charm. Via a smooth interpolation in the heavy quark mass from the easily accessible charm region to the asymptotic point, B-physics parameters are computed with a few percent (statistical + systematic) error using recently produced Nf=2 maximally twisted Wilson fermions data.


Introduction
Heavy flavour physics is a corner of the Standard Model where chances are higher to uncover signals of new physics [1,2]. However, to extract from experiments useful phenomenological information, it is mandatory to have an accurate knowledge of the relevant hadronic matrix elements of the effective weak Hamiltonian. For low mass states (up to around the charm mass) lattice QCD (LQCD) represents the ideal framework where such calculations can be performed with well under control systematic errors [3].
Due to present day computer limitations, it is not possible, however, to work directly with the heaviest quarks (as the b-quark) propagating on the simulated lattice. Various strategies, more or less inspired to the heavy quark effective theory (HQET) [4,5], have been devised to circumvent this intrinsic difficulty, which go from non-perturbative matching of HQET onto QCD [6] to finite size scaling methods with relativistic heavy quark(s) [7]. Relativistic heavy-quark actions designed (highly tuned) to have reduced cutoff effects [8] have also been employed for this purpose. Encouraging results have been obtained by several groups [9] though different is the level at which the various relevant systematic effects are controlled.
In this work we wish to present a novel approach to B-physics in which the b-mass point is attained by interpolating from the charm region to the asymptotic infinite mass regime suitable ratios of heavy-light (hℓ) meson masses and decay constants, computed at a number of pairs of quark mass values lying slightly below and somewhat above the charm mass. The key feature of the approach is the use of ratios of physical quantities which by construction have a well defined and exactly known infinite h-quark mass limit. Injecting knowledge of meson masses and decay constants at the charm region, their h-quark mass evolution can be computed by a chain of successive steps up to values as large as about twice the charm mass. The b-physics region is finally reached through an interpolation from the simulated points to the exactly known infinite h-quark mass value.
A first test of the viability of the method is presented here. It has been carried out by exploiting the unquenched N f = 2 data recently produced by the ETM Collaboration [10,11] which makes use of maximally twisted Wilson fermions [12]. The results obtained in this feasibility study are very encouraging and compare nicely with the unquenched determinations today available in the literature [13] as well as with PDG numbers [1]. We get 2) f Bs = 235 (12) MeV , (1.3) where, as indicated explicitly in eq. (1.1), the b-mass has been run in a world with two (active) flavours. The results in eqs. (1.1) to (1.3) represent a "first principle" determinations of B-physics parameters with errors whose magnitude can be systematically reduced. The quoted uncertainty will be discussed in sects. 2 and 3 for eqs. (1.1) and (1.2)-(1.3), respectively. A few observations are in order here. First of all we would like to remark that the results above are extracted from unquenched LQCD data where u and d light fermions are dynamical, while heavier quarks are introduced only as valence quarks. This scheme is what goes under the name of "partially quenched" setting (see ref. [14] for a discussion within the twisted mass regularization of QCD). Systematic errors due to partial quenching are not included in the figures quoted in eqs. (1.1) to (1.3). The second observation is that no complicated renormalization steps are required for the method to work, because, as noted above, the necessary inputs are (ratios of) physical quantities (hℓ-pseudoscalar meson masses or decay constants) evaluated at h-quark masses around the charm region which are extracted from the existing (large volume) lattice configurations produced for the study of pion physics.
The central value ofμ ) in eq. (1.1) may look somewhat higher (though still compatible within statistical errors) than the available phenomenological estimates of the MS b-quark mass at its own scale, which lie in the range 4.2-4.3 GeV. However, it is not unlikely that, when the quenching of quarks heavier than u and d will be removed with the inclusion of dynamical s and, possibly, c quarks, m b will receive corrections which one can argue will tend to make the quantityμ ) somewhat smaller than the number given in eq. (1.1) 1 .
1 Indeed, if we evolve the intermediate resultμ The results (1.2) and (1.3), instead, are only affected at a level of less than 1% by our present uncertainties in the b-quark mass, because the hu/dand hs-meson decay constants happen to have a rather mild dependence on the h-quark mass.
The content of this paper is as follows. In sect. 2 we discuss the theoretical basis underlying the strategy that we propose to extract the value of the bquark mass from present day LQCD data and we provide a rather accurate determination of it with controlled errors. In sect. 3 we extend the method to the determination of the f B and f Bs decay constants. We conclude in sect. 4 with a few words on how to improve the quality of the numbers (1.1) to (1.3), and how to extend the present method to other hℓ-physics quantities the large h-quark mass behaviour of which is known. We defer to an Appendix some technical details concerning the way chiral and continuum extrapolations of hℓ meson masses and decay constants are performed.

b-quark mass
In this section we present a simple strategy aimed at determining the value of the b-quark mass through a smooth interpolation of suitable ratios of hℓ pseudoscalar lattice meson masses from the well accessible charm region to the asymptotic (infinite mass) point where these quantities have an exactly known value. Inspired by HQET results we consider the lattice ratios (a = lattice spacing) In eq. (2.1) and in the following by a "hat" we denote quark masses renormalized at 2 GeV in the MS scheme. Byμ ℓ we indicate the renormalized light quark mass, whileμ are pairs of (renormalized) "heavy" valence quark masses lying around (from below to somewhat above) the charm mass. The function ρ(logμ h ) is the factor that "transforms" the renormalized MS quark mass at 2 GeV scale into the so-called "quark pole mass". In formulae (2.2) In continuum perturbation theory (PT) ρ is known up to N 3 LL (i.e. up to next-to-next-to-next-leading-log) order terms included [15,16]. Finally N is 2 GeV to the b-quark mass scale by using anomalous dimension and β-function of the N f = 4 (rather than N f = 2) theory, the value of the b-quark mass gets lowered by about 3% compared to the value we give in (1.1).
the number of h-quark masses at which the values of the hℓ pseudoscalar lattice meson masses, M L hℓ , are supposed to have been measured. The choice of the form of eq. (2.1) is suggested by the HQET (continuum) asymptotic equation [4,5] lim Although the above constant is known to be 1, its value is not really needed here.
In order to simplify our subsequent analysis we keep fixed the ratio between two successive values of the heavy quark masses in eq. (2.1). Calling it λ > 1, we set the mass renormalization constant factor, Z −1 P , cancels out 2 .
Ratios of the kind defined in (2.1) are introduced with the idea that they might have a smoother chiral (µ ℓ → µ u/d , with µ u/d the light quark mass that yields the physical value of the pion mass) and continuum limit than each of the individual factors. Setting where we have introduced the (continuum limit) shorthand notation we observe that (for all λ > 1) eqs. (2.2) and (2.3) imply the following exact property lim x→0 y(x, λ;μ u/d ) = 1 . (2.6) From lattice data the function y(x, λ;μ u/d ) can be determined at certain discrete values of x (x (n) , n = 2, · · · , N). In order to extend our knowledge 2 We recall that in maximally twisted LQCD the twisted mass renormalizes according toμ = Z −1 P µ. If standard Wilson fermions were to be employed, the quantitym h = Z −1 S 0 (m 0h − m cr ), should be used in place ofμ h in eq. (2.1).
outside these particular points, while at the same time fully exploiting the strong constraining power provided by eq. (2.6), we imagine proceeding in the following way. Suppose the perturbative expansion of ρ has been computed and resummed up to N P LL order. Then we can define a tower of y-ratios, y| p , p = 0, 1, . . . , P + 1, such that provided ρ in eq. (2.1) is correspondingly taken at tree-level in the case of p = 0, or to N p−1 LL order for p > 0. Then for sufficiently small values of x we parameterize y| p in the form 8) where the coefficients η j (log x, λ;μ u/d ), j = 1, 2, are p-dependent (though to lighten the notation we do not display explicitly this dependence in the following), smooth functions of log x which tend to zero as λ → 1 and to some fixed constant as x → 0 [4,5]. With the ansatz (2.8) and at any order where PT results for ρ are available, it is not difficult to determine the η j coefficients from lattice data, assuming that their log x-dependence can be ignored in the range of masses where the above formulae are used.
It is important to remark that the ansatz (2.8) is based on the same kind of assumptions under which HQET is usually employed in the study of heavy quark physics. A posteriori, we check that the best fit values taken by the coefficient functions η j come out of a reasonable order of magnitude. Indeed we find that η 1 r 0 and η 2 r 2 0 are O(1) quantities.

Implementing the method
Let us start considering for concreteness the case where ρ is taken up to LL order and subsequently compare the results we get in this way with what one would obtain taking ρ at NLL-order or at tree-level (at tree-level ρ = 1). In order to determine the coefficient functions η j we proceed as follows. Let us make for the smallestμ h value the choiceμ (1) h = 1.230 GeV (we recall, we are referring to the MS scheme at the scale of 2 GeV). Fixing λ = 1.278 (see below) and, in this exploratory study, N = 4, we shall successively consider the h-quark masseŝ  (3), 5.22(2) and 6.61 (3) are employed at β = 3.8, β = 3.9 and β = 4.05, respectively [17]. The overall scale and the light quark mass are set by the experimental values of f π and m π via chiral fits of the pseudoscalar meson mass and decay constant data in the light quark sector [18,17]. Here we useμ u/d = 3.6(3) MeV and r 0 = 0.433 (14) fm.
Actually at each lattice spacing we will be dealing with the dimensionless quantitiesμ (j) h r 0 , j = 1, 2, 3, 4, of which the numbers quoted in eq. (2.9) represent the central values in physical units. Uncertainties on r 0 /a and Z −1 P (present at the level of about 3%) will be taken into account in the final error analysis.
From the set of the ETMC simulation data [10] with parameters detailed in Table 1, we extract the values of the hℓ pseudoscalar meson masses that correspond to theμ h values listed in (2.9). With these masses we construct the lattice ratios (2.1) on which a combined continuum and chiral fit is performed. As we hoped, ratios appear to have a mild dependence on the light quark massμ ℓ and small cutoff effects, as seen for instance in fig. 1. This makes our continuum and chiral fit straightforward and numerically robust.
The red squares in fig. 2 represent the numbers y h ) computed at LL order (i.e. p = 1 in eq. (2.7)). The best fit through the red squares and the point at x = 0 determines the values of the η j coefficients (j = 1, 2) and yields the middle (red) curve in the figure. We note that a second order polynomial in x is necessary to get a good fit to the data (a straight line forced to pass through the point y = 1 ) r 0 µ l < < β=4.05 β=3.90 β=3.80 Figure 1: . The black square with its error is the combined continuum and chirally (μ ℓ →μ u/d ) extrapolated value. Here and in all the following figures uncertainties possibly affecting the value of the variable in the horizontal axis are propagated to the quantity plotted on the vertical axis. at x = 0 would have a very large χ 2 ). The quadratic fit gives for the quantities r 0 η 1 and r 2 0 η 2 numbers of order unity, in agreement with the standard assumptions underlying HQET.
At this stage, having in our hands the quantities y (n) 1 = y(x (n) , 1.278;μ u/d )| 1 for any n (actually for any x), the iterative formula 10) should be looked at as a relation between the mass of the hu/d-meson, , and the corresponding heavy quark massμ , which is fully explicit if the initial, triggering value M hu/d (μ (1) h ) is assigned. The latter can be accurately measured, asμ (1) h lies in the well accessible charm quark mass region. We show in fig. 3 the quality of the continuum and chiral extrapolation of the triggering mass lattice data. Once this number is known, determination of the b-quark mass is tantamount to find the value of K at   (1) h .
( 2.11) A few related remarks are important here. 1) It is not really necessary to have the lattice hℓ pseudoscalar meson masses computed at values ofμ h matched exactly as indicated in eq. (2.9). A µ h interpolation between nearby M L hℓ masses can be carried out if necessary. This is what we have actually done in the numerical study we present in this paper. 2) It is not a priori guaranteed that eq. (2.10) can be solved for an integer value of the exponent K. This is not a problem, however, as one can always retune the parameter λ (and at the same time readjust the values in the sequence (2.9)), so as to end up with an integer for K b (this is the reason why the peculiar value λ = 1.278 was chosen). Alternatively one could adjust the starting value of the heavy quark mass or both. 3) A detailed discussion of the numerical analysis will be given in a forthcoming publication [19]. Here we only mention that a simple h ;μ ℓ , a) lattice data extracted from the simulations detailed in Table 1. The black square with its error is the continuum and chirally extrapolated value giving M hu/d (μ (1) h ) = 1.89(10) GeV.
SU(N f = 2) chiral perturbation theory NLO-formula was used to model thê µ ℓ dependence of the triggering hℓ meson mass and y-ratios, while O(a 2 ) effects have been parameterized (at eachμ h ) byμ ℓ independent terms. A few further details on this point are given in Appendix A.
Following the procedure outlined above, one finds the result given in eq. (1.1), or equivalently the renormalization group invariant (RGI) valuê 12) where in the running only two flavours are assumed to be active and the conventions of ref. [20] for the RGI quark mass have been used.

Discussion and error budget
It is important to check the degree of reliability of the key smoothness assumption we have been implicitly making on the function y(x) and test the sensitivity of the procedure and its result (2.11) to the order of PT at which the expansion of ρ is truncated. To this end we have repeated the entire analysis above using for ρ both a lower (tree-level) and a higher (NLL) order per- turbative approximation in place of the previously employed LL order truncation. We recall that in the large h-quark mass limit y| 0 = y(x, λ;μ u/d )| 0 and y| 2 = y(x, λ;μ u/d )| 2 approach 1 (see eq. (2.7)) with corrections O(1/ log x) and O(1/(log x) 3 ), respectively.
One finds a (very) little shift in the value of λ necessary for the solution of the corresponding iterative equation (2.10) to be an integer (from 1.278 to 1.273) if we go from y| 1 to y| 0 . The shift is instead totally negligible within our statistical accuracy if we move from y| 1 to y| 2 .
The fits to y| 0 and y| 2 data are shown in fig. 2 (lower blue and upper green curve) together with the fit to y| 1 data (middle red curve). One clearly sees that the y| p curves tend to become flatter and flatter as we move from p = 0 to p = 2. As for the values of the b-quark mass, instead, the number extracted from y| 0 -ratios is only about 2% smaller than the one obtained using the y| 1 -ratios. The difference between the latter and the one extracted from y| 2 -ratios is smaller than 1%.
The stability of the value of the b-quark mass with varying p should not come as a surprise. It is enough to notice that, if the y| p -ratios were exactly known, for a generic value of p one would get (recall eq. (2.10)) 13) as all the intermediate ρ factors (except the first and the last) cancel out in the l.h.s. leaving behind simply the product of the y| 0 -ratios. The small p dependence we have found in the value of λ (hence inμ b ) is due to the slightly different level of accuracy by which the y| p -ratios (which instead significantly depend on p) can be described by a polynomial in x with the lowest order coefficient set to unity. In this respect, increasing p is expected to improve the quality of the ansatz (2.8) and reduce the systematic error associated to it. This is so until p becomes so "large" that the accuracy of the ρ estimate gets spoiled by the renormalon ambiguity in its perturbative expansion [21].
To account for the truncation to LL of the perturbative expansion for ρ we have conservatively decided to attribute to the b-quark mass value a systematic error of 1%, which is added in quadrature to the other errors discussed below, leading to the total error quoted in eq. (1.1) (and (2.12)).

Error budget
The total error we attribute to the b-quark mass results (1.1) and (2.12) takes into account a number of statistical and systematic effects which we now briefly illustrate. The relative error on the product of the continuum y-ratios in the l.h.s. of eq. (2.10) is only about 1%, whereas the pseudoscalar meson mass in the charm region (M hu/d (μ (1) h ) in the r.h.s. of eq. (2.10)) contributes a relative error of about 5%. These errors are the result of our statistically limited knowledge of hℓ-meson correlators, r 0 /a and Z P 3 as well as of a number of further systematic errors. Among the latter we mention those coming from the fit ansatz underlying the combined continuum and chiral (μ ℓ →μ u/d ) extrapolation, the error due to the x-interpolation to the bmass point, as well as the (tiny) error inherent the numerical solution of eq. (2.10) (giving K b = 6 and λ = 1.278). As we discussed above, the effect onμ b due to the truncation of the ρ perturbative series to order p is very small, not larger than 1%. Another .5-1% systematic error comes from the possible (neglected) logarithmic dependence of the η j , j = 1, 2 coefficients. The relative uncertainty on η j associated with these effects can be estimated to be O(α s (1/x)) ∼ 10 − 15%, a number which is never larger than the statistical errors on their best fit values. Finally cross-correlations between the different quantities (stemming from common ensembles of gauge configurations) are as usual taken into account by a bootstrap error analysis. Further technical aspects of the error analysis are deferred to ref. [19].
The information provided in figs. 1, 2 and 3 about the a 2 ,μ ℓ and x dependence of the intermediate quantities entering our analysis as well as about the precision in solving eq. (2.10) (see fig. 4) shows that the global systematic uncertainty is well within (or below) our present statistical errors.
We conclude by observing that, as expected, our results for the b-quark mass (and f B or f Bs discussed in the next section) do not significantly depend on the value of the intermediate quantity r 0 which is only employed to ease continuum extrapolations, while the physical scale is ultimately set by f π .

b and c quark masses
Although not necessary, the phenomenological value of the D-meson mass could have been used as a triggering mass. In this case µ (1) h would have to be identified with µ c .
We note in this context that, since we get (see the black square in fig in very good agreement with other estimates 4 . In closing this section we note that an independent determination of the b-quark mass can be obtained repeating the same analysis as before but using M Bs , instead of M B , and replacingμ u/d withμ s . By doing that we find a result which is fully consistent with the one in eq. (2.12). Alternatively, and perhaps more interestingly, one could use K b as determined from M B to predict M Bs , or better the ratio M Bs /M B , by the method we are proposing in this paper. Such an analysis is in progress and will be presented elsewhere [19]. 4 It is interesting to compare our unquenched result for the ratio of RGI masses m (13) (7) determined using the quenched data of ref. [22] (for the c-mass) and [23] (for the b-mass).
More recently the work of ref. [24] has appeared where the number

f B and f B s decay constants
A strategy very similar to the one outlined in sect. 2 can be employed to extract accurate values of the f B and f Bs decay constants from available lattice data. In analogy with what we have done before, one should now take with the (continuum limit) shorthand notation where the quark massμ ℓ must be extrapolated to eitherμ u/d or to the appropriate strange quark mass value,μ s , depending on whether one wants to compute f B or f Bs . ) r 0 µ l < < β=4.05 β=3.90 β=3.80 Figure 5: The black square with its error is the combined continuum and chirally (μ ℓ →μ u/d ) extrapolated value.
The form of the function z(x, λ;μ ℓ ) is dictated by the continuum asymptotic formula

2)
which follows by matching HQET to QCD [4,5]. The presence of the factor C stat A comes from the fact that in HQET the axial (and vector) current needs to be renormalized. The renormalization constant C stat A is known in PT up to three loops [25]. The ratio of ρ factors (raised to the appropriate power) is there to convert MS heavy quark masses to pole masses (see eq. (2.2)).
In analogy with what we did in sect. 2 in determining the b-quark mass, with the purpose of checking the robustness of the procedure, we shall take C stat A and ρ at increasing orders in PT, from tree-level up to NLL order, and construct z| p -ratios endowed with the asymptotic behaviour

5)
Just like in the case of the determination of the b-quark mass, the values of f B or f Bs that we shall extract will be almost independent of the PT truncation order.

Implementing the method. The case of f B
The z-ratios (3.1) have been evaluated at the reference h-quark masses of the list (2.9) for each of the lattice spacings and light quark mass values given in Table ( 1). When we perform the continuum and chiral extrapolation of the ETMC lattice data for the ratios (3.1) of hℓ pseudoscalar meson decay constants (again based on simple chiral NLO-formulae supplemented withμ ℓindependent O(a 2 ) corrections -see Appendix A), as hoped, a rather smooth behaviour is found since most of the a 2 andμ ℓ dependence gets canceled in taking the ratio. The observable dependence onμ ℓ and a 2 is mild and/or hardly significant within our present statistical errors (see e.g. fig. 5). From the structure of eq. (3.1) one derives the iterative formulae (analogous to eq. (2.10) with z (n)  Similarly to what we did in fig. 2, we collect in fig. 6 continuum and chirally extrapolated data for z| p , p = 0, 1, 2 and best fit curves through these data and the value at x = 0. Thus, for instance, the middle (red) curve is the parabola (eq. (3.4)) which best fits the values of z (n) 1 = z(x (n) , 1.278;μ u/d )| 1 , n = 2, 3, 4, at the heavy quark masses (2.9). The red vertical line marks the position x b which corresponds to the previously determined value ofμ b (eq. (2.12)) and crosses the curve at the point z With the help of this number and the values of z (j) 1 for 4 < j ≤ K b + 1, eq. (3.6) provides a determination of f hu/d (μ b ) in terms of f hu/d (μ (1) h ) (with LL-accurate fit for the z-ratios). As observed before, the latter does not necessarily has to be identified with the phenomenological value of f D for the method to work, as what we actually need to know is the dependence of f hu/d (μ h ) onμ h at around the charm mass. Nevertheless, since, as remarked in sect. 2.2, M hu/d (μ  (9) MeV given in ref. [11]. The latter was obtained in the standard way (see ref. [10]) from the same ETMC gauge configuration ensembles, but with a rather different analysis method where the meson masses rather than the renormalized quark masses were kept fixed as a → 0, resulting in somewhat different statistical (no use of Z P ) and systematic errors as compared to the present study. For the present computation of f hu/d (μ (1) h ) = f D the quality of the continuum and chiral extrapolation of our lattice data is shown in fig. 7. Taking as triggering value the continuum and chirally extrapolated value of the pseudoscalar decay constant computed at µ (1) h , we get which is precisely the result (1.2).
h ;μ ℓ , a) lattice data extracted from the simulations detailed in Table 1. The black square with its error is the continuum and chirally extrapolated value giving f hu/d (μ (1) h ,μ u/d ) = 211(9) MeV.

Discussion and f B error budget
To test the reliability of the interpolation of our trial functions, z(x), to b-mass point, we have explicitly checked the stability of the result (3.7) to the order of PT at which the expansion of ρ 1/2 /C stat A is truncated. For this purpose we have repeated the whole previous analysis employing values of ρ 1/2 /C stat A computed at tree-level (p = 0) and NLL (p = 2) order. Upon comparing with the decay constant values obtained in these other ways we see that numbers obtained using the z| 0 -ratios (upper blue curve in fig. 6) differ by less than 1% from the value one gets from the red data (LL z| 1 -ratios). If one employs the lower green data (coming from the NLL z| 2 -ratios) the difference with the previous determination is totally negligible (about 0.1%). This specific systematic effect on f B was hence conservatively estimated to be ∼ 0.5% and added quadratically to the full error.
As in the case of the determination of the b-quark mass, the remarkable numerical stability of f B with varying p can be traced back to the good quality of the interpolation ansatz (3.4) and the relation (again valid for exactly known z| p -ratios) . (3.8)

The f B error budget
The total error we attribute to f B in eqs (1.2) and (3.7) comes in almost equal parts from the product of z-ratios in the l.h.s. of eq. (3.6) and the value of f hu/d (μ (1) h ) and is a combination of statistical and systematic errors stemming from the same sources already illustrated in the case of the b-quark mass in sect. (2.2). As we saw above, the systematic error stemming from the truncation of the PT expansion of ρ 1/2 /C stat A has a negligible impact on f B . Another 1-2% systematic uncertainty comes from the possible (neglected) logarithmic dependence of the ζ j , j = 1, 2 coefficients. The relative uncertainty on ζ j associated with these effects can be estimated to be O(α S (1/x)) ∼ 10 − 15%, a number which, as in the case of the η j 's entering our analysis for µ b , is never larger than the statistical errors on their best fit values. In any case inspection of figs. 5 to 7 shows that all systematic uncertainties are smaller than our current statistical errors.

The case of f Bs
In order to come up with a determination of f Bs one has simply to repeat the whole procedure setting µ ℓ → µ s . With reference to the value ofμ M S,N f =2 s (2 GeV) = 99 (7) MeV given in [26], one finds (thanks to the equality ofμ (1) h with the charm quark mass) f Ds = f hs (μ (1) h ) = 252 (7) MeV (3.9) and the best fit z-ratio curves shown in fig. 8. We may quote as our final result (see also eq. (1. 2)) f Bs = f hs (μ b ) = 235 (12) MeV , (3.10) where errors are estimated as in the case of f B . We also note that the result (3.9) for the decay constant of the D s meson is in agreement with that of ref. [11] and contributes to further reduce the possible "tension" with the recent Cleo data reanalysis [27].
A more complete analysis of the many possibilities and refinements one can envisage will be presented in [19]. As is clear by comparing the results for f B and f Bs (eqs. (3.7) and (3.10), respectively), our method yields a significantly smaller error for the decay constant of the strange B-meson. The reason is that no large statistical fluctuations from the light (u/d) quark propagators nor (valence) chiral extrapolation uncertainties enter the computation of f Bs . In view of this observation we remark that, if one would know with high accuracy the ratio f B /f Bs , a more precise determination of f B could be obtained by simply multiplying this number by f Bs 5 . Actually the quantity f B /f Bs can be accurately computed by a simple generalization of the method discussed in this paper. It is indeed sufficient to consider the double ratio (3.11) and follow the procedure we described before starting from the triggering quantity [f hu/d /f hs ](μ (1) h ). This kind of analysis is under way and will be discussed elsewhere [19].

Conclusions and outlook
In this paper we have proposed a novel strategy to determine B-physics parameters from currently available Wilson fermion simulation data. As a first test of the method, we have computed in (the continuum limit of) QCD with N f = 2 light dynamical quarks the (renormalized) b-quark mass as well as the B-meson decay constants, f B and f Bs , employing the gauge configurations recently produced by the ETM Collaboration with maximally twisted Wilson fermion action [10].
The method provides rather accurate numbers with errors that are dominated by the uncertainties related to the limited statistical accuracy by which the (two-point) hℓ pseudoscalar meson correlators and the quark mass renormalization constant, Z −1 P , are evaluated. A better assessment of the systematic errors due to the limited knowledge of logarithmic corrections can only come from data taken at quark masses larger than the ones displayed in eq. (2.9).
In several respects the present feasibility study could benefit from the nice properties (particularly O(a) improvement [29,14]) of maximally twisted Wilson fermions. Indeed, an important feature of the present computation is the pretty good control we have of cutoff effects, which (judging from the spread between values at the coarsest lattice spacing and those at the continuum limit) are always smaller than 10%. This is so both for the triggering quantities at the charm mass scale and for the y-and z-ratios, which involve higher quark masses (up to twice the charm mass).
It is also interesting to note that the whole procedure only relies on the use of physical quantities that can be easily determined from lattice simulations, while the need for a renormalization step is limited to establishing the relation between the renormalized charm-like mass and the values of the triggering pseudoscalar meson mass. Fixing this relation requires the knowledge of Z P . No extra renormalization factor is needed for the calculation of the decay constants of interest if maximally twisted fermions are used as the charged axial currents are exactly conserved at finite lattice spacing.
There is a lot of room for improvement in the application of the method, like reducing the statistical error of the correlation functions, using several, suitably smeared meson sources, increasing the accuracy by which Z P is known and incorporating in the analysis the new ETMC set of data that are coming out at a finer lattice spacing (β = 4.2).
Needless to say, the method can be straightforwardly extended to LQCD computations with u, d, s and possibly c dynamical quarks where quenching uncertainties are virtually absent. In this respect we wish to note that in simulations with N f = 3 dynamical quarks, although low energy hadronic effects in the B-meson wave function are correctly treated, a conflict remains between the number (N f = 3) of dynamical quarks running in the loops (and thus relevant for the subtraction of UV divergencies) and the number (N f = 4) that instead should be used for continuum RG-evolution at scales above, say, 1.5 GeV. This problem and the related RG-uncertainties are completely removed if also the c quark is made dynamical.
Finally, we remark that the strategy we have outlined can be applied to any other hℓ physical quantity the large µ h behaviour of which is known (typically from large quark mass arguments).

Fit ansatz for hℓ meson masses and decay constants
We have modeled theμ ℓ -dependence of M hℓ and f hℓ in a form that is consistent with the known results of (NLO) SU(2) chiral effective theories for pseudoscalar mesons made up by a light plus a heavy quark. The form of the fit ansatz was chosen generic enough so as to encompass (see discussion below) the expectedμ ℓ -dependence both in the case when the heavy quark is treated as static [30,31,32] and when the latter is considered non-light but still relativistic [33,34].
In r 0 -units we write for masses and decay constants where the last terms, with parameters C L and D L , have been included to cope with the expected O(a 2 ) discretization effects. The fit parameters C ... 's and D ... 's in general depend onμ h , though the form of the fit ansatz (A.1) and (A.2) has been actually employed only for a fixed value of the heavy quark mass, namely forμ h =μ (1) h , when we evaluate the so-called triggering meson mass and decay constant, respectively.
We checked that, within the statistical accuracy of our data, noμ ℓ dependence is visible in the O(a 2 ) terms in eqs. (A.1) and (A.2). The fit ansatz for M hℓ r 0 does not include logarithmic terms. This is consistent with the results of the chiral effective theory for hℓ pseudoscalar mesons with a light plus a non-light and relativistic quark [33,34], but it can be equally well regarded as a simple Taylor expansion leading to a polynomial interpolation of data points with a very smooth dependence onμ ℓ (see fig. 3). In this sense the fit ansatz for M hℓ r 0 is also consistent with the spirit of the effective theory for static-light mesons [30,31,32], where theμ ℓ dependence of M hℓ is expected to be a tiny effect (as we indeed find).
The coefficient d 1 , multiplying the term ∼μ ℓ log(μ ℓ ) in eq. (A.2), was taken as a free fit parameter. Numerically we get for d 1 /D 0 atμ h =μ (1) h a value (−1.0 ± 0.4) which falls in between (and agrees within statistical errors with) what is expected from the arguments of refs. [30,31,32] (where the heavy quark is treated as a static source) and those of ref. [33] (where it is taken as a relativistic particle). The result we find for d 1 is not surprising as our heavy quark mass lies in the charm region.
The low energy constant f 0 and B 0 have been taken from recent ETMC analyses of light meson quantities [17]. which enter the ratios y and z at the variousμ h -values, we employed fit ansatz analogous to (and derived from) eqs. (A.1) and (A.2) above. They all are of the following form: a leading term plus a term linear inμ ℓ and another one proportional to a 2 (andμ ℓ -independent). Note that, if we assume that d 1 does not appreciably vary asμ h changes by the factor λ ∼ 1.27-1.28 (which is a natural expectation in any effective theory for hℓ mesons), one finds that in the ratios (A.3) the possibleμ ℓ log(μ ℓ ) dependence cancels at NLO, and is pushed to NNLO.