HQET at order 1/m: I. Non-perturbative parameters in the quenched approximation

We determine non-perturbatively the parameters of the lattice HQET Lagrangian and those of the time component of the heavy-light axial-vector current in the quenched approximation. The HQET expansion includes terms of order 1/m. Our results allow to compute, for example, the heavy-light spectrum and B-meson decay constants in the static approximation and to order 1/m in HQET. The determination of the parameters is separated into universal and non-universal parts. The universal results can be used to determine the parameters for various discretizations. The computation reported in this paper uses the plaquette gauge action and the"HYP1/2"action for the b-quark described by HQET. The parameters of the current also depend on the light-quark action, for which we choose non-perturbatively O(a)-improved Wilson fermions.


Introduction
Heavy quark effective theory (HQET) was developed already quite a while ago [1][2][3][4][5]. Still it is of considerable interest today, primarily for two reasons. First it describes the asymptotic expansion of QCD observables in the limit of a large quark mass, in particular the mass of the b-quark, m b . For this to be true, the observables have to be in the proper kinematical region. But in such a region, we expect the expansion in 1/m b to be valid also non-perturbatively if the parameters in the effective theory are determined non-perturbatively. Understanding QCD then includes understanding the HQET limit.
Steps towards establishing the expected equivalence on the non-perturbative level 1 have already been carried out [6,7], but we intend to go much further, in particular through a complete treatment of 1/m b corrections.
Second, flavour physics has become a precision field. No flavour physics observable has shown evidence for physics beyond the Standard Model with the presently available precision. In many cases the limiting part is the theoretical uncertainty, not the experimental one [8]. Reliable lattice computations are needed to make progress, but heavy quarks are difficult due to O((am b ) n ) discretization errors, where a is the lattice spacing. In particular, if am b is too large the expansion in a will break down altogether [9,10] 2 . Instead, in HQET the cutoff effects are O((aΛ QCD ) n ); they are thus much more manageable. Here we present the determination of the parameters of the effective theory in a formulation where all power divergences are subtracted non-perturbatively through a matching of QCD and HQET in a small volume [16] with Schrödinger functional boundary conditions. As illustrated in Fig. 1, we then continue to larger lattice spacings through a step scaling method. Section 3 explains the details, in particular how lattice spacing errors are removed in each step. The work of [17] is here extended by a determination of all parameters in the action as well as those for the weak current A 0 including the terms of order 1/m b .
Our numerical implementation is done in the quenched approximation. It represents a very useful test-laboratory, in particular for the first motivation presented above: the qualitative features of the 1/m b expansion will not depend on the fact that the light quarks are quenched. Applications of the parameters computed in this paper require the determination of energies and matrix elements in a large volume and will appear in separate papers, but already here we learn interesting lessons about the asymptotic convergence of the expansion. In the following section we recall the basic formulation of HQET, with a primary focus on defining the parameters of the theory and how they enter the computation of the spectrum and matrix elements. Section 3 gives a short but precise account of the strategy for the computation of the parameters. The technical Sections 4.1-4.2 discuss the numerical details of the different intermediate steps while Section 4.3 presents our main results for the HQET parameters.

HQET at order 1/m b
In this section we define HQET including terms of order 1/m b , in particular the parameters of the theory. We then show the expansion of a few observables as examples how the parameters can be used. We choose to regularise the theory on the lattice although our approach is in principle independent of the specific regularisation. Almost all formulae are established in [17]. They are repeated here for the benefit of the reader to keep the paper self-contained, but for details, such as the exact consequence of spin symmetry, the reader is referred to [17]. Terms of order 1/m 2 b are dropped without notice.

Lagrangian
The HQET Lagrangian, consists of the lowest order (static) term, and the first order corrections We use the backward covariant derivative D 0 as in [18] and the 4-component heavy quark field subject to the constraints P + ψ h = ψ h , ψ h P + = ψ h with P + = (1 + γ 0 )/2 and the discretized version σ·B = k,j σ kj F kj /(2i) , where σ kj and the lattice (clover) field tensor F are defined in [19]. The kinetic term D 2 is represented by the nearest neighbour covariant 3-d Laplacian. The normalization is such that the classical values of the coefficients are ω kin = ω spin = 1/(2m b ). A bare mass m bare has to be added to the energy levels (e.g. the B-meson mass) computed with this Lagrangian to obtain the QCD ones (up to 1/m 2 b ). At the classical level it is m b , but in the quantized theory, it has to further compensate a power divergence.

