Decay constants of charm and beauty pseudoscalar heavy-light mesons on fine lattices

We compute decay constants of heavy-light mesons in quenched lattice QCD with a lattice spacing of a ~ 0.04 fm using non-perturbatively O(a) improved Wilson fermions and O(a) improved currents. We obtain f_{D_s} = 220(6)(5)(11) MeV, f_D = 206(6)(3)(22) MeV, f_{B_s} = 205(7)(26)(17) MeV and f_B = 190(8)(23)(25) MeV, using the Sommer parameter r_0 = 0.5 fm to set the scale. The first error is statistical, the second systematic and the third from assuming a +-10% uncertainty in the experimental value of r_0. A detailed discussion is given in the text. We also present results for the meson decay constants f_K and f_\pi and the \rho meson mass.

1. Weak decays of heavy-light mesons with c and b quarks are interesting for studies of CP violation and determination of the CKM mixing angles. New experimental data on such decays are emerging (e.g. [1][2][3][4]) and their interpretation requires knowledge of hadronic matrix elements governed by the strong interaction. Lattice QCD allows one to calculate the strong matrix elements from first principles. However, if the heavy quark mass m Q is of the order of the inverse lattice spacing a, considerable discretization effects proportional to powers of am Q occur.
One possibility for coping with this problem is to use an effective theory such as Heavy Quark Effective Theory (HQET) [5] or Nonrelativistic QCD (NRQCD) [6]. These formalisms start from an infinitely heavy quark and consider corrections to this limit in the form of an expansion in the inverse of m Q . However, to study the charm quark in HQET or NRQCD requires a considerable number of correction terms, and one still has to worry about the uncertainty from the truncation of the 1/m Q expansion. A formulation for relativistic quarks where masses can be of O(1) in lattice units is the Fermilab approach [7] and modifications thereof, developed in [8,9] and [10]. One expects that the dominating discretization effects are then proportional to powers of momenta of O(Λ QCD ). While within HQET non-perturbative renormalization is possible [11], in many of the calculations using effective theories the renormalization constants are calculated only in perturbation theory (e.g. in Ref. [12]), leading to further uncertainties.
Another possibility is to simulate on very fine lattices, and this is the approach we have adopted in the present paper. We have performed a quenched lattice study of heavy mesons with a lattice spacing a of about 0.04 fm. On such a fine lattice a relativistic treatment of the charm quark should be justified and we expect that discretization errors are small compared to previous calculations on coarser lattices. We also make an attempt to study B mesons in our relativistic framework. Even on our fine lattice we cannot simulate B mesons directly, but the required extrapolation becomes relatively short-range. We expect that the resulting uncertainty is not much larger than the systematic error caused by the use of an effective theory. For example, recent unquenched calculations of f B and f Bs [12,13] employ NRQCD for the b quark and quote a ∼ 10% error based on perturbation theory and other systematic effects.
In this article we present results for the leptonic decay constants of the D s , B s , D and B mesons. We also evaluate light meson masses and decay constants to compare with previous quenched calculations of the light spectrum on coarser lattices and in order to be able to disentangle discretization and quenching effects.

