HQET at order $1/m$: II. Spectroscopy in the quenched approximation

Using Heavy Quark Effective Theory with non-perturbatively determined parameters in a quenched lattice calculation, we evaluate the splittings between the ground state and the first two radially excited states of the $B_s$ system at static order. We also determine the splitting between first excited and ground state, and between the $B_s^*$ and $B_s$ ground states to order $1/m_b$. The Generalized Eigenvalue Problem and the use of all-to-all propagators are important ingredients of our approach.


Introduction
In recent years, there has been significant progress in determining the spectrum of hadrons containing a b quark, both experimentally [1][2][3][4][5], and on the lattice [6][7][8][9][10]. A comparison of theory and experiment is of considerable interest for these hydrogen like systems, in particular since Heavy Quark Effective Theory (HQET) [11][12][13][14] is applicable and is at the same time an important theoretical tool to isolate new physics in the flavor sector [15].
The extraction of information on excited states from lattice simulations is a difficult problem, since the excited states appear as subleading exponentially decaying contributions to the lattice correlators, which at large times are dominated by the ground state and at small times by the combined contributions of arbitrarily highly excited states, and on top of this are affected by noise. A number of different methods [16][17][18][19] have been proposed for overcoming this challenge; recently we have shown that with an efficient use of the generalized eigenvalue problem (GEVP), rigorous statements about the systematic error caused by the admixture of other states can be proven [20]. In particular, we have shown that for a suitable choice of Euclidean times t and t 0 , the systematic error decays exponentially with an exponent given by an energy gap that can be made large by an appropriately chosen variational basis.
Since their Compton wavelength is much shorter than any realistically achievable lattice spacing, b quarks cannot be simulated as relativistic quarks on a lattice. A description of heavy-light mesons that is suitable for use in the context of lattice QCD simulations is given by HQET with non-perturbatively determined parameters [21,22]. The parameters necessary to match HQET to QCD at order O(1/m b ) have been determined in the quenched approximation by our collaboration [23], and the non-perturbative determination for QCD with N f = 2 dynamical quark flavors is far advanced [24].
In this paper, we present our results for the heavy-light meson spectrum from HQET using the GEVP method. In section 2, we give a brief review of our methods. Details of the simulations and data analysis procedures are given in section 3, and our results for the quenched heavy-light spectrum are presented. Section 4 contains our conclusions. Some technical details of our implementation of all-to-all propagators are relegated to the appendix.

Non-perturbative HQET
HQET on the lattice offers a theoretically rigorous approach to the physics of B-mesons since it is based on a strict expansion of QCD correlation functions in powers of 1/m b around the limit m b → ∞. Subleading effects are described by insertions of higher dimensional operators whose coupling constants are formally O(1/m b ) to the appropriate power. This means that HQET can be renormalized and matched to QCD in a completely non-perturbative way [25], implying the existence of the continuum limit at any fixed order in the 1/m b expansion.
To fix the notation we write the HQET action at O(1/m b ) as where ψ h satisfies 1+γ 0 2 ψ h = ψ h . The parameters ω kin and ω spin are of order 1/m b , and δm is the counter-term absorbing the power-divergences of the static quark self energy. The signal-to-noise ratio of large-distance correlation functions is significantly improved by replacing the link U (x, 0) in the backward covariant derivative [26].
Exponentiating the action of eq. (2.1) would give (non-renormalizable) NRQCD [27]; in order to retain the renormalizability of the static theory, we treat the theory in a strict expansion in 1/m b , where the O(1/m b ) parts of the action appear as insertions in correlation functions. For the expectation value of some operator O this means (ignoring the possibility of explicit 1/m b operator corrections, which do not affect the energy levels [23]) where O stat denotes the expectation value in the static approximation.
To fully specify HQET, the parameters δm, ω kin and ω spin must be determined by matching to QCD. In order to retain the asymptotic convergence in 1/m b this matching must be done non-perturbatively, since for a perturbative matching at loop order l, the O(g 2l ) truncation error of the static term is much larger than the power corrections to the static limit: A fully non-perturbative determination of the parameters of HQET has been carried out in [23]. Here we employ the same discretization of QCD and HQET and in particular the determined values of ω kin and ω spin . For further details of the matching and discretization, the reader is referred to [23].