Weak axial current
For many applications, the weak, left handed, heavy-light current is needed. We here consider just the time component of the axial-vector part, A 0 . The other components can be treated analogously, but so far we restricted ourselves to A 0 , which is sufficient for a computation of pseudoscalar decay constants.
At the lowest order the current is form-identical to the relativistic one. At first order it is corrected by two composite fields of dimension four. The explicit form is where all derivatives are taken to be the symmetric nearest neighbor ones, By considering Z HQET A to be a function of m l /m b with the light quark masses m l , the above set of operators is complete after using the symmetries and the equations of motion. Since m l ≪ m b and such effects are further reduced by a factor of the coupling constant α(m b ) we ignore this dependence and determine Z HQET A with the light quark mass set to zero. Note that A (2) 0 (x) does not contribute to correlation functions and matrix elements at vanishing space-momentum. It is not needed for the computation of decay constants and has not been written down in [17]. We will also not determine its coefficient c (2) A here. A short note on discretization errors is useful at this point. Remaining in the static approximation (ω kin = ω spin = 0), the Lagrangian has been shown to be automatically O(a) improved by studying its Symanzik expansion. Therefore, all energy levels are as well [19]. The zero space-momentum matrix elements of the current A  [18,19]. 3 Including the 1/m b terms, linear terms in a remain absent, except for those accompanied by a factor 1/m b [16]. So leading discretization errors are O(a/m b ) and O(a 2 ). 3 At zero space momentum, partial summation can be used to bring A (2) 0 into the form used in [19].

Correlation functions
Observables of interest are obtained from Euclidean correlation functions in large volume. As an illustration we consider the QCD correlator of the heavy-light axial current in QCD, A µ = ψ l γ µ γ 5 ψ b (and A † µ = −ψ b γ µ γ 5 ψ l ). It is normalized by Z A to satisfy the chiral Ward identities [20,21]. Ignoring renormalization for a moment, an expansion in 1/m b reads where the notation indicates C In HQET the path integral weight is expanded in 1/m b and one obtains the fully renormalized expansion in terms of bare static expectation values, for example (2.12)

Spectrum and matrix elements
The spectral representation yields the large time behaviour in QCD where ∆ is a gap in this channel ( ≈ 2m π in large volume) and A = B(p = 0)|A 0 (0)|0 in the non-relativistic normalization of states B(p)|B(p ) = 2δ(p − p ) . Of course A = f B m B /2 is a phenomenologically interesting quantity. Again a naive expansion reads (2.14) Using the transfer matrix of the static theory, the HQET correlators (static and beyond) are easily analysed. The result is for the static term and for example for the kinetic correction. 4 A comparison to eq. (2.13) yields where r 0 is an arbitrary length-scale of the theory. Completely analogous formulae hold for the excited state masses and matrix elements. The necessary bare quantities such as E kin can be efficiently computed with the generalized eigenvalue method [22]. We have here chosen to write the expansion of ln(A r 3/2 0 ) instead of A itself, since in this way one explicitly avoids terms quadratic in 1/m b , while e.g. in eq. (2.10) terms such as Z AA are understood to be dropped to remain consistently in order 1/m b . This rule is necessary for a correct renormalization of the theory.

Properties of non-perturbative HQET parameters
Before entering the presentation of a computation of the HQET parameters, we here mention some of their properties.

Scheme dependence
Since we are working in the framework of an effective field theory, which has non-trivial renormalization, the parameters have to be determined by a matching performed at finite value of the expansion parameter 1/m b . The truncation of the effective theory then introduces an ambiguity of order of the left out terms.
As a result the parameters such as ω kin and Z HQET A have a dependence on the choice of the matching condition. This is a scheme dependence, analogous to the one of the usual perturbative expansion in the renormalized QCD coupling. For an explicit example consider ln(Z HQET A ). Evaluated in the static approximation we will denote it by ln(Z stat A ) and including all 1/m b terms it is denoted as ln(Z HQET ). The scheme dependence of ln(Z HQET A ) is then of order 1/m 2 b , while the one of ln(Z stat A ) and ln(Z ) individually is of order 1/m b . Note that in all this discussion we are working non-perturbatively in the QCD coupling but order by order in 1/m b . 4 An explicit expression for E kin is where the state |B is an eigenstate of the static transfer matrix.