2.
Our results are based on the analysis of 114 quenched Wilson gauge configurations simulated at the coupling parameter β = 6.6 with a mixed heatbath and microcanonical overrelaxation algorithm using the publicly available MILC code [14]. The lattice volume is 40 3 × 80, i.e. our lattice extends over 40 points (∼ 1.59 fm) in space and 80 points in time. The lattice spacing is determined using the Sommer parameter r 0 = 0.5 fm. This choice is motivated by a previous calculation [15] which used r 0 to determine the lattice spacings and found that the results for f Ds from a quenched lattice and a lattice with N f = 2 agreed (a ≈ 0.1 fm in these calculations). From the interpolating formula given in [16], one finds for our lattice a −1 = 4.97 GeV.
For the quarks we use the O(a) improved clover formulation [17], with the nonperturbative value of the clover coefficient c SW = 1.467 determined in Ref. [18]. We work with seven quark masses corresponding to three "light" hopping parameters κ = 0.13519, 0.13498, 0.13472 and four "heavy" hopping parameters, κ = 0.13000, 0.12900, 0.12100, 0.11500. Statistical errors are estimated by means of a bootstrap procedure using 500 bootstrap samples. For the central values we take the median. The error bars are calculated including 34% of the sample values below and above the median, respectively. Since the upper and lower error bars are found to be quite symmetric for most of our data, we just quote the larger of the two. The autocorrelation times for the pseudoscalar meson propagator appear to be small. In the worst case we studied, the autocorrelations decay after a distance of one configuration.
To extract the decay constants we follow the procedure described in Ref. [19]. For light and for heavy-light mesons we calculate the correlation functions where A 4 is the local axial vector current operator, P the pseudoscalar density which can be local (i = L) or Jacobi smeared (i = S), and V is the spatial lattice volume.  Masses and amplitudes are determined from fits of the correlation functions with where E is the ground state energy. In Table 1 we give the raw data for the light pseudoscalar meson masses, determined from C SL P P , and for light vector meson masses from smeared-local correlation functions of the spatial components of the vector currents.
To determine the bare quark masses, we calculate κ crit , the κ value corresponding to massless quarks, from a fit of the squared mass of a pseudoscalar meson ("pion") consisting of quarks with mass parameters κ 1 and κ 2 as a function of the averaged O(a) improved quark mass We use the non-perturbative value of −0.6636 for the improvement parameter b m using an interpolating formula from Ref. [20]. The fit includes all data with κ 1,2 ≥ 0.13472, where we find the improved quark masses to just lie on a straight line. We find κ crit = 0.135472 (11). The hopping parameter corresponding to the average u and d quark mass, κ ℓ , is determined by setting m P S on the left hand side of Eq. (4) equal to the physical pion mass, m P S = 138 MeV. We find κ ℓ = 0.135456 (10).
We parameterize the quark mass dependence of light meson decay matrix elements with hopping parameters κ 1 and κ 2 by fitting them to a function of the form The light quark mass dependence of masses and decay matrix elements of heavy-light mesons is parameterized using a linear fit as in Eq. (5), withm q being the light quark mass instead of the average quark mass. We also calculate the vector ("ρ meson") mass. The fit and the chiral extrapolation assuming a quark mass dependence as in Eq. (5) are shown in Fig. 1. At κ ℓ we find 846(37) MeV (the error is statistical), which is roughly a 10% (2σ) discrepancy with experiment. We compare our result to other recent quenched calculations in Table 2. Within errors our result agrees with Ref. [21], where a continuum extrapolation from coarser lattices with O(a) improved clover fermions is performed. We also list studies employing chiral lattice fermions where smaller quark masses can be reached while coarser lattices are used [22][23][24]. Ref. [23] quotes results from two lattice spacings. In Table 2 we present the results from their finer lattice. To determine the strange quark mass parameter κ s , we interpolate the vector meson mass to the physical φ meson mass, M φ = 1.01946 (19) GeV. We find κ s = 0.13502 (6). This is our "method I" for determining the κ value corresponding to the strange quark mass. Using Eq. (4) and setting m 2 K to the experimental value for (m 2 K + + m 2 K 0 )/2 gives a value in very close agreement: κ s = 0.134981 (9). The raw data for the heavy-light meson masses are given in Table 3. To find the physical values of the heavy-light meson masses, we extrapolate for each heavy quark mass linearly in m q , see Eq. (5). The fits are shown in Fig. 1. The quark mass dependence is linear to very good accuracy. This is in contrast to the findings of, e.g., Ref. [25]. In the final step, the calculation of the decay constants, the physical values of the c and b quark masses will be reached by interpolating or extrapolating the heavy-light meson mass to the D or B mass and the heavy-strange meson mass to the D s or B s mass.
In a quenched calculation, different methods to choose the input for determining physical parameters may give different answers. In order to investigate the influence of this arbitrariness we also use the heavy-light spectrum to determine κ s and call this procedure "method II". We consider the splitting between mesons with a heavy quark and a strange quark (generically denoted by M s ) and a meson with a heavy quark and a quark with the u, d quark mass (denoted by M ℓ ). In our data, as well as in experiment, the M s −M ℓ mass difference is fairly independent of the heavy quark mass. To fix κ s in method II, we choose a heavy quark close to the charm mass from our simulation points, namely κ = 0.129, and set the splitting between the M s and the M ℓ masses equal to the experimental value for the D meson, m Ds − m D = 98.85(30) MeV. The corresponding value for the strange hopping parameter is κ s = 0.134929 (15).
We calculate the pseudoscalar decay constants from the improved axial vector current A I where   [27]. Using a boosted coupling g 2 2143 which is close to the result one finds using the tadpole-improved scheme of [28]. A non-perturbative determination of b A on coarser lattices (β ≤ 6.4) [29] also gives values in agreement with boosted perturbation theory within errors.
The meson matrix elements of the currents are related to the amplitudes by where M denotes the meson mass. The convention for the factors of √ 2 corresponds to the normalization where f π ≃ 130 MeV.
3. The fit of the light meson decay constants according to Eq. (5) is shown in Fig. 2. Mesons with degenerate and nondegenerate quark masses fall on the same straight line. For f π , the value at the physical u, d quark mass, and f K , the value extrapolated to the averaged strange and u, d quark mass, we find This result for f π agrees well with the value of 137(2) MeV determined by [21] using an extrapolation to the continuum from coarser lattices. Both values are slightly larger than the experimental value of f π + = 130.7(4) MeV. Our value for f K is 6% or 2σ lower than the result from [21]     The SU(3) flavor breaking ratio of the light decay constants in our calculation turns out to be relatively small. We find f K /f π − 1 = 0.088 (12) .
Our number is substantially lower than the experimental value of 0.222, but is consistent with a recent quenched calculation using overlap fermions [30], which finds f K /f π − 1 = 0.09(4) using the same scale setting with r 0 = 0.5 fm. It is also consistent with other quenched determinations (see [31]).