The Generalized Eigenvalue Problem
In this section, we recall the relevant contents of [20]. Starting from some fields O i (x) localised on a time slice and their momentum zero projection a 3 x O i (x) =Õ i (x 0 ), a matrix of Euclidean space correlation functions, provides the basis for the GEVP Effective energies for the n-th energy level are given by where a is the lattice spacing. Provided that t 0 > t/2, the effective energies converge to the exact energy levels as [20] E eff In HQET, all correlation functions are computed in an expansion in a small parameter, ω ∝ 1/m b . Correspondingly, the energy levels expand as where x ∈ {kin, spin}, with the behavior at large time t ≤ 2t 0 , 14) with (v stat m (t, t 0 ) , C stat (t 0 ) v stat n (t, t 0 )) = δ mn . We note that the GEVP is only ever solved for the static correlator matrices.

Mass splittings to O(1/m b )
With the static Lagrangian eq. (2.2), all HQET energies satisfy exactly E n = E n | δm=0 + 1 a log(1 + aδm), and the power divergent δm drops out in energy differences. Since we only consider these in this paper, we need just ω kin and ω spin of [23].
The excitation energies at the static order of HQET are given simply by the differences of the static energies, and at order 1/m b , the excitation energies of pseudoscalar B s states become At the static order, the masses of pseudoscalar and vector states are degenerate due to spin symmetry [13]. This degeneracy is lifted

Lattice parameters
We use three quenched ensembles of 100 configurations each, generated using the Wilson gauge action with parameters β = 6.0219, 6.2885 and 6.4956. The physical volume was kept constant at L ≈ 1.5 fm, leading to L/a = 16, 24, 32 for the three ensembles. We used time extent T = 2L throughout. The static quark is discretized on each ensemble with both the HYP1 and HYP2 [26,28,29] actions. The light valence quark is discretized by a non-perturbatively O(a)improved Wilson action [30,31], and its mass was fixed to the strange quark mass, giving κ s = 0.133849, 0.1349798, 0.1350299, respectively [32]. A summary of the simulation parameters used, with the corresponding values of ω kin and ω spin , is given in Table 1.

Measurements of correlation functions
The strange quark propagators are computed through a variant of the Dublin all-to-all method [34]. We use approximate instead of exact low modes (the method remains exact) and employ even-odd preconditioning in order to reduce the size of the stochastically estimated inverse of the Dirac operator by a factor of 2. The reader is referred to appendix A for details.
The interpolating fields are constructed using quark bilinears where the gauge fields in the covariant Laplacian ∆ are first smeared with 3 iterations of (spatial) APE smearing [36,37] to reduce noise. At β = 6.2885, we use R k = 0, 22, 45, 67, 90, 135, 180, 225 with κ G = 0.1. At the other values of β, we rescale the values of R k used so that the physical size r phys,k ≈ 2a √ κ G R k of the wavefunctions is kept fixed; in particular we keep r max = r phys,7 ≈ 0.6 fm.
For these bilinears, we compute the following correlators:

Determination of energies
The correlator matrices of eq. An alternative procedure, used in [20] is "pruning" [38], where the N × N matrices are formed by projection onto the subspaces spanned by the lowest N eigenvectors of C(t i ) for some small t i . We decided not to use this version since it introduces a dependence on the relative normalization of the fields O i . Moreover, with thinning we found a somewhat faster convergence to the plateau for the ground state energy compared to the pruning version [39] when the normalization of O i in [20] is used. We note that the fact that we found the same low-lying energy levels with thinning as with pruning is a n 2 3 4 5 6 r 0 (E stat n − E stat 1 ) 1.50(5) 2.7(1) 4.0 5.0 6.0 Table 2: Rough estimate of the static spectrum. For n ≤ 3, the gaps come from the continuum limit; for n ≥ 4, the rough estimate used to stabilize the fit is quoted.
further confirmation that the GEVP is quite robust against changes of the variational basis employed. For each of the resulting N × N correlator matrices, we solve the static GEVP and compute the static and O(1/m b ) energies. This gives a series of estimates E eff,stat n (N, t, t 0 ), E eff,kin n (N, t, t 0 ) and E eff,spin n (N, t, t 0 ) with associated statistical errors, which we determine by a full Jackknife analysis.
To arrive at final numbers for E n we also need to estimate the size of the systematic errors coming from the higher excited states. To do this, we first perform a fit of the form to the GEVP results for E eff,stat n (N, t, t 0 ), fitting the data at N = 3, . . . , 5, 1 2 r 0 < t 0 ≤ 6a, t 0 ≤ t ≤ 2t 0 and n = 1, . . . , 6 simultaneously. The stability of the fit is enhanced in the following manner: First we perform an unconstrained fit to extract E stat n −E stat 1 for n = 4, 5, 6 for each lattice spacing and action and compute a rough average of r 0 (E stat n −E stat 1 ) for these values of n. Then we repeat the fit, constraining r 0 (E stat n − E stat 1 ) for n ≥ 4 to the previous average (this renders the fit linear). For n ≤ 3 this is unnecessary as these levels are well determined at each individual lattice spacing. We list the extracted values of r 0 (E stat n − E stat 1 ) in Table 2. Finally, using the values of E stat n and β stat n,N determined from this fit as (fixed) input parameters, we fit E eff,kin n (N, t, t 0 ) and E eff,spin n (N, t, t 0 ) by in the same manner. While the fitted results are quite stable, we consider them as rough estimates only, since our fits include only the leading exponential correction, and there are systematic effects from higher corrections to the GEVP. We therefore employ the fits only to estimate the size of the leading corrections. For a reliable estimate of the energy levels, we calculate plateau averages from t = t min ≥ t 0 on at each N and t 0 , and take as our final estimate that plateau average for which the sum σ tot = σ stat +σ sys of the statistical error σ stat of the plateau average and the maximum systematic error σ sys = (t min ) becomes minimal, subject to the constraint that σ sys < 1 3 σ stat . We impose the latter constraint in order to ensure that the total error is dominated by statistical errors.
An illustration of the more problematic plateaux is shown in Fig. 1. It is rather clear that without some analysis of corrections due to excited states it is very difficult to locate a safe plateau region at least for n = 3.
Our results are given in Table 3; besides the plateaux found by the method described in the preceding paragraph, we also show the results of the global fit, which in almost all cases agrees very well with the final plateau value.