Difference to previous renormalizations in the static approximation
In order to avoid a misunderstanding we also point out the difference to the renormalization carried out for example in [23,24]. As long as one is working just in the static approximation, the weak currents do not mix with operators of different dimensions. It is then also consistent to renormalize them perturbatively. Both to avoid ambiguities in renormalization schemes and to be able to profit from the high order continuum perturbative results [25,26], it is advantageous to first introduce the RGI current, The renormalization constant Z stat A,RGI has been determined non-perturbatively in [23] for the quenched case, and in [27] for two flavours of dynamical quarks.
However, here we determine Z stat A by a non-perturbative matching to QCD, and it is different. In fact the correspondence is in terms of the matching function C PS introduced in [24]. In the "old" strategy of for example [24,28], functions such as C PS are determined from high-order perturbation theory, while here we evaluate the full factor Z stat A non-perturbatively. We collect all HQET parameters into one vector ω with components ω i , i = 1, . . . , 5 listed explicitly in Table 1. In static approximation the parameters are ω 1 , ω 2 . When 1/m b corrections are included, the additional parameters in the action are ω 4 , ω 5 ; moreover the previous parameters change by (partially power divergent) terms of order 1/m b . The situation with ω 3 is more intricate. It is needed for O(a)-improvement of the static approximation and for genuine 1/m b -terms at order 1/m b . We now turn to an explanation of the various steps involved in the determination of the ω i . Our strategy is illustrated in Fig. 2.

Computation of HQET parameters
HQET QCD match a ω ω Figure 2: Illustration of the strategy. Our numerical application uses much finer resolutions a/L than those shown here. For each step i a series of simulations Si is necessary, they are described in the text and in Table 2.