4.
Next we consider the heavy-light decay constants.
To determine values at the physical quark masses, we extrapolate or interpolate the decay constants separately in the light and the to the c quark mass we use a formula motivated by HQET (see e.g. [32]). In the heavy quark limit, matching of the decay matrix element in the effective theory to the matrix element in full QCD introduces logarithmic corrections in the heavy quark mass which have to be resummed. In addition, power corrections in 1/m Q have to be added. Since the extrapolation to the b mass in our case is rather short, the precise form of the extrapolation formula is not important. We use only the lowest order running for α s , and take the heavy-light meson mass M ℓ (and not the quark mass m Q ) as an expansion (scale) parameter: Here γ 0 = −4 is the leading order anomalous dimension of the axial vector current, and b 0 = 11 is the leading coefficient of the QCD β function for zero dynamical flavors. The fits and the interpolated values are shown in Fig. 3. Decay constant ratios (18)(20) 1.080 (28)  The systematic errors are discussed in the text.
Our final results for the ratios of heavy-light decay constants are presented in Table 4, and the heavy-light decay constants are given along with a comparison in Table 5.
Estimation of systematic errors is notoriously difficult. One source of uncertainty concerns setting the quark masses to their physical values. For the strange quark this can be estimated by comparing the results from our methods I and II and suggests an error of 4 MeV for f Ds and f Bs . For the u, d quarks a chiral extrapolation is required. The corresponding error is difficult to estimate. Our data are consistent with the simplest linear chiral extrapolation. Quenched chiral perturbation theory provides a more sophisticated formula. However, it is not clear if it is applicable to our data. The uncertainty in fixing the heavy quark mass can be estimated by comparing the difference between the mass fixed from quarkonium and from the heavy-light meson system. Since the η c meson mass using the charm quark hopping parameter determined from the D s meson agrees with the physical value, we assume that this uncertainty is rather small in our calculation. In addition, for the B system there is an uncertainty from the extrapolation in the heavy quark mass. The difference between a quadratic fit to the matrix elements f √ M and a quadratic fit to Φ is very small and changes the values for the decay constants by less than 1 MeV. If only the three lighter heavy quark masses are included in the extrapolation to the b mass, f Bs (f B ) changes by −3 (+1) MeV.
Since we have results only from one lattice spacing, we cannot perform a continuum extrapolation from our data alone and have to estimate the discretization effects as a systematic error. Leading discretization effects are O(a 2 ). A rough estimate of them can be obtained by squaring the O(a) corrections appearing in the Symanzik improvement program. For the charm quark, the correction proportional to c A is small, around 2%, while the term proportional to the quark mass and b A is around 10% of the size of the matrix element itself. The square of the sum of these variations is around 1%, which we take as our estimate for the discretization error of f D and f Ds . A similar consideration for the B and B s systems results in an estimate of a discretization error of roughly 12%. For the error in the renormalization constants we use the estimate given for Z A in Ref. [26] of 1%. Since the heavy-light meson masses in lattice units in our simulation increase up to values of ∼ 0.8 one might be concerned about cutoff effects in the dispersion relation for the heavy-light meson. We therefore compare the kinetic mass [7] M kin calculated from where M is the rest energy of the meson. Results for the light quark mass close to the strange quark mass are shown in Fig. 4. We find that discretization errors in the dispersion relation are smaller than the statistical errors. The finite volume effects of the ratio f Bs /f B have been investigated in the framework of heavy meson chiral perturbation theory [33]. For quenched lattices of spatial extent 1.6 fm and pseudoscalar meson masses around 500 MeV (which corresponds to the smallest quark mass used in our simulation) they are quoted to be around 1% or smaller. It is plausible that the finite volume effects for D and D s mesons are of similar size.
We estimate the total systematic error due to discretization effects, errors in Z A , finite volume effects and ambiguities in fixing the physical quark masses by collecting all contributions and adding them in quadrature. It is given as the second error in Table 5.
The ratios are less sensitive to some of the systematic effects. The dominating ones for f Ds /f Bs and f D /f B are the discretization effects. For f Ds /f D and f Bs /f B we find a variation depending on how the strange quark mass is set, while the estimated discretization effects are smaller than the statistical errors. The uncertainties from fixing the physical quark masses and the discretization errors (added in quadrature) are given as the second error in Table 4.
We have used the value of the Sommer parameter r 0 = 0.5 fm to set the scale in physical units. This choice allows for a direct comparison with previous lattice determinations (see below) but is not universally accepted. With a different value of the Sommer parameter our results have to be modified accordingly. The variation of the decay constants if r 0 is changed by ±10% is given as third error bar for our results in Table 5 1 .