Continuum limit
We now turn to the continuum extrapolation of the level splittings. Using the fact that the static actions employed are discretizations of the same continuum theory, we perform a combined continuum limit by fitting a function of the form (k = 1, 2 for HYP1, HYP2 actions) Φ i,k (a/r 0 ) = Φ i + c i,k (a/r 0 ) s i     at each β and take the continuum limit of the combination. The linear term in a which is present in the combined data, is suppressed by a factor 1/m b . Given in addition the flatness of the data in a (see Fig. 2) we just use s 4 = 2 for the combination, assuming that this term dominates. Figure 2 shows the approach to the continuum limit for the static and full HQET energy splittings. The results for HYP1 and HYP2 are distinguished by color; the static splittings are shown as circles, whereas the full HQET splitting ∆E HQET 2,1 is shown as diamonds. Also shown are the fits to the continuum limit together with their error bands. We see that the approach to the continuum limit is rather flat in particular for ∆E 2,1 , and that the O(1/m b ) corrections constitute only a small shift of the energy splitting between the first excited and ground states. In Fig. 3 we show the approach to the continuum limit for the spin splitting ∆E P−V in the same fashion.

Results
Our findings for the continuum values of the level splittings are summarized in Table 4. (diamonds). Shown are the results for both HYP1 (red, shifted to the left) and HYP2 (blue, shifted to the right).
For the static spectrum, we obtain giving ∆E P−V = −29.8(3.2) MeV via r 0 = 0.5 fm. This is to be contrasted to the experimental value [40] of m Bs − m B * s = −49.0(1.5) MeV. We note that while O(1/m 2 b ) effects may not be entirely negligible for such a small splitting, the difference is too large to be entirely attributed to these. Instead, a genuine quenching effect is involved.
The spin splitting between the first radial excitations in the pseudoscalar and vector channels is r 0 ∆E P −V = −0.056 (27) , (3.12) which is compatible with the spin splitting between the ground states.
Our results confirm the expectation that the O(1/m b ) contributions are small compared to the static results. To get an idea of the size of higher order corrections in 1/m b , we have considered the dependence of the energy level splittings on the matching conditions chosen in the determination [23] of the HQET parameters ω kin and ω spin . We find that the dependence is very minor, with a maximum deviation of δ|r 0 ∆E P−V | = 0.005, less than the statistical errors, and δ|r 0 ∆E 1/m 2,1 | = 0.002, much less than the statistical errors. Comparing this to the naive power counting estimate of |O(1/m 2 b )| ∼ (1/(r 0 m b ) 2 ∼ 1/100, we see that the tested 1/m 2 b terms are as small as expected or even smaller. Note that 1/r 0 ≈ 400 MeV is indeed a typical non-perturbative QCD scale.

Conclusions
In this paper, we have reported results from quenched lattice QCD for the spectroscopy of the low-lying excited states of the B s and B * s systems. An application of the generalized eigenvalue method with all-to-all propagators to non-perturbative HQET at O(1/m b ) allows us to extract precise results for the energies of the lowest-lying radial excitations as well as for the B s − B * s splitting. However, we emphasize again that a careful analysis of systematic errors due to excited state contaminations is necessary.
A first relevant observation to be pointed out concerns the renormalizability of HQET. Unlike for QCD on and off the lattice, there is no proof of renormalizability of the theory to all orders of perturbation theory. However, we find that in our nonperturbative computations the divergences cancel after proper renormalization of HQET [23]. The left over lattice-spacing dependence in Fig. 2, Fig. 3  as seen in Table 3 and weaker but still very prominent divergences are present in E stat n . In other words we find strong numerical evidence for the renormalizability of the theory; in fact also the universality of the continuum limit is demonstrated in the figures. It is also worth emphasizing that the present demonstration is the first time the continuum limit is taken for mass splittings in HQET.
We find the physical O(1/m b ) corrections to be small throughout. The precision attained, in particular when taken together with the relative smallness of the O(1/m b ) effects, indicates that non-perturbative HQET combined with the use of the GEVP for data analysis is a reliable method for determining B meson spectra. We intend to apply it to the N f = 2 case in the near future. In this context one should remark that we were able to achieve good precision using only 100 configurations in our quenched study. Therefore we do expect to be able to decrease the errors for dynamical fermions. However the influence of topological modes being updated only slowly [41] needs to be controlled or better algorithms with a faster decorrelation need to be used. A promising proposal has been made in [42].

A All-to-all propagators
In this appendix, we explain the details of our implementation of all-to-all propagators, which follows the idea of [34] with some useful improvements.

A.1 Even-odd preconditioning
To reduce the computational effort and storage requirement for the matrix inversions, we consider even-odd preconditioning of the (hermitian) Wilson-Dirac operator Q = 2κγ 5 D. With even/odd ordering of the sites one has a block structure where M ee (M oo ) differs from unity by the clover term on the even (odd) sites, and M oe (M eo ) is the hopping term. Defining the preconditioned matrix B † QB is block-diagonal and the propagator can be factorized as

A.2 Approximate low modes and a stochastic estimator
We consider an orthonormal basis {v i : i = 1 . . . N L } of an N L dimensional subspace (N L ≥ 0) of all fermion fields which live only on odd sites. Defining the projectors A natural choice forv i are approximate eigenvectors of the low-lying eigenvalues ofQ oô with v i = 1 andv † ir k = 0. Then, the partQ −1 oo P L in (A.2) is expected to approximate the long-distance behaviour of the propagator [34,43], and the inversions 1 needed in (Q −1 oov i ) are cheap. 1 Since we do not explicitly useQ −1 oovi ≈ λ −1 iv i, the errors ri in (A.3) are allowed to be large. In practice, we require ri ≤ 0.001 · |λi|, and take λ −1 iv i only as start vectors for the inversion.
On the other hand, we can introduce a stochastic estimator for P H . We take random vectors η i with where α, β denote combined (color, Dirac, and site) indices, and . η is the average over η. The relations (A.4)-(A.6) hold, for instance, in the case of a U(1) noise where φ i,α are independently uniformly distributed in [0, 2π) (while for Z 2 noise the average in (A.6) would not vanish for i = j). Thus, the second term in (A.2) can be written asQ −1 and the estimator ofQ −1 oo can be written as a sum of dyadic productŝ The full propagator Q −1 is then obtained from (A.1). Since B connects only adjacent time slices, the blockQ −1 ee does not contribute to the propagator between sites with time separation |x 0 − y 0 | > 2a. In this case, we can simply write Even-odd preconditioning can be seen as a form of dilution [34], since there are only half as many components of the noise field η in the even-odd preconditioned case as without preconditioning. Note, however, that unlike other dilution schemes, even-odd preconditioning does not increase the number of inversions needed.

A.3 Time dilution
In addition, we use the more conventional time dilution scheme. It is implemented by replacing P H in (A.2) by P H t P t where P t projects on the components corresponding to (odd) sites with time coordinate t. Then, an independent stochastic estimator is introduced for each term where the noise vectors η ti have non-vanishing components only for (odd) sites on timeslice t.
Note that due to the hopping term in B the full propagator (A.9) from time slice x 0 to y 0 receives contributions from noise vectors η ti on three time slices, t = x 0 , x 0 ± a, i.e. three inversions are required for the propagator from one time slice x 0 . However, a total of T inversions is sufficient, and hence no extra effort is required, if one computes the propagator for all x 0 , as we do in our measurements.
Analysing the variance of a heavy-light two-point correlator as described in [44], one sees that the variance with time dilution decays roughly as e −(x 0 −y 0 )mπ , while the expression without time dilution contains pieces independent of x 0 − y 0 . This renders time dilution very profitable.