Finite volume observables
Five observables Φ i are required to determine the HQET parameters ω i . We choose Φ i (L, M, a) defined from Schrödinger functional correlation functions by forming suitable renormalized combinations. They are universal, which means in particular that their continuum limit Φ i (L, M, 0) exists.
As variables we have chosen the box size L (which plays the rôle of a kinematical variable), the RGI mass of the heavy quark, M , and the lattice spacing a. Equivalent but dimensionless variables are a non-perturbative running couplingḡ(L) (e.g. in the Schrödinger functional scheme), the combination z = M L and the resolution a/L. The light quark masses are assumed fixed (in the numerics they are set to zero).
Most of the details of our choice follows closely [17]; for others we refer the reader to App. A. Here we note only a few properties. As L becomes large, Φ i , i = 1, 2 tend to the B-meson mass and the logarithm of the B-meson decay constant, up to kinematical constants. They are thus mainly determining ω i≤2 . Correspondingly Φ 4 and Φ 5 allow easy access to the 1/m b -parameters in the Lagrangian and Φ 3 to the correction term of the current with coefficient c (1) A . By a choice of notation (e.g. considering ln(A) instead of A), we arranged our finite volume observables to be linear in ω, (3.1) By construction φ is a block diagonal matrix, which is explicitly given in terms of (bare) correlation functions in HQET. The inhomogeneous pieces η i involve just correlation functions in the static approximation. Their continuum limits η i (L, 0) exist for i > 2, while for i = 1, 2 additive renormalizations are necessary. They are contained in the second term of eq. (3.1). In the case i > 2, the two terms of the right hand side of eq. (3.1) are then computed separately, see section 3.5.
The HQET parameters ω are defined by matching at a certain value of L = L 1 the HQET observables Φ i to the corresponding values Φ QCD i computed in QCD. These QCD observables are first extrapolated to the continuum limit (indicated by S 1 in Fig. 2 (3.2)

Matching
The matching scale L 1 must be chosen such that 1/L 1 m b to allow for a precise expansion in 1/m b . Moreover, we want lattice spacings of order 10 −2 fm in order to keep am b < 1/2 while performing the continuum extrapolation eq. (3.2). These constraints lead to a box size L 1 ≈ 0.4 fm [17,24].
The HQET parameters are then defined by imposing for any value of the lattice spacing. By inverting eq. (3.1) one obtains Through our choice the matrix φ is of the form where C = diag(L, 1). The 3 × 3 sub-matrix A has a further upper triangular structure. It is written explicitly in App. B, together with its inverse.

Step scaling
For reasonable resolutions L 1 /a ≥ O(10), only lattice spacings a 0.05 fm are accessible, while for standard large volume HQET computations we would also like larger a.
In order to obtain ω also there, we first need Φ(L, M, 0) at larger L = L 2 . We simply use the valuesω(M, a) of eq. (3.4) and determine the continuum limit of the HQET observables at L 2 (indicated by S 3 in Fig. 2) This can be done as long as the lattice spacing is common to the n 2 = L 2 /a and n 1 = L 1 /a-lattices and is kept at a fixed, small, ratio. 5 We choose s = 2 in the numerical application. Typical resolutions are L 1 /a = O(10), see Sect. 4. This procedure which takes us from L to sL with a finite scale factor s is called step scaling [29]. The explicit forms of the step scaling functions are given in App. B.

HQET parameters
The parameters ω i for use in large volume (see Sect. 2.3.2) are finally obtained from

Splitting lowest order and first order in 1/m b
In Sect. 2.4.1 we discussed that the splitting of a prediction into different orders of the 1/m b -expansion is not unique. Nevertheless, it is of interest to organize the calculation into a static one and the remainder. In this way one can judge on the generic size of 1/m b corrections and thus get some indication of the asymptotic convergence of the series. A second reason is that the static theory is O(a) improved when ω 3 is included and then represents the improvement coefficient ac stat A = ω 3 [19]. Thus the by far dominating part of the result can be extrapolated quadratically in a and only a small correction has to be extrapolated linearly to the continuum limit. Since a nonperturbative determination of c stat A has not been carried out, we will here use its one-loop perturbative value [30]. We anyway carried out quadratic extrapolations in a and then studied the effect of incomplete improvement. We will see in Sect. 4 that c stat A is of little relevance. Therefore we profit from the O(a)-improvement of the static theory, irrespective of the precise value of c stat A .
To summarize, the static approximation is defined by using exactly the same formulae as in HQET with 1/m b but setting All quantities are determined once in this static approximation and once with 1/m b terms included and then the pure 1/m b -corrections are given by For any new observables, in particular in large volume, the static results are given by inserting ω stat into the expressions such as eq. (2.19) and the 1/m b -correction by inserting ω (1/m) instead. Since, as we discussed earlier everything is to be linearized in the ω, the "full" result up to 1/m 2 b is obtained from summing static and 1/m b -correction or by using directly ω. We finally note that even though the "full" result contains O(a) discretization errors, it appears justified to extrapolate numerical data with a leading correction term ∝ a 2 since the linear terms are suppressed by a small 1/m b factor, which one would estimate e.g. as 1/(m b r 0 ) ≈ 1/10.

Numerical results
The numerical computations have been performed in the quenched approximation. For the QCD part we used non-perturbatively O(a)-improved Wilson fermions [31], see [17] for details of the Schrödinger functional implementation. Two discretizations of the static quarks have been considered: the so-called HYP1 and HYP2 actions [18,32]. In order to reduce discretization effects, we have implemented tree level improvement of our QCD observables and of the HQET step scaling functions (see App. D).
Our observables depend on three periodicity angles θ 0 , θ 1 , θ 2 , (see App. A). As in [17] we considered all combinations of those with θ i ∈ {0, 0.5, 1.0} , θ 1 < θ 2 , i.e. nine different matching conditions. For our discussion of the results and the extraction of the final HQET parameters we chose a "standard set" because this yields overall the smallest statistical errors in the HQET parameters. We will comment on the spread of results with other choices for θ i as we go along. We summarise the different simulations needed in our strategy in Table 2. Since a size L = L 1 is used for the matching (eq. (3.2), eq. (3.4)) both QCD (S 1 ) and HQET (S 2 ) are simulated in that volume. The volume of space extent L 2 was simulated (in HQET), with two different sets of lattice spacings: • First (S 3 ) with the same set of lattice spacing as used in S 2 , to evaluate eq. (3.6).  • Then (S 4 ) with lattice spacing of order 0.05 . . . 0.1 fm, in order to compute ω(M, a) given by eq. (3.8).
The choice of the simulation parameters is described in detail in [17]; Table 3 of that paper lists those for S 1 , Table A.1 of [16] the parameters of S 2 and S 3 and Table 6 of [18] the parameters of S 4 . Here we just note that L 2 = 2L 1 and L 1 is fixed by the Schrödinger functional couplingḡ 2 (L 1 ) = u 1 = 3.48. Since this condition was implemented only within a certain precision, a small mismatch ofḡ 2 (L) =ũ 1 used in S 1 andḡ 2 (L) = u 1 used everywhere else has to be taken into account. This is done in complete analogy to appendix D of [17]. Without discussing the details we will quote the small corrections below. We recall that Φ 1 is simply a finite volume pseudo-scalar meson mass (up to a normalization) and determines the quark mass. In Fig. 3, we show its continuum extrapolation for our chosen values of the heavy quark mass, corresponding to z = L 1 M = 10.4, 12.1, 13.3. We have taken into account the errors coming from the relation between the bare quark mass and the RGI quark mass. In particular, a part of this error is common to the data at all lattice spacings. It is included after the continuum extrapolation; the final error bar is shown on the left side of this plot (see also [17]). Since O(a) improvement is implemented, this as well as all other extrapolations of eq. Out of the O(z 0 ) observables, Φ 2 determines the normalization of the axial current. Its continuum extrapolation, eq. (3.2), is also illustrated in Fig. 3.
Some of the static quantities η i (L 1 , a) have a well defined continuum limit (see appendix B). In these cases, η i (L 1 , a) is replaced by η i (L 1 , 0) in eq. (3.4). We show examples of the determination of the continuum limit η 3 (L 1 , 0) and η 4 (L 1 , 0) in Fig. 4. The graphs use the one-loop values of the improvement coefficient [30] c stat A = 0.0029 g 2 0 for the HYP1 action .  (L1, a). The three different choices (θ1, θ2) = (0, 0.5), (0.5, 1), (0, 1.0) are represented by squares, circles, and diamonds respectively. In each case two different static actions, HYP1 (the red points, slightly shifted to the right for visibility) and HYP2 (the blue points) are used. The continuum limit is obtained by a constrained fit ηi(L1, a) = ηi(L1, 0) + ci,j a 2 /L 2 1 , with j = 1, 2 for the two different actions.
We turn now to Φ 4 and Φ 5 , needed for the determination of the 1/m b parameters ω kin and ω spin . Their continuum extrapolations are shown in Fig. 5. Due to the exact spin symmetry of the static effective theory, Φ 5 has no static contribution. Note that both Φ 5 and the pure 1/m b part of Φ 4 , obtained after subtraction of η 4 (L 1 , 0) (see Fig. 4), are an order of magnitude smaller than Φ 4 . As expected Φ 4 − η 4 and Φ 5 are decreasing functions of the quark mass.

Observables for
The step to observables in the larger volume is described by eq. (3.6). It can be broken up into several step scaling functions defined in appendix B, which individually have a continuum limit. Discussing them one by one would be too lengthy and is also not very illuminating. We follow exactly Sect. 3 and insert eq. (3.4) into eq. (3.6) lattice spacing by lattice spacing and then extrapolate a/L 2 → 0.
Two examples are shown in Fig. 6. We observe that having the data for two different static actions is very useful to constrain the continuum limit, particularly so for the 1/m b parts. Also the resolution L 2 /a = 32 which is in addition to those of [17] helps a lot (the reader is invited to compare Φ (1/m) 1 in Fig. 6 to Fig. 5 of [17]). In all cases the 1/m b corrections are much smaller than the leading terms, suggesting a good asymptotic convergence of the 1/m b expansion. The precision of the 1/m b term is worse than the static one for Φ 2 , partly since the latter can be extrapolated quadratically in a to the continuum. For Φ 1 this is different because the overall error contains a large piece coming from the renormalization factor determining the RGI quark mass in QCD.
For reasons of numerical precision, Φ 5 (and only Φ 5 ) is not computed exactly as described in Sect. 3. Its definition (App. A) involves the propagation of a heavy quark over a distance T = L, introducing significant statistical errors for large L/a in the effective theory. These become unpleasantly large in eq. (3.6), more precisely in φ 55 (L 2 , a). We therefore replace Φ 5 (L 2 ) byΦ 5 (L 2 ) differing only by the choice T = L/2. An obvious question is why this is not done already for L = L 1 . The reason is thatΦ 5 (L 1 ) turns out to be quite a bit smaller than its natural order of magnitude of O(1/z). In such a situation O(1/z 2 ) terms may be numerically comparable and the 1/m b expansion may be compromised in the matching step. We therefore chose the described solution, even if it is lacking elegance. 6 The continuum extrapolation ofΦ 5 (L 2 ) is shown in Fig. 7.
Since the computed Φ(L 2 , M, 0) may be used in the subsequent step eq. (3.8) also with lattice discretizations which differ from ours we list Φ i (L 2 , M, 0) in Table 3 and Table 4. Starting from these numbers, the remaining computations to obtain the HQET parameters for a different lattice action require a very modest effort.

Renormalization group invariant b-quark mass
For each value of z = M L 1 , eq. (3.8) yields the desired HQET parameters. However, for future use we want to list them for M = M b . This saves space here and somebody doing a computation in the future does not have to fix the right quark mass. The b-quark mass was already obtained in [17] from the experimental spin-averaged B s meson mass and r 0 = 0.5 fm. Here we repeat its determination using the mass of the pseudoscalar B s meson as experimental input. This is the natural quantity since we also start from the effective mass defined from the f A correlator to fix the b-quark mass in the small volume. Moreover, we have improved some of the necessary steps through the finer resolution a/L 2 = 1/32, and the use of tree level improvement (see App. D). With our new, improved determinations [33] using the GEVP method [22]   These results are in perfect agreement with the ones obtained when using the large volume numbers E kin , E stat of [17]. In the following we use eq. (4.3) and the knowledge of L 1 /r 0 and interpolate all results quadratically in z to z b = 12.30 (19) in the static approximation and to z b = 12.48 (20) at first order in 1/m b .

Bare parameters ω i
Our ω i determined from eq. (3.8) and simulations S 4 are listed in Table 5 for our standard θ−combinations. The errors take all sources into account (through a jackknife analysis incorporating all steps), but one has to be aware that there are very significant correlations in the parameters. We discuss these correlations in App. C and provide them in tables available on the Internet. In static approximation a small shift proportional toũ 1 − u 1 (see the discussion at the beginning of this section) is applied. We refer to Appendix D of [17] for the case of the b-quark mass, and with similar considerations we have evaluated the effect on the current renormalization, ln(Z stat A ). We found 8  Generically the bare parameters ω i are completely non-universal and depend on all details of the action. However here we are working in the quenched approximation. In this situation the HQET action can in principle be determined independently of the light quarks. Thus ω 1 , ω 4 and ω 5 can be used for any light quark action also different from our specific one.
We illustrate the cutoff-dependence of ω 1 = m bare as a function of L 1 /a = L 1 Λ cutoff . In Fig. 8 we show the static bare quark mass and its 1/m b contribution in units of L 1 , using L 1 /a = L 2 /(2a) and the data of Table 5. In addition, smaller lattice spacings are covered by including the numbers forω 1 from eq. (3.4). The two sets of bare parameters ω i and ω i differ by cutoff effects: the former are determined directly for L = L 1 and the latter after a step scaling to L = L 2 . Indeed these cutoff effects are visible but 7 The reader might wonder why for r0M b we don't obtain a better precision compared to the result quoted in [17]. The reason is that the uncertainty which affects the b-quark mass is dominated by the renormalization of the quark mass in QCD [17]. 8 These shifts have to be added to the raw numbers computed with our simulations S1, . . . , S4 but are already included in the numbers given in eq. (4.3) and in Table 5.   not dominating. The largest part of the variation with a is due to the divergences. In the static case, the divergence is known to one-loop order from [18] (see table 1 of that paper). It is plotted for HYP1 and HYP2 using both standard (dotted lines) and boosted (dashed lines) perturbation theory, fixing the constant piece at the smallest lattice spacing. The non-perturbative results agree qualitatively with the perturbative approximations. We show the equivalent plots for ln(Z stat A ) and ln(Z ) the theoretically expected 1/a divergence is not clearly visible. It either has a rather small coefficient or it is masked by terms with positive powers of the lattice spacing. For illustration we nevertheless fit a linear behaviour to ln(Z (1/m) A ) at the smallest three lattice spacings and extend this curve to larger ones. For ln(Z stat A ) we show the leading logarithmic divergence, ln(Z stat A ) ∼ ln(a/L 1 ) g 2 0 /(4π 2 ) + constant, in the graph. Again the constant is adjusted to the data point at smallest a and replacing g 2 0 with a boosted coupling defines the boosted perturbation theory expression.
Before closing this section we would like to add a remark concerning different matching conditions. By construction, the observables Φ depend on the values of the angles θ i . Nevertheless, since the HQET parameters ω i are bare parameters of the HQET Lagrangian and fields, they are θ independent up to truncation corrections of order (1/m b ) n . As one can see from eq. (3.1), this implies that the θ-dependence of Φ is absorbed by the ones of η and φ. In practice this means that a parameter ω computed in the static approximation, such as m stat bare or Z stat A exhibits a small θ-dependence, but once the 1/m b corrections are added, this dependence has to be further suppressed to order (1/m b ) 2 . We checked that we obtain this behaviour in our numerical simulations. In the static case, we also show a comparison with standard and boosted perturbation theory [18], as described in the text. Data for ω1 from Table 5 (large volume simulations S4) are represented by circles (the three points on the left), whileω1 obtained from the small volume simulations S2 are represented by squares. We use the colour blue for HYP2 and red for HYP1 (also slightly shifted to the right). The results are shown for the central mass and for the standard set of θi. It is illustrated for the case of ω 4 in Table 6.

Outlook
Our non-perturbatively computed HQET parameters show a qualitative agreement with perturbation theory as far as the leading divergence in each of the parameters is concerned (see figures 8,9), but this does not extend to the quantitative level needed for precision flavour physics.
Fortunately, in the quenched approximation, we now have the full set of parameters for HQET spectrum calculations as well as for (zero space-momentum) matrix element of A 0 including the terms of order 1/m b -all of them are known non-perturbatively. A detailed test of the 1/m b expansion is thus possible. In companion papers we are carrying this out for the examples of the B s decay constant and for some mass splittings in the B s system; some preliminary results are described in [33]. The parameters do of course depend on our choice of discretization, but starting from Tables 3 and 4, the effort to compute them for a different regularization is very modest. Note that, as long as one remains in the quenched approximation, the parameters for the HQET action do not need to be recomputed if one uses just a light quark action differing from ours.
The small volume simulations for a determination of the parameters for N f = 2 are already far advanced, see [7] for a recent account. Therefore, the situation will soon be similar with 2 flavours of dynamical fermions. R A = ln(f A (T /2, θ 1 )/f A (T /2, θ 2 )| T =L is easily seen to have a sensitivity to the coefficient ω 3 which is approximately proportional to θ 2 − θ 1 : one just has to note that the covariant derivatives in eq. (2.5) acting on the quark fields are proportional to their momentum. This free theory argument is valid qualitatively in small volume. In the same way one sees that the combination ) has a sensitivity proportional to θ 2 2 − θ 2 1 to ω 4 , while ω 5 does not contribute due to spin symmetry.
These qualitative considerations combined with some numerical experiments lead us to introduce

B Explicit form of step scaling functions
From our definitions we find immediately  By choosing only L 1 and a as arguments we have assumed that L 2 /L 1 = s is fixed (typically to s = 2) which means D = diag(s, 1). The continuum limit of each element of the step scaling functions exists and the above split into two pieces is suggested by the fact that the limit a → 0 of η i (L, M, a) may be performed for i ≥ 3, while for i ≤ 2 there is an additive renormalization which only cancels in eq. (B.6). Splitting η accordingly, η = η a + η b with η a i≥3 = 0 , η b i≤2 = 0, we can also unify the step from L 1 to L 2 to (note Σ(L, a) i≥3 = 0) Let us now list explicit expressions for the various matrices. The expansion of the observables in HQET follows directly from the expansions detailed in section 2.3 of [17]. The inhomogeneous part is given by the static quantities Note that f stat A does not contain the improvement term proportional to c stat A . The improvement term is included as explained in section 3. From the definition eq. (B.6), one has (B.14) The matrices A, B which make up φ as in eq. (3.5) are In order to avoid unnecessary repetitions, we have used the shortcut y ∈ {kin, spin} .
(B. 25) in the previous formulae. From our definitions we obtain and In static approximation the only non-vanishing elements are Σ stat 13 = Σ 13 , Σ stat 23 = Σ 23 .

C Covariance matrices
When our parameters are used in a subsequent computation, one has to note that their errors are correlated. We therefore supply covariance sub-matrices as far as they appear necessary in practice. They are still too large to be printed (and typed in). We therefore supply them for download together with the paper under http://arxiv.org/ and at http://www-zeuthen.desy.de/alpha/public tables/. All parameters are given in lattice units a = 1, and for our standard choice of θ angles. As explained in the text, we give the HQET parameters and their covariance matrices for c stat A fixed to its one-loop value [30] and fixed to 0.
We normalize the covariance matrix between two quantities O α and O β in the following way where represents the average over N jack = 100 jackknife samples. The covariance matrix is such that C βα = C αβ and C αα = 1 (no summation).
The second subset is the one relevant for a computation of the pseudoscalar decay constant. There the indices α and β take the value α, β = l + 5(i − 1) + 15(k − 1) . . while i, k have the same meaning as before. Thus we obtain a 30 × 30 matrix C b αβ , and give its 435 relevant components. The heavy quark mass is fixed by matching to the spin averaged B s -meson mass. This matching is done at the static order for ω stat and including the 1/m b terms for ω (1/m) .

D Tree level improvement
Perturbative improvement of the observables [36] has been proven to be effective in reducing cutoff effects for example in the case of the step scaling function for the coupling [37,38], where it has been pushed to the two-loop order. The idea is to remove all the O(( a L ) n ln( a L ) m ) at a given order in theḡ 2 expansion of the renormalized lattice observable. At tree level, as done here, this produces classically perfect observables in the sense of [39], by removing all the O(( a L ) n ) effects. The tree-level improved observable is either defined as Here O is the non-perturbative observable evaluated by Monte Carlo, O impr is the treelevel improved one and O tree is the same observable evaluated for g 0 = 0. Generically we use eq. (D.1), but when O tree (0) vanishes eq. (D.2) is appropriate. We also choose the latter for the step scaling functions Σ since it is more natural given their form eq. (B.14), eq. (B.15). We apply this improvement to all QCD and HQET observables used for the matching in small volume and to all HQET step scaling functions.
At tree level there are some useful relations among the Schrödinger Functional correlation functions, which can be used in the perturbative computation [40]: Since in our QCD simulations L is fixed throughḡ 2 (L), each value of z determines a value of M . In the corresponding tree level computation the improved subtracted quark mass 9 m q = m q (1 + b m am q ) enters. At the considered order we could set m q to m(µ) at an arbitrary scale µ in an arbitrary scheme. We decided to choose m q = m , where m ≡ m univ (m ) and m univ (µ) is the running mass computed from the RGI mass (but) using only the universal part of the β (i.e. the coefficients b 0 and b 1 ) and τ (i.e. the coefficient d 0 ) functions. This choice is scheme independent and it is based on the expectation that cutoff effects are dominated by scales around µ = m . The explicit relation between z and m , is then solved for m . In our non-perturbative computation, the Φ QCD i have been computed for z = M L = 10.4, 12.1 and 13.3 withḡ 2 (L) =ũ 1 (see the beginning of Sect. 4) and correspondingly Λ SF L = 0.195 (16). Our tables, provided under http://www-zeuthen.desy.de/alpha/public tables/, are for these specific values.
As an example, the non-perturbative Φ QCD 4 with (right) and without (left) tree level improvement of the observable are compared for (θ 1 , θ 2 ) = (0.5, 1) in Fig. 10. Turning now to HQET, the correlation functions are defined in [17] and useful relations at tree level are The last relation holds when the magnetic background field vanishes, which is the case in our application. It is instructive to consider now the step scaling function Σ 44 , which has also been discussed in [17] 10 , but without tree level improvement. Explicitly with Σ tree 44 (0) = 0.5 . The first factor is responsible for cutoff effects linear in a, which appear because we have not taken care to O(a) improve the theory at order 1/m b . The correction δ is then large and tree-level improvement leads to a strong reduction of the cutoff effects in the non-perturbative Σ 44 , as illustrated in Fig. 11.
For future applications, this discussion suggests to eliminate linear (tree level) a/L effects such as those in Fig. 11 from the start. This is achieved by implementing O(a)improvement at tree-level in the order 1/m b terms of the action and the effective fields. It just means to add terms to action and currents. The normalization is chosen such that a coefficient ω kin will implement tree-level improvement, A A A A 0 + aω kin A 0 ] , (D.14) 10 There it is called Σ kin . and similarly for the action. The effect of these terms is that the kinetic operator D 2 is inserted only with weight 1/2 on the initial and final time-slice of the correlation functions such as eq. (2.12) or f kin A , f kin 1 , k kin 1 . After the Wick-contraction (in a given gauge background field) this corresponds to a standard discretized representation of the integral over the time position of the insertion of D 2 in the static propagator. The spin operator O spin can be treated in complete analogy, but due to eq. (D.8) this has no effect at tree-level. We provide the tree-level improvement coefficients under http://www-zeuthen.desy.de/alpha/public tables/.