5.
Finally, we compare our results to other lattice calculations of decay constants. There exist recent quenched results for f Ds from nonperturbatively O(a) improved clover fermions [34][35][36] Ref.  [34] 0, clover, r 0 = 0.5 fm 243(2)( 03 24 ) 222(3)( 04 33 ) [35] 0, clover, r 0 = 0.5 fm 252(9) [36] 0, clover, r 0 = 0.5 fm 225(6) [9] 0, mod. Fermilab, r 0 = 0.5 fm 237(5) [25] 0, overlap, f π 266 (10) ) and for the B system (lower table). For the lattice calculations, the number of flavors in the simulation (N f ), the heavy quark (HQ) action, and the quantity used to set the scale are also indicated. The first error bar is the statistical, the second (where given) the systematic error except for the uncertainty in r 0 . For our work we quote a third error assuming a ±10% uncertainty in the physical value of r 0 . For the result from [36] we quote the value from the finest lattice instead of the continuum extrapolated result.
for a range of lattice spacings (0.03 ≤ a ≤ 0.1 fm) as well as for overlap quarks [25]. The comparison with the clover results is particularly interesting because it sheds some light on the discretization effects and might indicate the possibility of a joint continuum extrapolation. We plot the clover data in Fig. 5 as a function of the squared lattice spacing together with the overlap data. First we notice that on coarser lattices there is a discrepancy between the clover data of Refs. [34] and [35]. The discrepancy corresponds roughly to the difference one obtains when c A values from different nonperturbative calculations for a meson mass > 2.4 GeV are used, as discussed in [34]. Furthermore, the work [35] uses a nonperturbatively determined value for b A [37]. On the finer lattice of Ref. [34] (β = 6.2) the value used in Ref. [35] is about 6% larger than the perturbative number, which according to our estimates would affect the decay constants by at most 2%. At β = 6.0 the difference is even smaller. On the finest lattice used by [35] (β = 6.45) the difference between the perturbative and nonperturbative values of b A is ∼ 7%, which translates on a fine lattice into only a very small difference in the decay constants. In a more recent calculation [29] the nonperturbative value at that β value has come into agreement with perturbation theory, as mentioned in Section 2. Our data is in good agreement with the value obtained by Jüttner on his finest lattice [36]. The overlap value from [25] is on the other hand substantially larger. Being determined on a relatively coarse lattice it might be affected by discretization errors. It is important that all data shown in Fig. 5 come from lattices with similar spatial extent between 1.5 and 1.6 fm. So, finite size effects can be expected to be roughly the same in all calculations.  , diamonds), (ALPHA [35], squares), (Jüttner [36], circle), and overlap quarks (Ref. [25], triangle). The error bars show statistical and fitting uncertainties only. The scale is set using r 0 = 0.5 fm with the exception of [25] where f π is employed for the conversion of the decay constant to physical units. If r 0 = 0.5 fm is used instead, their lattice spacing decreases by 12%, which would increase their result for the decay constant even further.
Our result for f Bs is consistent with the quenched calculations of [38,39], but considerably lower than the nonrelativistic (but unquenched) calculation of [12]. The fit to the standard model gives a value with a relatively small error in between these two numbers.
For f D and f B , the values obtained from lattice calculations are consistent with the experimental results. Since the experimental errors are still large, this comparison is not conclusive, however.
6. Let us summarize our main findings. We have calculated decay constants of heavy-light pseudoscalar mesons on a very fine quenched lattice using clover fermions. Our extrapolations to the b quark mass appear reasonable. Nevertheless, from a comparison of the results at the charm mass to data obtained on coarser lattices we obtain the impression that discretization errors with the relativistic formalism adopted here are still significant for the b sector, unless the inverse lattice spacing becomes larger than ∼ 10 GeV.
Our results and those of Ref. [36] for f Ds are 10 − 15% smaller than the central values quoted for other recent lattice calculations, and roughly 20% smaller than recent experimental values. Eventually one would like to determine the decay constants to an accuracy of a few percent. Our work and the result of Ref. [36] indicate that discretization errors for the clover results on lattices with a −1 ≤ 2 − 3 GeV are too large to reach this precision, and that even a continuum extrapolation from a set of coarser lattices has a large uncertainty for heavy quarks. On the other hand, we do not find any source of large systematic errors, other than quenching, that could affect our calculation. It seems, therefore, that the new lattice results on fine lattices (this work and [36]) indicate a relatively small value for f Ds from lattice QCD. Quenching effects are notoriously difficult to estimate. However, since in previous calculations with a −1 ≈ 2 GeV [15] it was found that the quenching error is insignificant with our choice of lattice parameters, we expect that they will not be too large.
The systematic uncertainties on our results are larger for the B system than for the D system and more difficult to estimate reliably. Our results are in agreement with several other recent lattice calculations, but smaller than the values from recent unquenched calculations using nonrelativistic methods.
We find a rather small SU(3) symmetry breaking ratio of the heavy-light and light decay constants compared to experiment and also to several recent unquenched lattice calculations. The difference between our numbers and the unquenched results may be partially due to the use (see e.g. [42]) of a chiral extrapolation formula for the unquenched data which is inspired by chiral perturbation theory and predicts a particularly strong decrease of the decay constant at lighter quark mass values than are accessible in the simulation. This is in contrast to the use of a simple linear extrapolation in our calculation.