{'fraction_non_alphanumeric': 0.03595360229215768, 'fraction_numerical': 0.017453055440023416, 'mean_word_length': 4.3085289066971955, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Properties of B-Mesons in Lattice QCD

The results of an extensive study of B-meson properties in quenched lattice QCD are presented. The studies are carried out in the static quark limit where the b-quark is taken to be infinitely massive. Our computations rely on a multistate smearing method introduced previously, with smearing functions generated from a relativistic lattice quark model. Systematic errors arising from excited state contamination, finite volume effects, and the chiral extrapolation for the light quarks are estimated. We obtain continuum results for the mass splitting M_{B_s}- M_{B_u} = 86 (+/-)12(stat) {+7/-9}(syst) MeV, the ratio of decay constants f_{B_s}/f_{B_u} = 1.22 (+/-)0.04(stat) (+/-)0.02 (syst). For the B-meson decay constant we separately exhibit the sizable uncertainties in the extrapolation to the continuum limit a ->0 and higher order perturbative matching. We obtain f_{B} = 188 (+/-)23(stat) (+/-)15(syst) {+26/-0}(extrap) (+/-)14 (pert) MeV. ----- [Postscript version of paper available by anonymous ftp at fncrd6.fnal.gov. The file is fb.ps in subdirectory theory.]


INTRODUCTION
Heavy-light mesons are ideal systems for lattice QCD studies. For hadrons that contain a single heavy quark, the dynamics of QCD simplifies as m Q → ∞. In this limit, the heavy quark is on shell and at rest relative to the hadron [1] and the QCD dynamics becomes independent of m Q [2]. Heavy quark mass dependence can be extracted analytically to produce an effective action for a static quark and the remaining light degrees of freedom [3,4].
In the continuum, the resulting heavy quark effective theory (HQET) makes transparent the symmetry [5] and scaling relations between systems which differ by heavy quark spin or flavor [5,6,7]. Furthermore for finite heavy quark masses, the effective action can be improved order by order in 1/m Q . The m Q → ∞ limit may also provide other insights into QCD dynamics. In particular, heavy-light mesons are physical systems with a single light valence quark and fully relativistic QCD dynamics. It is likely that much can be learned about constituent quark ideas in this simple setting [8].
These theoretical developments have immediate physical applications for B physics.
Since the b quark mass is significantly heavier than the other mass scales (the QCD scale and the light quarks masses) which enter into the dynamics of B hadrons, it is likely to be a good approximation to treat the b quark in the m Q → ∞ limit within B hadrons.
Recent developments in lattice gauge theory have led to the possibility of calculating the masses and decay constants of B q (q=u, d, s) mesons from first principles (QCD) with enough accuracy to be of both phenomenological and theoretical interest. In particular, we consider the mass difference M Bs − M Bu , the decay constant f B u,d and the ratio f Bs /f B d .
We will only consider the lattice action appropriate to the static limit (m b → ∞) and hence all the results reported have corrections of order Λ QCD /m b . A variety of other methods have been developed to study B mesons on the lattice. Methods for treating the b quark using nonrelativistic actions [9], a non-zero velocity formulation [10], and a generalization of the usual Wilson action which is not constrained to quark masses less than the inverse lattice spacing [11] are being actively pursued.
The present study encompasses a systematic analysis of data for M B and f B at four different lattice spacings a (with associated beta values β = 5.7, 5.9, 6.1, and 6.3) and a variety of physical volumes (in lattice units 12 3 , 16 3 , and 20 3 ) for one fixed spacing (β = 5.9).
For each case above at least four light quark mass values (kappa values corresponding to pion masses in the range 300 -800 MeV) were studied. A uniform and consistent fitting scheme was used in all cases. This is particularly important for an accurate extrapolation to results at physical light quark masses (κ = κ u,d and κ = κ s ). and for assessing the a-dependence of the results.
In the static approximation, the heavy quark propagator is reduced to a straight timelike Wilson line, making it possible to calculate correlators of spatially smearedQq operators without having to compute light quark propagators from smeared sources. Because of this simplifying feature, the heavy-light meson system is an ideal place to develop sophisticated operator smearing techniques. Such techniques are indispensable for the accurate calculation of f B and other B-meson properties. Until recently, most such calculations have relied on a more-or-less ad hoc choice of smearing functions (e.g. walls [12], cubes [13], or exponentials [14], or Gaussians [14]).
In the present study, we have applied a multistate smearing method introduced previously in Ref. [16]. This analysis provides some significant improvements over previous investigations. First, we have made a serious effort to construct smearing functions which closely resemble the actual Coulomb gauge wave functions of the valenceQq system as measured on the lattice. As reported in a previous paper [17], the heavy-light wavefunctions from lattice QCD are reproduced with remarkable accuracy by a simple relativistic quark model (RQM) Hamiltonian which contains the static QCD potential extracted from Wilson lines in Coulomb gauge. In addition to being an interesting statement about QCD dynamics, the success of the RQM has a practical consequence which we will exploit here. The RQM Hamiltonian provides a simple and precise way of constructing orthonormal sets of realistic smearing functions for lattice heavy-light calculations. For each value of β and lattice size, the static Wilson potential is calculated from the gauge configurations and used in the RQM to generate heavy-light smearing functions. The only tuneable parameter in this procedure is the light-quark constituent mass parameter µ in the kinetic term of the RQM Hamiltonian. In practice, this parameter was initially selected by measuring the lattice QCD ground state wavefunction and adjusting µ to give the best fit for the RQM ground state. In some cases, after an initial multistate fit to the heavy-light propagators, it was found useful to iterate the procedure with a more finely tuned value of the constituent quark mass µ, using the more accurate wave functions obtained from the multistate fit.
In addition to this method for constructing smearing functions, another important inno-vation introduced in the present study is the fitting procedure used to extract information from the heavy-light correlators. Starting with the wavefunctions from the first M S-wave states of the RQM, it is relatively easy to construct the entire M × M matrix of correlators among the corresponding smearedQq operators, as well as the "smeared-local" correlators between each of these operators and the localQq source. Along with the local-local correlator, these form an (M + 1) × (M + 1) matrix. This matrix contains far more information than just the smeared-smeared and smeared-local correlators of any single smearing function. In particular, the matrix contains information about excited states, which, when properly exploited, allows an accurate extraction of ground state properties even at very short time separations, where excited state contributions are still large. The method we introduce to accomplish this employs a χ 2 minimization procedure to simultaneously fit the (M +1)×(M +1) matrix of correlators to a sum of M exponential (pole) terms, representing the contribution of the M lowest lying heavy-light eigenstates. (In all fits, we exclude the local-local correlator, which, at short time separations, is not well fit by a few low lying states.) The matrix coefficient (residue) of each pole term factorizes and can be written in terms of an M + 1 component vector whose entries represent the vacuum-to-eigenstate matrix element of each smeared operator. In practice we have found an M = 2 fit to yield fairly accurate results for ground state properties. For all of the fits used, the χ 2 per degree of freedom was less than 1.3. The multiparameter fits were carried out using the CERNLIB minimization routine MINUIT.
The improved control over systematic errors gained from the multistate fitting method allows us to better address a number of issues. In particular the dependence of heavy-light meson parameters on both the light quark mass and the lattice spacing are examined in detail. One of the difficulties with previous analyses which prevented accurate chiral and a → 0 extrapolations was in the arbitrariness of the smearing procedure. It is clear that any ad hoc smearing function will have a substantial overlap with excited states. Typically one tries to optimize the smearing function (e.g. by adjusting the size of the cube) and to go far enough out in time that excited states have died away. The approximate equality of smearedsmeared and smeared-local effective masses, combined with some indication of an effective mass plateau, are the main criteria of success in this procedure. Unfortunately, the procedure is somewhat subjective and it is difficult to rule out large systematic errors due to excited state contamination. An attempt to reduce these errors by extracting results from larger time separations leads to a rapid deterioration of statistics. Moreover, because of the wellknown signal-to-noise difficulty for the heavy-light propagator in the static approximation, [15,16], the problem of isolating the ground state becomes more difficult at smaller lattice spacing. As a result, extrapolation to a = 0 is particularly problematic. Furthermore, use of any fixed smearing function at different light quark masses introduces a significant systematic error in the extracted kappa dependence. This is important in the determination of results for f Bs /f Bu and M Bs − M Bu . The multistate fitting procedure effectively deals with these difficulties, greatly reducing our errors.
Extrapolating to the continuum limit (a = 0) we obtain the ratio f Bs /f Bu = 1.22 ± .04(stat) ± .02(syst) and M Bs − M Bu = 86 ± 12(stat) +7 −9 (syst) MeV for the mass difference. For these quantities, only a slight dependence on the lattice spacing is observed, and the systematic errors associated with the a → 0 extrapolation (included in the above) are small.
The situation for the decay constant is more complicated. We find a significant lattice spacing dependence for the ground-state pseudoscalar decay constant f B . The results for the four β values are consistent with either a linear or quadratic dependence on the lattice spacing a. The linearly extrapolated result in the a → 0 limit is f B = 188 ± 23 ± 15 MeV This result is notably smaller than previous estimates of f B in the static approximation.
The primary reason for this is the a → 0 extrapolation. The quadratic extrapolated result in the a → 0 limit is f B = 214 ± 13 ± 17 MeV. This fit reflects the fact that our results at β = 6.3 are not inconsistent with those of Ref. [12], [18], [19] (the latter two results using a clover action for the light quarks). In both cases, the first error is statistical while the second includes systematic errors, which we discuss in Section 5. Our final result for the B u,d meson decay constant (in the heavy quark limit) f B = 188±23(stat)±15(syst) +26 −0 (extrap)± 14(pert) MeV explicitly separates out this theoretical uncertainity associated with the a → 0 extrapolation as well as our estimate of the uncertainity associated with higher order perturbative matching corrections.
We will give a complete discussion of the perturbative matching of lattice to continuum results in Section 2. We include a discussion of the heavy quark mass renormalization in relation to the residual mass parameter (Λ = M B − m b (quark)) of HQET. In Section 3 we discuss our analysis procedure. The details of the multistate smearing technique and the construction of the smearing functions from a relativistic quark model (RQM) are presented.
In Section 4 the numerical lattice results at each β are presented. The statistical and fitting errors associated with our final physical results are determined. The light meson results used to set light quark masses are contained in Appendix A. The discussion of systematic errors associated with excited state contamination, finite volume, nonzero light quark masses, scale uncertainities, and the extrapolation to zero lattice spacing are all discussed in Section 5. A study of the time evolution of the wavefunctions for heavy-light states is presented in Appendix B. These results provide an independent check that our multistate smearing analysis has removed excited state contamination. In the Section 6 we present our final results, compare them with other recent calculations, and discuss upcoming studies.

Extracting Properties of Heavy-Light Mesons from LQCD
In this section, we will focus on short-distance corrections to the results obtained from lattice QCD. These corrections are common to both traditional techniques for extracting meson properties and to the multistate smearing method employed here. Explanation of the details of the multistate smearing method are deferred to the next section.
The corrections to the matrix elements of time component of the heavy-light axial current, J 05 (n), are computed by demanding that the ratio of the current renormalized with some continuum regularization scheme and the lattice-regularized current be unity. To be a little more precise, one computes this ratio, Z, using some matrix element, and from then on one multiplies any result obtained using the lattice-regularized current by Z. The states used to determine Z can be chosen for calculational convenience since the ratio is independent of the choice of states. Although there is no choice of states for which the numerator and denominator of the ratio are separately calculable, because the operators only differ at scales on the order of the cutoffs of the two regularizations, and at these scales QCD is perturbative, the ratio can be calculated in perturbation theory.
The procedure is actually slightly more complicated than explained in the preceding paragraph. Because the lattice QCD calculations are done in the heavy quark effective theory (a theory which does not have the same particle content as the full standard model), it is necessary to compute an additional ratio, Z cont . This is the ratio of the axial current renormalized at the scale m * b in the standard model to the axial current renormalized at a scale q * in the heavy quark effective theory.
The calculations of Z and Z cont are the subject of the following two subsections. An analysis of the heavy quark mass renormalization follows thereafter. In the last subsection, we summarize the various constants used in this study.

Calculation of Z
The ratio Z introduced above was calculated some time ago, but there is substantial uncertainty in these calculations. Tadpole-improved perturbation theory, as formulated by Lepage and Mackenzie [20], promises to reduce these uncertainties below the ten per cent level at one-loop. The application of the tadpole-improvement program to heavy quark effective theory has recently been discussed by Bernard [21], and calculations have been performed by Hernández and Hill [22]. In this subsection we summarize the calculation of Z within the framework of tadpole-improved perturbation theory. Hernández and Hill considered both the zero-separation and unit-separation point-split heavy-light axial currents. We will restrict our attention here to the zero-separation heavy-light current used in our Monte Carlo calculations. We will further restrict our attention to the case of Wilson fermions with r = 1. The use of tadpole-improved perturbation theory results in a substantial reduction in our best estimate of the central value and uncertainty for f B .
The root of the tadpole-improvement program is a nonperturbative renormalization of the basic operators in the lattice action. These redefinitions absorb the large renormalizations arising from lattice tadpole graphs. A related additional part of the Lepage-Mackenzie prescription is the use of a larger perturbative coupling. If one uses β to determine the perturbative coupling, α lat , one-loop perturbative corrections are consistently underestimated.
These perturbation theory problems are due to the fact that α lat is a poor choice of expansion parameter. For example at an inverse lattice spacing of 2 GeV, the tadpole-improved expansion parameter is α V = 0.16, which is twice as large as α lat . Lepage and Mackenzie argue that the best way to arrive at α V is from a non-perturbative lattice determination of a perturbatively calculable quantity, such as the gauge field plaquette expectation value.
Using tadpole-improvement of the Wilson action for quarks on the lattice as a guide, one can perform tadpole-improvement of the heavy quark action, and this has been done in Ref. [21]. Instead of discretizing where u 0 is defined as The combination U µ (x)/u 0 more closely corresponds to the continuum field (1 + iagA µ (x)), than does U µ (x) itself. With the tadpole-improved action, there is an additional factor of 1/u 0 for each gauge field link in the product. Thus the Green's function of two heavy-light currents separated by n 0 lattice spacings in the time direction satisfies, The B meson decay constant f B is usually extracted from numerical simulations by Thus we see that the tadpole improvement procedure has no effect on the fitted value of f B . Its only effect is the change that is, a linearly divergent mass renormalization.
So far we have seen that tadpole improvement does not affect the extraction of f B as it is generally done in lattice Monte Carlo calculations. However we must still take into account the effect of tadpole improvement of the light quark action, and this will involve some additional factors.
As conventionally defined in lattice Monte Carlo calculations the lattice operator J 0 5 involved in calculating f B is renormalized by a factor √ 2κ c Z, where κ c is the critical value of the hopping parameter for the light quarks. The tree level value of κ c is 1/8. Lepage and Mackenzie advocate a reorganization of perturbation theory such that a factor of √ 8κ c is included inZ and the renormalizing factor becomesZ/2.
Let us see what this factor does at one-loop. Calculations of Z have been carried out to one-loop order and the result is of the form, This definition of g(q) (and similar ones for h(q), j(q), and k(q) which will be introduced below) follow the definitions in Ref. [22]. A one-loop calculation of 8κ c has been performed, and the result is expressible as The relationship betweenZ and Z is The one-loop expression forZ is thereforẽ We continue with the application of the Lepage-Mackenzie prescription to determine the Λ-value of the coupling and the scale q * at which it is evaluated. The prescription for fixing the value of the coupling, α V is to extract it from a non-perturbative calculation of the 1 × 1 Wilson loop (i.e., the expectation value of the plaquette, U plaq ). Once the coupling is known at some scale (alternatively, once the value Λ V is known), it can be run to any other scale. The formula which relates α V to the non-perturbatively determined (lattice Monte Carlo) plaquette expectation value is, The coefficient of n f is the one appropriate for Wilson fermions with r = 1. In the quenched It remains to fix the scale q * at which α V (q) is evaluated in the expression forZ.
Lepage and Mackenzie propose to do that by calculating the expectation value of ln q 2 in the one-loop perturbative lattice correction. The formulae determining this scale are: Hernandez and Hill quote the errors on the numerically evaluated values of Y and Y ln < q 2 > as order 1 in the last decimal place.
Using the two-loop formula for the running of α V with zero quark flavors one obtains α V (q * ). The explicit dependence on the value of a −1 drops out of the ratio q * /Λ V . Hence, the only way lattice Monte Carlo results have been used so far is for the expectation value of the plaquette; the determination of a −1 has not yet entered. The results for Z at various values of β are summarized in Table 1.
This completes our calculation ofZ. It remains to multiplyZ by the continuum running and matching factor, Z cont . We leave this for the following subsection.

Continuum Running and Matching Factor
For consistency, the one loop computation ofZ should be combined with a two loop running in the continuum effective theory and a further one loop matching between the continuum effective and full theories. This produces the continuum running and matching factor Z cont which multipliesZ to give the full perturbative correction. It will turn out that the lattice to continuum matching factorZ is most significant, while Z cont produces only a small additional change.
In the previous section, the coupling α V was determined in a no-flavor (quenched) lattice theory. We must now match onto a four-flavor (or five-flavor depending on the value of q * ) continuum theory. In order for the infrared behavior of lattice and continuum theories to match exactly, we could use a continuum coupling whose value is equal to α V (q * ) (although differences between couplings are higher order effects). This implies we should choose a continuum scale q cont according to α cont (q cont ) = α V (q * ). In practice we ignore this criterion and simply set q cont = q * , using α V (q * ) everywhere in the matching. We then run in the continuum theory using four or five flavors depending on whether q * is greater or less than the b-quark threshold mass.
The b-quark threshold itself is determined as follows. We assume the pole mass m b pole is known and relate it to the MS running mass m b (µ) according to [23] .
, we fix the threshold mass m * b by solving We use the usual two loop result for α cont : Here n f is the number of light flavors and we take Λ = 175 MeV for five light flavors [24].
For other n f , Λ is fixed by demanding that α cont be continuous. Applying this procedure with m b pole = 4.72 GeV [25] gives Now that q * and m * b are fixed we use them in the combined two loop running plus one loop matching formula for Z cont . The one-loop anomalous dimension calculation [26] for the heavy-light axial current was extended to two loops by Ji and Musolf [27] and the two loop result was confirmed by Broadhurst and Grozin [29]. The same authors [30] also confirmed the one-loop matching calculation of Eichten and Hill [31]. The result for Z cont is: .
In this equation γ 0 and γ 1 come from the anomalous dimension γ of the heavy-light axial current in the effective theory, while c 1 comes from the effective to full theory matching at one loop in the continuum, obtained by the method described at the beginning of section 2.1. This matching produces a contribution to Z cont of Note that Ji and Musolf [27] quote c 1 with a sign error in the term which differs for vector  Table 1: Length scales and renormalization constants used in this paper. The values of a −1 for = 5.7, 5.9, and 6.1 are taken from Ref. [43]. For β = 6.3, a −1 is estimated from that at 6.1 by one-loop asymptotic freedom. Plaquette expectation values are taken from [20].Z is the lattice to continuum renormalization factor for the axial current, Z cont is the factor relating the continuum heavy quark effective theory to full QCD, and Z A is the overall renormalization factor used in previous discussions. Using equation (4) we determine the Z cont values given in Table 1. These are then combined with theZ's obtained in the previous subsection to give the overall perturbative corrections Z listed in Table 1. Here, to conform with our previous notation [28] we also list the quantity

Renormalization of the Heavy Quark Mass
Before leaving the subject of short-distance perturbative correction we will compute the tadpole-improved estimate of heavy quark mass renormalization.
In remains. The precise definition of the renormalized mass parameter in the dimensionally regularized heavy quark effective theory has been discussed by Falk, Neubert, and Luke [32].
Since the heavy quark mass renormalization is a linear divergence it vanishes in dimensional regularization with minimal subtraction.
On the lattice, the heavy quark mass is renormalized in the static limit. This mass shift, δm, is proportional to 1/a with a perturbatively calculable coefficient. Hence the mass, M ef f , calculated for the ground state B meson in this lattice theory can be expressed as in terms the continuum residual massΛ and the mass shift δm. It is clear from Eq.(6) that the mass M ef f is linearly divergent as the lattice spacing a → 0. However, it would appear that by measuring M ef f and removing the tadpole-improved one-loop mass counterterm, we have a determination ofΛ. Hence we could obtain m b as defined in Ref. [32] from m B . Unfortunately, Bigi et al [33] argue that non-perturbative effects ruin the preceding connection. Even if that is so, it is still possible to verify that the linearly divergent piece of M ef f is correctly determined by the perturbative calculation of δm.
At one loop, the lattice heavy quark mass renormalization is of the form, The one-loop correction to the self-energy is linearly divergent and positive. The coefficient, δm, of the mass counterterm is negative. The integrand k(q) is given by [34], As already noted the effect of tadpole improvement on the fitted value of the heavy quark mass is to change This non-perturbative reduction of the mass is accompanied by a reduction in the counterterm: wherek (q) = k(q) + j(q) As in the calculation of the scale q * , we need to compute the expectation vlaue of ln q 2 in the integral ofk(q). The formulae determining this scale are: The calculational technique used is exactly that of Ref. [22], and the errors on the numerically evaluated values of X and X ln < q 2 > are order 1 in the last decimal place. The values of α V (p * ) at the β values used in this study and a δm + ln u 0 are presented in Table 2.
Two methods for the determination ofΛ are possible. The first of these is to do the tadpole-improved perturbative subtraction just described. This results in a statistically independent result forΛ for each β value. The second method is to use the fact that the subtracted term is proportional to 1/a, and the physical value desired is independent of the lattice spacing. A two-parameter fit of the mass over the four values of β has two degrees of freedom, and one of the parameters is the one we desire. The latter procedure ignores the running of the coefficient of the 1/a term.

Comparison and Summary
It is worthwhile at this point to compare the results in Table 1 with the widely used value of Z A of 0.8, which does not benefit from tadpole improvement. Consider the results at β = 6.1. From the table, we findZZ cont / √ 8κ c to be 0.68. Consequently, tadpole-improved perturbation theory results in a reduction of the physical value of f B by a factor of 0.68/0.8, i.e., a reduction of 18%.
As an aside, we note that for the clover action-also termed "improved", but in the sense that order a effects rather than tadpole graphs are being incorporated-the axial current has been renormalized [35], but tadpole improvement has not been applied to this operator.
Now that perturbation theory has been reorganized to include tadpole corrections to all orders, we expect our one-loop calculation of the renormalization factor to be accurate to about 7%. This estimate of the magnitude of the two-loop corrections is obtained simply by squaring the largest one-loop correction for the various values of r quoted in Ref. [22], for both discretizations of the axial current considered there.
While the values for Z cont depend on the estimate of the lattice spacing, this dependence is weak; an increase by 10% in the estimate of the scale results in a reduction of at most 1.5% in the value of Z cont . Thus this source of uncertainty in the Z factor is negligible compared to both the estimated size of the two-loop corrections or the direct dependence of f B on the lattice spacing.
3 Analysis Procedure for Multistate Smearing

Relativistic Quark Model for Heavy-light Systems
The rapid deterioration [15,16] of the signal to noise quality of Euclidean correlators of heavy-light mesons at large Euclidean time makes the choice of an efficient smearing scheme essential if we wish to extract accurately the properties of low-lying heavy-light systems.
In the multistate smearing approach previously introduced [16], the coupling of smeared bilocal Coulomb gauge operators to higher meson states was reduced by using smearing wavefunctions obtained from a relativistic quark model (RQM). The basic features of such a model are (a detailed examination of the connection of such a model with the full field theory in the case of the t'Hooft model can be found in [36];see also, [37]): (a) the use of a relativistic kinetic term p 2 + µ 2 (with µ a constituent quark mass) for the kinetic piece of the Hamiltonian, and (b) a static confining potential V ( r), which can be chosen to be the static interaction energy obtained from correlators of temporal Wilson lines in lattice QCD.
The importance of relativistic kinematics in determining the shape of meson wavefunctions (with a light quark) was already implicit, if not clearly recognized, in the puzzling persistence of purely exponential falloff (≃ exp −Cr) of hadronic wavefunctions, instead of the more rapid falloff one might naively expect in a confining model (≃ exp −Cr 3/2 for a nonrelativistic particle in a linearly rising potential). This exponential falloff is due to the nonlocal character of the kinetic part of the RQM Hamiltonian which implies that Ψ n (r) satisfying cannot fall exponentially faster than e −µr (if it did, the integral over r ′ in (12) would be dominated by r ′ ≃ 0, giving an asymptotic behavior ≃ e −µr for the kinetic term, in contradiction with the assumed asymptotic behavior of (E n − V (r))Ψ(r) for any V (r) with power growth). In other words, irrespective of the power rise of the confining potential, relativistic kinematics automatically smears out the wavefunction of a light quark over the Compton wavelength corresponding to the constituent quark mass.
In the static limit in which the heavy quark mass is taken to infinity, the relativistic Schrodinger equation (12) gives a single parameter fit (the constituent quark mass µ is the only adjustable parameter after the static potential has been measured on the lattice) to a complete set of orthogonal spin-independent wavefunctions corresponding to arbitrary radial and orbital excitations of the heavy  The procedure used for generating lattice smearing functions from the RQM is as follows.
We wish to obtain orthonormal lattice wavefunctions which are eigenstates of a lattice RQM Hamiltonian defined on a L 3 lattice (with r, r ′ lattice sites): Such an eigenstate, in a channel of given orbital quantum number (S,P,D etc), will correspond to a pole of the resolvent applied to a source wavefunction Ψ (0) of the same orbital symmetry: For the starting source functions Ψ (0) , one may take for example a monopole localized at the origin for S-states, a dipole for P-states, etc. After the energy E is tuned close to an eigenvalue E n (until R is at least 10 3 larger than the background value), a smearing eigenstate Ψ In most cases we have found it adequate to fit the RQM constituent mass by matching the 1S wavefunction generated by the above procedure to the Coulomb gauge Bethe-Salpeter wavefunction obtained at a roughly fixed Euclidean time (corresponding to time slice 4 at β=5.9). For the particular case of β=5.9, κ=0.159, on a 16 3 lattice, a more detailed fitting procedure was used to determine the optimal choice of quark mass µ in order to fit the meson wavefunction at various time slices. The mean square deviation of Ψ (1S) smear from the measured quenched wavefunctions at various times T , for various µ, is displayed in Table 3. Table 3: Mean square deviation of RQM smearing and LQCD heavy-light wavefunctions (x10 6 ), for β = 5.9, κ = .159 as a function of time. µ is the constituent quark mass parameter in the RQM Hamiltonian. here reduces to a very small level (4 percent or less) the coupling of the exact ground state to the higher smearing states, but will not turn out to be essential to the extraction of accurate masses and couplings for the ground state. A detailed discussion of the dependence of the results on the RQM mass parameter chosen for smearing is given in Section 5.1.
As described in a recent article [17], the RQM gives a single-parameter fit to all the excited radial and orbital meson wavefunctions of our heavy-light system. After fixing µ by a match to the 1S wavefunction, we have found remarkable agreement with measured excited state wavefunctions (for example, the 1P state, cf [17]). This agreement suggests that this Ansatz accurately describes at least the valence quark sector of the full mesonic bound-state.

Multistate Smearing
Our object in this section is to outline a general procedure for extracting the maximum usable information from the multistate correlator matrix: where q(Q) are light (heavy) quark operators, and the Ψ Defining states we have where the states | n > are exact eigenstates of the lattice Coulomb gauge transfer matrix with eigenvalues e −En . The remainder term of order e −EM+1T will of course be small at large Euclidean time, but in addition should have a small prefactor to the extent that our smearing functions Ψ Next, define mixing coefficients (in our case, they are real): Neglecting the exponential contamination of order e −EM+1T , we see that the multistate coupling matrix can be fit to an expression of the form Of course, we cannot hope to extract M independent time-dependencies with N < M smearing wavefunctions, so only N ≥ M will be considered. Typically we shall extract the maximum information from the lattice data by picking N = M + 1 (the extra source function being the local current needed for the extraction of f B ).
The fit is performed by a chi-square minimization of with respect to the fitting parameters {v a n , E n }, over a fitting range T < ≤ T ≤ T > in Euclidean time. The fit is performed on an ensemble of N c jack-knife coupling matrices obtained by replacing each in turn of the N c coupling matrices by the average matrix and reaveraging. We have chosen | Φ N >≡ J axial (0) | 0 >, so the parameters v N n should be interpreted as lattice f-parameters for the ground and excited meson states, E n as the corresponding masses, and v a n (a=1,2,..,M ) as mixing coefficients indicating the degree of overlap of the exact meson states with our RQM smeared states | Φ a >. The sum over a, b in (22) does not include the local-local correlator a = b = N , which is not well described by a sum over a few low-lying states. Note that this fitting procedure automatically gives the lattice f-parameters without the need to divide by the square-root of the smearedsmeared correlator as in the usual approach. Moreover, the ensemble of N c parameter sets {v a n , E n } obtained in this way can be subjected to a straightforward statistical analysis to determine the error in each of these parameters separately, correlations between parameters (e.g. between masses and f-parameters), and so on.
Once the overlaps < Φ a | n > have been estimated by a best fit of C ab (T ), a smearing operator can be constructed which is guaranteed to contain at most one of the first N exact meson states, thereby removing any other exponential time-dependence to the e −EN+1T level. Specifically, if ǫ a1a2..aN is the totally antisymmetric symbol in N dimensions, the smeared state is guaranteed (to the extent that we have accurately extracted the mixing coefficients v a n ) to contain only the exact meson state | A >, together with contaminations from the (N + 1)'th excited state and higher. An effective mass plot of the usual kind can then be obtained for the A'th state by displaying (we use a smeared-local correlator to minimize noise) where Of course, this plot will be most flat for the ground state A = 1, where the relative expo- For each β the time window for the multistate fit was chosen to be over approximately the same interval in physical units, viz. about 1 3 f m to 1f m. Noting that the splitting between the ground state and the second excited state is found, in our multistate calculations, to be around 800 to 900M eV , the choice of 1 3 f m for the lower end of the time window should provide an exponential supression of excited states by at least a factor of 3. In addition, our careful tuning of the smearing functions should produce a relatively small coefficient for the higher excited states. The equality of the smeared-smeared and smeared-local effective masses exhibited in the plots confirms our choice of fitting interval. A more complete discussion of systematic errors due to excited states is given in Section 5.1.
The fitting formula (22) is easily generalized to allow a global fit to the data at various κ values (for fixed β): this is essential in order to take into account correlations between the coupling matrices at different κ values, which would affect our estimate for the error of the results when linearly extrapolated to κ c . Namely, (22) is replaced by

Lattice Results for Heavy-Light Mesons
To extract results for masses and decay constants we have used the set of gauge configurations and light quark propagators enumerated in Table 4. The light quark action we use is not O(a) improved. The four columns in Table 4   The multistate smearing analysis outlined in Section 3 provides a powerful method for extracting heavy-light meson parameters. Unlike single or double exponential fits to singlechannel "smeared-smeared" and "smeared-local" correlators, the fitting of the N ×N matrix of correlators to an expression of the form (21) is highly constrained. As we will show, this method allows a determination of f B and other heavy-light parameters which is less prone to systematic errors than previously applied methods. In this section, we present our results.
Most of the results presented here were obtained from the fitting procedure discussed in In the multistate fitting procedure, E 1 is the energy associated with the ground state contribution to the correlator. This parameter represents a combination of the binding energy of the B-meson ground state plus a divergent mass shift of the heavy quark. Recall that only the bare mass of the heavy quark is removed in reducing to the effective static theory. The mass shift induced by QCD is therefore measureable on the lattice as a non-scaling piece in the parameter E 1 . We measure the dependence of E 1 on both the light-quark hopping parameter κ and the lattice spacing a. The graph in Fig. 4 shows the κ-dependence of E 1 for the four values of β studied. In each of the four data sets, the dependence on κ −1 is quite linear, allowing an accurate extrapolation to the chiral limit κ c . Numerical results are tabulated in Table 5. and aE 1 at β = 5.9 listed in Table 5 are the infinite volume values obtained by fitting all three box sizes to a Luscher finite volume formula, as described in Section 5.2 (except for κ = .156, which has only been done on the 16 3 box). The results for the three box sizes are listed separately in Tables 7 and 8. A measure of the overlap between the true ground state and the RQM wavefunction smeared operators is given by The measured values of this overlap for each of the fits is recorded in Table 5.
In addition to analyzing the data at each value of κ separately, we have also performed a simultaneous fit to all κ values for a given β by allowing the multistate fitting parameters v a n to depend on κ, as described in Section 3. Using this procedure on jackknifed subensembles ). Results at κ = κ c are from the multi-κ fits as discussed in Sec. 3.2. The column labeled overlap is a measure of the total overlap between the true ground state and the RQM wavefunction smeared operators used.
Let us consider E 1 as a function of the lattice spacing a and of the naive light quark mass (see Appendix). The linear dependence on κ −1 becomes The quantity M ef f defined in Section 2, Eq. (6) is just E 1 evaluated at κ = κ c , The slope parameters E ′ 1 (0, a) are obtained by a correlated fit to all κ-values, as described in Section 3. The results for the slopes and intercepts are given in Table 6. Table 6: Slopes and intercepts for aE 1 andf B as a function of κ −1 − κ −1 c . For a given β,values are obtained from a simultaneous two-state fit to all values of κ over a time window ∆T . (11) .236 (25) .564 (29) .424(72) 5.9 3-10 .638 (9) .350 (25) .250 (13) .444 (38) 6.1 4-12 .544 (12) .450 (47) .138 (13) .387 (46) 6.3 5-14 .499 (9) .376(56) .099 (7) .290 (56) In a similar way, we obtain a linear fit (see Fig. 4) to the κ −1 dependence of the groundstate pseudoscalar decay constant f B . Define a quantityf B which is just the matrix element parameter v N 1 in the multistate fit (overlap of the ground state with the δ-function source), with a normalization factor √ 2κ for the light quark included. In the scaling limit, the physical value of f B is related tof B by the following multiplicative constants, where Z A is the renormalization factor associated with matching the full theory with the effective static theory on the lattice, as discussed in [34], [42], and in Section 2.
The lattice spacing dependence of E 1 in the chiral limit m q = 0 is plotted in Fig. 4.
The results are consistent with a linear a dependence, with The first term E 1 (0, 0) is the linearly divergent (i.e. O(1/a)) term in the heavy quark mass shift.
The tadpole-improved estimate of aδm was discussed in Section 2. The corresponding quantity is a δm + ln u 0 , and was tabulated in Table 2. Comparing the ground state effective mass with the tadpole improved 1-loop result at each lattice spacing one finds about a 30 % discrepancy in the singular part of the mass shift, which can easily be accounted for by higher-loop and/or nonperturbative contributions. In fact, ordinary (non-improved) one-loop perturbation theory [34] gives Simply identifying this value with the extrapolated lattice result gives α s ≃ 0.162(6) which is in reasonable agreement with other determinations of α s in the range of lattice spacings considered here [20]. So the entire "discrepancy" can be removed by a reasonable redefinition of the perturbative coupling being used.
Since the divergent self-mass of the heavy quark is independent of light quark mass, it will cancel in the mass difference ∆M Bs , and the latter should therefore scale properly with a.
For each value of β, we use the determination of κ s and κ u discussed in the Appendix, along with the observed κ dependence of E 1 to determine ∆M Bs . The results are shown in Fig.   4.
A linear extrapolation of the mass difference ∆M Bs to a = 0 gives Notice that the results for ∆M Bs shown in Fig. 4 Table 6. The ratio f Bs /f Bu is plotted in Fig. 4. Again, the a-dependence of the ratio is much weaker than that of f B itself. Extrapolating to a = 0, we obtain It is worth emphasizing here that our ability to control excited state contamination depends crucially on the use of the multi-state fitting procedure. Although the 1-state fits were all found to lead to consistent results after extrapolation to t = ∞, an accurate extrapolation would not have been possible without an independent determination of the splitting ∆, which is only obtainable from the 2-state fit.

Finite volume corrections
Using the scales in Table 1 due to finite volume effects and describe our method for estimating these effects. Although our overall conclusion is that these effects are negligible on the lattices considered here, the estimates discussed in this section may be useful for selecting β's and lattice sizes in subsequent studies.
An extensive theoretical study of finite volume effects on field theoretic calculations has been carried out by Luscher [45]. Consider, for example, the effect on the mass of a particle m P . For large enough volume, the leading effect is due to the propagation of the lightest mass meson (e.g. pion) "around the world," leading to the expression for m(L), the particle mass in an L × L × L box, where the exponent λ is determined by the mass of the pion, and A is given in terms of the on-shell πP P coupling. This finite volume correction can be interpreted as the effect of squeezing the pion cloud surrounding the particle. A somewhat different situation takes place when the particle P is a loosely bound state of constituents. In this case, the finite size effect is caused by the squeezing of the bound state wave function [46]. As pointed out by Luscher [47], this situation falls into the same general framework as that which led to Eq.(39), except that, in this case, the particle that travels around the world is one of the constituents of the bound state. In fact, for a nonrelativistic bound state in a non-confining potential, the finite volume effect assumes exactly the same form as (39), but in this case, the exponent λ is related to the binding energy (and hence to the spatial extent of the bound state wave function).  Tables 7 and 8. They are seen to be the same, within errors, in all three size boxes, and thus, no significant finite size effect is observed. In order to determine an upper limit on the finite volume corrections to our results, we will make the assumption that these effects can be parametrized in the Luscher form (39).   (20) .154 .716 (9) .718 (7) .719 (6) [   (20) .154 .341 (21) .347 (13) .346 (14) .157 .299 (22) .303 (14) .298 (15) .158 .284 (24) .288 (15) .278 (16) .159 .279 (35) .271 (18) .252 (17) κ c .261 (29) .260 (18) .245 (18) calculated on 12 3 , 16 3 , and 20 3 boxes. (The relativistic Van Royen-Weisskopf formula [37] was used to obtain f B from the RQM wave function.) For the lightest quark mass studied, to the value of λ, we also fit the data using λ/a = 0.9 GeV, which yields an upper bound on the 16 3 box about 50% smaller than the λ/a = 0.7 GeV fit. In Table 9, we give the results of fitting the 12 3 , 16 3 , and 20 3 LQCD Monte Carlo results to the finite volume formula Eq. From Table 9 it is seen that, in all cases, the estimated finite volume effect on both aE 1 and on f B is smaller than our statistical error by more than a factor of two on the 16 3 lattice and by more than an order of magnitude on the 20 3 lattice at β = 5.9. To determine the size of finite volume effects on the quantities M Bs − M B d and f Bs /f Bu , we also need to estimate the error on the slope parameters E ′ 1 (0) and f ′ B in Table 6. From the κ-dependence of the finite volume fit parameters, we estimate an approximate upper bound on the finite volume error for the slopes to be (∆E ′ 1 )/E ′ 1 < .03 and (∆f ′ B )/f ′ B < .05 for the 16 3 box, and by .719 (5) .678 (6) .664 (7) .645 (9) .638 (9) [.719 (7) .346 (11) .300 (11) .283 (12) .257 (14) .246 (14) [.348 (16) Recently, it has been argued [48] that, in intermediate ranges of volume where the asymptotic behavior predicted by Luscher's volume formula has not yet set in, the volume dependence might be expected to exhibit a power law dependence of the form instead of the exponential falloff of Luscher's result. This power law form is also found by the authors of Ref. [48] to fit better to their data on light hadron masses (in full QCD). If we assume a similar power-law dependence for the heavy-light data, we obtain extrapolated infinite volume results and bounds on finite volume corrections which differ from those obtained with Luscher's form. The values in Table 9 which are enclosed in square brackets are the results obtained by assuming a power law dependence of the form (42). Notice  Table 10. The range of kappa values used in the fit is indicated in the first column using the abbreviated notation κ = 0.15x → x (thus: 467 indicates that the extrapolation to κ c was made using κ values 0.154, 0.156, and 0.157).
Referring to Table 10 .643 (10) .319 (22) .464 (29 We may conclude from the preceding that, as in the case of finite volume corrections, nonlinearities in the chiral extrapolation are not an important source of systematic error in our results. In order to arrive at an actual estimate of the chiral extrapolation contribution to the total systematic error we have taken the variation in the 3-kappa fits in Table(10) (i.e rows labelled 467,678, and 789) which have a reasonable lever arm in 1/κ, and fairly small statistical errors, as an indication of the extrapolation error to κ c (i.e. in f Bu and M Bu ).
As we measure quite close to κ s (at κ =0.157), there is effectively no extrapolation error in the strange quark quantities. We assume that the chiral extrapolation at β=5.9 is typical of other β values. In this way a chiral extrapolation part of the total systematic errors quoted in Section 6.1 can be obtained.

Scale Errors
In order to quote physical values of masses and decay constants, one must select a particular dimensionful quantity to define the scale. In our discussion, we have taken the values of a −1 obtained from the 1P-1S charmonium splitting [43] at β = 5.7, 5.9, and 6.1. Our choice of a −1 at β = 6.3 is obtained by evolving from β = 6.1 via one-loop asymptotic freedom.
(The same value of a −1 = 3.08 is also obtained from our value for m ρ .) Other possible choices for the scale-defining parameter include string tension, rho mass, and f π . Since the quoted values of the decay constant f B include a factor of a 3 2 , it is particularly important to estimate the possible systematic error in our results arising from uncertainty in the overall scale at each β value. In Fig. 5.4 we have plotted the scale obtained from m ρ (circles and filled circles), f π (squares), and string tension (diamonds and filled diamonds) relative to the scales chosen in this paper (Table 1). Also included on the plot are points (filled squares) obtained from lattice calculations of the deconfinement temperature T c [38]. Since the experimental value of T c is not known, these calculations only give a relative determination of the scale at different β's. (The absolute scale for these points has been chosen to be equal to that in Table 1 at β = 6.3, which corresponds to a deconfinement temperature of kT c = 264 MeV.) The values for m ρ are from GF11 [40] (filled circles) and from our data (Table A)(circles), while the string tension is a combination of our results at 5.7, 5.9, and 6.1 (diamonds) and those of Ref.
[39] at 6.0, 6.2, and 6.4 (filled diamonds). For the latter points, the charmonium scales were estimated from Table 1 by linear interpolation in ln a.
The values for f π are taken from Table A. The trend exhibited by the data in Fig. 5.4 indicates a significant scale discrepancy in the range β = 5.7 to 6.0, with the m ρ scale being about 10-20% higher and the stringtension scale about 10-15% lower than charmonium. For β ≥ 6.2, the scales appear to converge to much better agreement, with deviations of ≤ 5%. This suggests that much of the discrepancy at lower β is due to finite lattice spacing effects, as opposed to being an effect of the quenched approximation. (Discrepancies which do not go away in the scaling limit can be ascribed to the neglect of closed quark loops.) The data shown in Fig. 1 illustrates are therefore subject to a 5% and 7% error (resp.), assuming that continuum extrapolated objects are determined primarily by the larger β values where the scale discrepancy is small.
The larger error for f B arises from the fact that the quantity computed on the lattice scales like a 3/2 . These issues are particularly important in the case of heavy meson decay constants. We find that the lattice spacing dependence for f B in the static limit is considerably stronger than for f π . This is illustrated in Fig. 5.5. Although our calculations have led to reasonably precise results at finite a, quoting a systematic error on the continuum extrapolated result clearly requires an investigation of the variation induced by alternative fitting procedures. We have investigated this variation by fitting the f B results at finite β to various functional forms for the finite a correction. In Table 11 Table 11 as an estimate of the systematic error associated with extrapolating to a = 0. In our final result for f B , we quote the continuum value obtained from the linear fit with asymmetric errors of +26 and -0 to allow for this extrapolation uncertainty.

Extrapolation to the continuum
6 Physical Results and Discussion

Final results
In this final section, we collect our results for f Bs /f Bu , M Bs − M Bu , and f Bu at each lattice spacing calculated and then present the results extrapolated to the continuum. The summary of results for each of the four lattice spacings studied are collected in Table 12. A comparison with other recently reported results follows immediately. In all cases the first error quoted is a statistical error obtained by the jackknife procedure described previously.  (7)) is very small, so we have not included a continuum extrapolation error. The result is For the mass splitting M Bs − M Bu the systematic errors include finite volume effects, the chiral extrapolation (to determine M Bu ), and an estimate of the scale error. For the continuum extrapolation, we quote the result obtained from the linear fit in a, and take the difference between the two fits shown in Table 11 as an estimate of our extrapolation error (+0 to −6). Combining the extrapolation error with our other systematic errors, we obtain the result Finally, for f B itself, there is, in addition to the usual finite volume and chiral extrapolation errors, a substantial (not necessarily linear) lattice spacing dependence, as well as the possibility of sizable renormalization corrections to Z at the two and higher loop level.
These additional systematic errors are quoted separately in our final result:

Comparison with other results
There have been several recently reported studies of the heavy-light meson system in the static limit of quenched LQCD. In this Section we will assess the results reported here in comparison with these other studies. We will consider four quantities, E 1 , the heavy-light First consider the ground state energy E 1 . In Section 4.1 it was shown that the κ and a dependence of E 1 is most easily described by introducing the naive quark mass The values of κ c for β = 5.74 and 6.26 are taken from [14] while that at β = 6.0 is taken from [12]. The κ c values for β = 5.7, 5.9, 6.1, and 6.3 are from our own data. The values for the scale a −1 are subject to somewhat more uncertainty.
For 5.7, 5.9, 6.1, and 6.3, the scales we have used are given in Table 1. In order to have a reasonably self-consistent set of scales, the remaining values for a −1 are obtained from those at 5.7, 5.9, and 6.1 by a simple linear interpolation or extrapolation in ln a. Note that we are only using the value of a −1 here to define the naive quark mass m q . The numerical values of E 1 andf B are obtained directly from the lattice propagators and do not depend on choice of scale. Thus, our comparison of data is insensitive to an overall, uniform change of scales. The values for the ground state energy E 1 (in lattice units), extrapolated to κ c , from the various studies are plotted in Fig.6.2. We conclude from these plots that all the data are in reasonable agreement, both on the magnitude of E 1 and on its a dependence. This is not surprising since, for any reasonably well-chosen smearing function, the value of However, the results are basically consistent within errors. For this quantity, the lattice spacing dependence is very mild.
Results from various groups for the ratio of decay constants f Bs /f Bu are compared with ours in Fig 6.2.
In contrast with the situation for E 1 or M Bs − M Bu , there is significant disagreement on the value off B among the various studies, even after correcting for the different values of κ and β. Results from references [44,41,12,14,49] are compared with our data in figure 6.2. This is a plot off B a −3/2 against lattice spacing a, so that we can compare 'raw' lattice results without renormalisation factors. The results of [12] (at β = 6.3) and [49] are in good agreement with ours, while [44,41,14] report substantially larger values.
We believe that the discrepancies can be traced primarily to the incomplete isolation of the ground state with the sources chosen. Observe in particular that the result in [49], where a variational method is used to isolate the ground state, is in good agreement with us. Sensitivity of the apparent value of f B to the form of the source has been emphasized in [50].

Future Studies
We have presented results for the decay constant f B and for masses of low-lying heavy-light states in the static approximation. The analysis procedure introduces several improvements over previous smearing methods. First, the success of the RQM in reproducing the measured lattice wave functions is exploited by using the RQM to construct not only an accurate ground state wave function, but also a set of orthonormal excited state smearing functions.
Second, we make full use of the information contained in the matrix of smeared-smeared and smeared-local correlators, including both ground state and excited state smearing functions at each end. Our method provides much greater control over systematic errors from higher state contamination, because of the fact that the source smearing functions are tuned directly to the lattice wave functions, without regard to the behavior of the effective mass plots. The appearance of long plateaus in the SS and SL plots at the same value of effective mass is thus strong evidence that the systematic error from higher states has been largely eliminated.
We will report the application of our methods for heavy-light mesons to the spectrum of radial and orbital excitations for heavy-light systems and the B parameter in forthcoming papers.
We expect to be able to improve the accuracy of the present results for the M Bs − M Bu and f Bs /f Bu by using larger ensembles. There are two other sources of systematic uncertainty in our results for f B . Use of the Wilson action for the light quarks implies lattice spacing corrections in O(a) and the large one loop renormalization for the axial current suggests that the two loop correction may be sizable. Study of these effects will be required to substantially improve the error on f B .

A Light Hadron Results
In this Appendix we report the results for the light quark sector which we have used to set the parameters for the heavy-light meson properties. First consider the light-light pseudoscalar meson mass as a function of the light quark hopping parameters. This serves to determine the critical hopping parameter at each value of β studied. Using the physical K mass, we also establish the correct hopping parameter for the strange quark. (The small effect from the nonzero up or down quark mass has also been included in our results.) In addition to determining the light quark hopping parameters from the pseudoscalar meson masses, the meson propagators have also been analyzed to obtain the values of f π and m ρ . All of the results discussed here for the light hadron parameters have been extracted from local qq operators. An analysis of light meson parameters using smeared sources in Coulomb gauge is in progress and will be reported elsewhere.
Let the hopping parameters for the two valence quarks in the light meson be denoted by κ 1 and κ 2 , and denote the pseudoscalar meson mass by m P . All of our data is consistent with a linear dependence of m 2 P on κ −1 1 and κ −1 2 . The results can be expressed in terms of the parameters C and κ c of the linear fit (m P a) 2 = C κ −1 The results for the four β values studied are given in Table 13. Using the scales in Table 1, the naive quark masses, defined by are also listed in Table 13 Table 13: Light quark hopping parameters obtained from light pseudoscalar meson mass measurements. Values are shown for each run, as well as the infinite volume extrapolated value at β = 5.9. C is the slope of the pseudoscalar (mass) 2 as a function of κ −1 . The last two columns give the strange and up quark masses ×2a. Errors are statistical.
The uncertainty arising from the scale determination, along with other sources of error, will be discussed in Section 6.

B Multistate Extraction of Meson Wavefunctions
A particularly graphic illustration of the power of the smearing technique in reducing the contamination of higher states is obtained by examining the time development of the Coulomb gauge Bethe-Salpeter wavefunction of a static-light meson beginning with either (a) the RQM smeared source defined in equation (17), or (b) a cube smeared source [41].
In the first case, we extract is shown for a source smeared over a cube of width 7 lattice spacings. Here the convergence is much slower, with the pointwise convergence near the origin particularly tardy.
Although the use of a single smearing function obtained from the RQM is adequate to the task of extracting the ground-state Bethe-Salpeter wavefunction, even the improved smearing given by the RQM is not sufficient if we wish to do the same for the higher excited states in a given channel. In Section 4 we showed how to define optimized smeared states |Φ A > in which the admixture of all but one (the A'th) of the first M meson states in a given channel is tuned to zero. For example, taking M =3 at β=5.9, κ=0.159 on a 16 3 lattice, one finds that the choice of smearing function (cf. Sect 3.2,Eq (23)) produces an optimized first excited state in the sense that admixtures of the ground and second excited state are tuned out (based on a fit of the form (22)) with T < =2,T > =7). The use of such an optimized smearing is crucial if we wish to extract the correct Bethe-Salpeter wavefunction of the first radial excited state near the origin. Any sizable admixture of the exact meson ground-state will otherwise dominate the small r region of the wavefunction at large Euclidean time, before the higher (2nd, 3rd,etc) states have decayed away. In Figure  B we show this phenomenon with a Bethe-Salpeter wavefunction Ψ BS (R, T ) defined as Ψ BS (R, T ) ≡< 0 | q( R, T )Q(0, T ) | Φ (2) , 0 > The wavefunction (renormalized to unit norm) is plotted for Euclidean times T =1,2,4, and 6. There is a steady upward drift of the wavefunction at the origin as T increase to 4, but by time slice 6 the influence of the ground state is clearly apparent as the latter begins to dominate the evolved meson state. It would clearly be very difficult to draw any firm conclusions about the behavior of the excited state wavefunction close to the origin from these measurements.
On the other hand, using the optimized smearing found above (55), and computing Ψ opt BS (R, T ) ≡< 0 | q( R, T )Q(0, T ) |Φ (2) , 0 > (57) one finds (Fig.B shows the optimized wavefunction for times T =1,2,3) a rapid convergence to an asymptotic shape by the third time slice, giving a value for the wavefunction at the origin ≃0.44, as compared to a maximal value ≃0.36 obtained at T =4 from Eq (56) before convergence is lost. Eventually, of course, the ground state will dominate in this case also, but by using the optimized state, we correctly extract the exact Bethe-Salpeter wavefunction for the excited state before the ground state contamination has a chance to become sizable.        Error bars are statistical only. Solid lines are obtained from a simultaneous two-state fit to all kappa values for a given β as described in the text.        Table 1. Our data is denoted by open symbols. Figure 16: f B at κ = κ c as a function of lattice spacing. The scales and renormalization constants used are given in Table 1. The solid line is the best linear fit. The dashed line is a quadratic (a 2 ) fit used to estimate systematic error in a → 0 extrapolation.
(See Section 5.5).   Points at the same a have been slightly displaced for readability. Points at the same a have been slightly displaced for readability. Figure 21: m 2 π vs. naive quark mass for β = 5.7, 5.9, 6.1, and 6.3. The vertical scales for 5.9, 6.1, and 6.3 are offset by multiples of 0.2 for display purposes. Points labeled 5.9e and 5.9u represent mesons with equal and unequal quark masses, respectively.
Values with unequal quark mass are plotted at the average mass (m 1 + m 2 )/2.