Parameters of Heavy Quark Effective Theory from Nf=2 lattice QCD

We report on a non-perturbative determination of the parameters of the lattice Heavy Quark Effective Theory (HQET) Lagrangian and of the time component of the heavy-light axial-vector current with Nf=2 flavors of massless dynamical quarks. The effective theory is considered at the 1/mh order, and the heavy mass mh covers a range from slightly above the charm to beyond the beauty region. These HQET parameters are needed to compute, for example, the b-quark mass, the heavy-light spectrum and decay constants in the static approximation and to order 1/mh in HQET. The determination of the parameters is done non-perturbatively. The computation reported in this paper uses the plaquette gauge action and two different static actions for the heavy quark described by HQET. For the light-quark action we choose non-perturbatively O(a)-improved Wilson fermions.


Introduction
Particle physics enters very exciting times with the start of collecting data at the Large Hadron Collider. Two ways are explored to probe New Physics (NP): either performing a direct search of new particles (e.g. at ATLAS and CMS) from the electroweak scale up to the TeV scale, or studying rare decays (e.g. at LHCb). The latter give rise to a very rich set of constraints on NP scenarios because they are either mediated by quantum loops, in which high energy particles circulate (that is the case, for instance, in flavor changing neutral currents), or by new current structures (in decays occuring at tree level). B-mesons offer a highly interesting and rich set of (rare) decay channels, such as b → sγ, which arises in the standard model only from penguin diagrams, and hence puts strong bounds on NP scenarios. The analysis of inclusive decays such as B → X s γ, in particular in the framework of a heavy quark expansion, strongly depends on the knowledge of the b-quark mass, m b . In tests of the CKM mechanism, some tensions exist at the moment between sin 2β, obtained from the golden mode B → J/ΨK s , and the CKM matrix element V ub extracted from the B → τ ν leptonic decay [1]. The theoretical input of the latter is the B-meson decay constant f B , whose uncertainty is ∼ 10%. Currently, there is also a 3-σ discrepancy between V B→τ ν ub and V B→πlν ub . Although it would be surprising if this difference were due to a significant underestimate of f B , reducing its error is an important task of lattice QCD.
Lattice QCD will enable us to compute m b and f B with a precision comparable to the one of forthcoming experimental measurements from high luminosity collisions. Nevertheless, a delicate issue for those extractions is how to get a satisfying control on the cut-off effects. Indeed, the Compton length of the b-quark is smaller than the typical finest lattice spacing of simulations in large volumes (used to compute hadronic quantities). Several strategies have been explored in the literature to circumvent this two-scale difficulty [2][3][4][5][6][7][8], see [9] for a recent review. The ALPHA collaboration has proposed to use the framework of Heavy Quark Effective Theory (HQET) [10,11] with a non-perturbative determination of the couplings [12]. Its implementation at the first order in 1/m h has successfully been applied in the quenched approximation [13][14][15][16]. In this paper we will report on our effort to realize our HQET program for N f = 2 flavors of dynamical quarks, leaving the phenomenological results for forthcoming papers. Here, in particular, we present our determination of the couplings of the effective theory regularized on the lattice.
We use the same notations as in [14]. The HQET Lagrangian density at the leading (static) order is given by (1.1) A bare quark mass m stat bare has to be added to the energy levels E stat computed with this Lagrangian to obtain the physical ones. For example, the mass of the B-meson in the static approximation is given by At the classical level m stat bare is simply the (static approximation of the) b-quark mass, m b , but in the quantized lattice formulation it has to further compensate a divergence, an inverse power of the lattice spacing. Including the 1/m h terms, the HQET Lagrangian reads 1 At this order, two other unknown parameters appear in the Lagrangian, ω kin and ω spin . Our normalization is such that the classical values of the coefficients are ω kin = ω spin = 1/(2m h ). In order to compute the decay constant of a heavy-light meson, one needs the time component of the axial-vector heavy-light current A 0 . At the lowest order of the effective theory, the current is form-identical to the relativistic one. At the 1/m h order, it is enough to add only one term to the static current (because we are only interested in zero-momentum correlation functions, see [17]) A , ω kin , and ω spin at values of the lattice spacing relevant for the computation of hadronic observables. These parameters allow us to compute, for example, the spectrum of heavy-light mesons and heavy-light decay constants. As we explain in detail in the remainder of this paper, the basic idea is to match, in a small volume, a few observables expanded in the effective theory at finite lattice spacing to their non-perturbative continuum values determined in QCD at a given renormalization group invariant (RGI) heavy quark mass M . Collecting the observables in a vector Φ, the matching equation reads where L is the space extent of the lattice and a → 0 taken in QCD. Expanding the left hand side of eq. (1.6) at a given order of the inverse heavy quark mass defines a set of HQET parameters. We remind the reader that such an order-by-order treatment [12] is part of the very definition of HQET. The remainder of the text is organized as follows: in Sect. 2 we summarize the strategy of the computation, in Sect. 3 we give the details of its implementation, the final results can be found in Sect. 4 and Sect. 5 contains our conclusions. Some definitions are relegated to Appendix A, in Appendix B we explain how we tuned the parameters of the simulations and how we performed the renormalization of the QCD quantities, whereas Appendix C contains more technical information about the dynamical fermion runs.

Strategy
The computation reported here is done along the lines of [14]. For completeness, we repeat here the basic ingredients but refer the reader to this work for more detailed explanations. We start by the simulations of QCD in a small volume with space extent L 1 ≈ 0.4 fm at four different values of the lattice spacing. We consider N f = 2 dynamical light quarks that we tune to be massless and a quenched heavy (valence) quark that we simulate at nine different mass values, such that the lightest mass is around the charm mass and the heaviest mass is above the b-quark mass. We compute five renormalized observables that we extrapolate to the continuum The definitions of these observables can be found in [14] and in Appendix A. Here we just mention that Φ 1 and Φ 2 are finite volume versions of the heavy-light meson mass and the logarithm of the decay constant, respectively (up to kinematic factors). Φ 3 is sensitive to the correction of the heavy-light current. Finally, Φ 4 and Φ 5 are proportional to the kinetic and magnetic corrections, respectively. At the next-to-leading order of the effective theory (i.e. keeping the static and the 1/m h terms) these observables can be expressed in terms of the five parameters discussed in Sect. 1, which we cast into a vector A , ω kin , ω spin t . (2.2) More precisely, we write the 1/m h expansion of the observables in the following way 2 : 3) The entries of the five-component vector η and the five-by-five block-diagonal matrix ϕ are computed by performing a series of numerical simulations of HQET at fixed L = L 1 and for various lattice spacings a. A more explicit form of eq. (2.3) can be found in Appendix A.
As anticipated in the introduction, the matching condition that we impose is eq. They are the bare couplings of the theory and, as such, can be determined by finite volume matching conditions. By imposing eq. (2.4), the parametersω i become functions of M , but this heavy quark mass dependence comes entirely from Φ QCD . We then perform another set of simulations of the effective theory in a larger volume of space extent L 2 = 2L 1 . The observables in this volume are then simply obtained by taking the continuum limit of eq. (2.3) in which we insert the parametersω(M, a) computed in the previous step: Our formulation of the theory and the non-perturbative determination ofω guarantees that all divergences, including those of order 1/a , 1/a 2 are cancelled, and that the limit a → 0 exists. Next, the parameters at larger lattice spacings (to be used in large volume) are obtained by inverting eq. Finally, in the last step we perform an interpolation (or, in one case, a slight extrapolation) in the inverse bare coupling β = 6/g 2 0 and obtain ω(M, a) at exactly those values of the lattice spacing used in our large volume simulations [18].

Continuum extrapolation of the QCD observables
For the QCD simulation, we use the plaquette gauge action and non-perturbatively O(a)improved clover fermions [19] for N f = 2 flavors of dynamical quarks with Schrödinger functional (SF) boundary conditions. We have four different values of the lattice spacing ( 0.02 fm), β varying in the range 6.16 − 6.64, and the number of lattice points per space direction being L 1 /a = 20, 24, 32, 40. The physical volume is kept fixed by imposing for the Schrödinger functional coupling the valueḡ 2 (L 1 /2) = 2.989, which corresponds to L 1 ≈ 0.4 fm. For each resolution L 1 /a we tune the dimensionless RGI heavy quark mass such that With this choice, M varies approximately from 2 GeV to 10 GeV. More details about their renormalization and about the tuning of the bare coupling and quark masses can be found in Appendix B. For the run parameters we refer the reader to Appendix C. Finally, we also implement tree-level improvement, following exactly the procedure described in Appendix D of [14]. Our strategy for the continuum extrapolation differs somewhat from our previous work. The discretization effects can be important for our heaviest masses. In particular, for the simulations where L 1 /a ≤ 24, we might have noticeable contributions of order (a/L) n with n > 2. We take advantage of the fact that we have various (and quite different) heavy quark masses: since the cut-off effects are smooth functions of a/L 1 and z, we perform a global fit of the form using only the data points such that aM < 0.7, as motivated in [20,21]. Note that the two last terms in eq. (3.2) are proportional to (a/L 1 ) × (aM ) and (aM ) 2 . For each observable, Φ QCD i , the parameters of the fit are the nine different values Φ QCD i (L, M, 0) and A i , B i , C i . As an illustration we show the results of the fit of Φ 1 and Φ 2 in Fig. 1. We have checked that different fit ansätze (e.g. adding cubic lattice artefacts) give consistent results for both the central values and the errors, and that results are also compatible with the standard approach where the slope in (a/L) 2 is not constrained.

Subtraction of the static part
The effective theory is first simulated in the small volume of space extent L 1 , which is tuned to be the same as in QCD 3 . We use five different lattice spacings, such that L 1 /a = 6, 8, 10, 12, 16 (but for the continuum extrapolation we discard the coarsest point). The corresponding βvalues lie in the range 5.26 − 5.96. For the static quark we use two different lattice actions, HYP1 and HYP2, which are known to lead to small statistical errors [22]. Once again the light quarks are tuned to be massless, both in the valence and in the sea sectors. More details on the simulations can be found in the appendices. This set of simulations serves to compute the quantities η(L 1 , a) and ϕ(L 1 , a) in eq. (2.5). For i = 3, 4, 5, the static contributions of Φ i are η i . They have a well defined continuum limit and η 5 = 0. Thus we compute and show the continuum extrapolations in Fig. 2. They are done linearly in (a/L 1 ) 2 , which is justified by two reasons: (i) The quantity η 3 contains the time-component of the static axial current; we improve this current using the 1-loop value of ac stat A computed in [23]. This is sufficient since the improvement term has almost no numerical influence: even setting c stat A = 0 gives compatible results. (ii) The quantity η 4 is constructed from SF boundary-toboundary correlators, which are already O(a)-improved in our setup, cf. Appendix C. (L1, M, 0). As explained in the text, these quantities are expected to vanish in the static limit 1/z = 0.
Then it is interesting to define the 1/m h contribution to Φ 3,4,5 , 4) and to study the mass behavior of these quantities, where Φ

Matching in L 1
With the set of simulations described in the previous section, we have computed η(L 1 , a) and ϕ(L 1 , a). Thus the HQET parametersω(M, a) can be obtained from eq. (2.5), viz.
However, in practice we split this equation in the following waỹ where we have used eq. (3.4). Whether one uses eq. (2.5) or eq. (3.5) to defineω(M, a) only affects the way we treat the lattice artefacts, but one expects a better precision in the second case. The reason is that the quantities η 3 and η 4 are by far the dominant part of Φ 3 and Φ 4 (and Φ 5 vanishes in the static approximation). As explained in the previous subsection, they are static quantities which extrapolate to the continuum with O(a 2 ) corrections. On the contrary, η 1 and η 2 are divergent and have to be kept in the combination as in eq. (3.5). They are extrapolated linearly in a/L 1 . At this point we would like to comment about the separation of the different orders in the effective theory. The situation for the first two observables is different from the others, since η 1 and η 2 do not have a continuum limit. However, the whole matching procedure can be carried out at static order as well. In that case,ω 3 = c (1) A is just the improvement coefficient ac stat A which is approximated by perturbation theory. Moreover, the parameters ω 4 andω 5 are of order 1/m h , and therefore are set to zero. This setup defines the static approximation of the first two observables, Φ stat 1 and Φ stat 2 , with ϕ = diag(L 1 , 1). Performing the matching only for these observables (instead of Φ HQET i with i = 1, . . . , 5) allows us to determine the two parameters m stat bare and ln(Z stat A ) at static order.

Evolution to a larger volume L = L 2
We then consider a set of simulations of HQET in which we use the same parameters as in the previous step (HQET in volume L 1 ), but where we double the number of points in each space-time direction. There we compute the quantities ϕ(L 2 , a) and η(L 2 , a) with L 2 = 2L 1 .
We now have all the ingredients to compute Φ(L 2 , M, 0) using eq. (2.6), but for the reason given above, we use a slightly different version: At the static order, the matching procedure at L 1 and the evolution to L 2 can be done in the same way as with the five-component vector Φ. The result defines the quantities Φ stat i (L 2 , M, 0) for i = 1, 2, and their continuum extrapolation is shown in Fig. 5. When the matching is performed at the next-to-leading order, we have Φ stat i (L 2 , M, 0) = η i (L 2 , 0) for i = 3, 4, 5 and define the 1/m h -contributions as Their continuum extrapolations are shown in Fig. 6 and Fig. 7. The parameters ω(M, a) are computed by the relation analogous to eq. (3.5) in the volume L 2 In eq. (3.6) we have chosen to write the evolution of the observables to the volume L 2 in terms of the HQET couplingsω i . Equivalently we can introduce a matrix of step scaling functions. In order to do that, one substitutes in eq. (3.6) theω i by the matching eq. (2.5). One obtains an equation of the following form: The explicit definitions of D, Σ and Σ can be found in [14]. We have implemented the treelevel improvement of these step-scaling functions in order to obtain a smoother approach to the continuum limit. Our results are obtained in this way, but tree-level improvement actually only has a small influence on our results after extrapolation.

HQET parameters to be used in large volume simulations
As we have already mentioned in the text and explained in detail in the appendices, the HQET parameters are obtained at five values of the bare coupling g 2 0 , which are such that the renormalized couplingḡ 2 (L 2 /4) is kept constant for the different ensembles. This is done by setting β to some precise values: β = 5.2638, 5.4689, 5.619, 5.758, 5.9631. In order to be able to use the HQET parameters in large volume simulations -and make phenomenological predictions -we interpolate (or, in one case, we extrapolate) them to β = 5.2, 5.3, 5.5, 5.7, using a quadratic polynomial in β. We show this interpolation/extrapolation in Fig. 8. The results are reported in Table 1 for the action HYP2. All errors quoted in this table are statistical (obtained with a standard jackknife procedure), including errors coming from the QCD renormalization constants at finite lattice spacing. There is still an overall relative error contribution of about 0.9% from the quark mass renormalization in QCD, i.e., the factor h in eq. (B.5) of Appendix B, which relates the RGI mass to the SF running mass in the continuum limit. Since in practice this error will only become relevant when the z dependence of the HQET parameters is actually applied to interpolate to the B-meson scale and to extract physical quantities, it is enough to include it then. The complete set of results can be found on our website http://www-zeuthen.desy.de/alpha/, together with the relevant error-correlation matrices. From preliminary studies [18,24], we know that the physical b-quark mass corresponds approximately to z = 13. For this value of z we display our results for am bare rescaled by L 1 /a for each value of the lattice spacing in Fig. 9. The parameter m bare absorbs the power divergences present in the binding energy of an heavylight meson in HQET: a 1/a divergence at the static order, and a 1/a 2 divergence at the 1/m h order. As one can see from the figure, m bare is dominated by these power divergences. Clearly, in order to guarantee the existence of the continuum limit, these divergences have to be removed non-perturbatively, which is one of the benefits of our approach.
Finally, we would like to comment on the size of 1/m 2 h terms. In our setup, the fermion fields are periodic (in space) up to a phase. As a consequence, the static quantities depend on the choice of one angle, that we call θ 0 , and the quantities computed at the next-to-leading order of HQET depend on three angles: θ 0 , θ 1 , θ 2 . So far in this work we have considered only what we call the standard choice of θ angles [13,14], namely (θ 0 , θ 1 , θ 2 ) = (0.5, 0.5, 1). When we change the values of θ i we change the set of observables, and thereby the matching conditions. The static quantities computed for different values of θ are thus expected to differ by terms of order 1/m h . Once the 1/m h corrections are added, the difference should be of order 1/m 2 h . In general we find no significant θ dependence for all the HQET parameters ω i at the next-to-leading order of the effective theory, meaning that the 1/m 2 h terms are not visible within our statistical precision. As an illustration, we show the spread of our results for ln(Z stat A ) and ln(Z HQET A ) in Table 2 for the discretization HYP2, β = 5.3, and for z = 13. Due to statistical correlations, some of the errors on these 1/m 2 h terms are significantly smaller than our errors for ln(Z stat A ) and ln(Z HQET A ) themselves, but still no significant 1/m 2 h term is found.

Conclusions
We have reported on the first unquenched determination of a set of HQET parameters. While other approaches treat the dependence of the effective parameters on the renormalized coupling constant in a perturbative fashion, our approach is entirely non-perturbative in the QCD coupling constant and avoids the difficult-to-estimate uncertainties of perturbation theory for this case [17,25] 5 . Furthermore, our matching procedure takes into account non- 5 Ref. [25] computes the matching of the static-light currents including O(α 3 s ) and comments on a bad behavior of perturbation theory. In Section 3.3.2 of [17], it is shown that different choices for the renor-     Table 1: HQET parameters as a function of the bare coupling for the action HYP2 at z = 11, 13, 15. The central value z = 13 is close to the physical b-quark mass. The first two parameters, m stat bare and ln(Z stat A ), result from the matching at static order, while the remaining entries are the parameter set resulting from the matching of HQET at next-to-leading order.     perturbatively the power-law divergences of the effective theory which can be numerically dangerous if simply addressed by relying on perturbation theory. We have shown that these divergences get absorbed into the effective parameters of the theory. The results presented here can be combined with hadronic matrix elements and energies computed in large volume simulations (some preliminary results have been reported in [18,24,26]). We look forward to presenting our determination of the b-quark mass, of the B-meson spectrum and of heavylight decay constants at the 1/m h order of HQET in the two flavor theory, which will use the parameters computed in this work. Concerning the decay constant, the precision of the current matching is about 0.5% in the static approximation and 3% to first order in 1/m h . The latter puts a mild restriction on the available precision for the decay constant. It is reassuring to see that in the cases where we can check for contributions of 1/m 2 h terms explicitly, these are considerably below our precision. This is in line with the strong hierarchy that we observe between different orders in the HQET expansion. We finally note that our numerical computations have been carried out on apeNEXT computers, decommissioned by now. A future application of our matching strategy with three or more flavors can be expected to reach a further improved numerical precision on present or future hardware. HE 4517/2-1. N. T. acknowledges the MIUR (Italy) for partial support under the contract PRIN09. We thank NIC/DESY and INFN for allocating computer time on the apeNEXT computers to this project as well as the staff of the computer centers at Zeuthen and Rome for their valuable support.

Appendix A Observables
For completeness, we remind the reader of the definitions of the observables Φ introduced in [14], where the quantities are built from the (renormalized) SF boundary-to-bulk correlation function f A of the pseudoscalar channel and the boundary-to-boundary correlator f 1 and k 1 of the pseudoscalar and vector channel, respectively. The definitions are such that, at the 1/m h order of the heavy quark expansion, the observables assume the following form which is just an explicit version of eq. (2.3). The precise definitions of the HQET quantities can be found in [14].

Appendix B Tuning of L 1 and renormalization in finite-volume QCD
As explained in the main text, the basic element of our non-perturbative strategy to compute the HQET parameters consists in imposing matching conditions between a set of renormalized finite-volume observables in QCD, extrapolated to the continuum limit, and their counterparts in HQET, expanded up to order 1/m h . Since in this way the HQET parameters get determined by feeding results from non-perturbatively renormalized QCD into the effective theory, we summarize here how the renormalization in QCD (i.e., of the gauge coupling, the quark masses and the relevant composite fields) is performed. We work in a small volume of linear extent L 1 ≈ 0.4 fm, where the SF setup [27,28], with T = L and θ = 0.5 as the periodicity angle of the sea quark fields, serves as our finite-volume renormalization scheme. The physical volume stands in one-to-one correspondence to the SF gauge couplingḡ 2 (L) running with the scale L [29,30]. We thus define L 1 by the condition The known step scaling function of the SF coupling in two-flavour QCD [30] implies According to eq. (B.1), β is determined by requiringḡ 2 (L 1 /2) = 2.989 for given resolutions 2a/L 1 . This peculiar value of the SF coupling was indicated by an initial simulation with L 1 /(2a) = 20, while for 10 ≤ L 1 /(2a) ≤ 16 additional simulations and interpolations in β, based on the known dependence of the SF coupling and the sea quark mass on the bare parameters available from the data of [30], were employed for fine-tuning to that target. Employing the non-perturbative β-function of the SF coupling and estimates from our quenched calculation on the effect of propagating uncertainties inḡ 2 (cf. Appendix D of [13]), we can assess an uncertainty of about 0.05 inḡ 2 to translate via L 1 into an uncertainty in the b-quark mass of at most 0.5%, which is negligible compared to the present direct uncertainty in the quark mass renormalization discussed below.
In the quark sector, the sea and light valence quark masses are taken to be the same; in the numerical computation they are actually tuned to zero. The condition m l (L 1 /2) = 0 is met by setting κ l to the critical hopping parameter, κ c , which is determined by the vanishing of the PCAC mass of the sea quark doublet, defined as in eq. (2.17) of [30] through the O(a)-improved axial current in the SF setup with boundary field "A" [29]. Again, κ c was estimated and partly fine-tuned on basis of the data published in [30], whereby for the improvement coefficient of the axial current, c A , 1-loop perturbation theory [32] and the non-perturbative estimates of [33] (after they had become available) were used. A slight mismatch of | L 1 2 m l (L 1 /2)| < 0.05 of this condition is tolerable in practice. The resulting triples (L 1 /a, β, κ l ) are collected in the three leftmost columns of Table B.1. Note that the mentioned 1-loop value for c A was only used in the preliminary determination of κ l . All PCAC masses listed here are computed with the non-perturbative c A of [33].
It remains to fix the renormalized mass of the heavy valence quark to a sequence of values, fairly spanning a range from around the charm to beyond the bottom quark mass. To this end we choose the dimensionless variable z ≡ LM , with M being the RGI mass of the heavy valence quark, because the latter is related via  Table B.1: Bare parameters (L/a, β, κ l , κ h ) used in the computation of the heavy-light QCD observables for L = L1. As explained in the text, they are fixed by the renormalization conditionsḡ 2 (L/2) = 2.989, L 2 m l (L/2) ≈ 0 and LM = z for the SF coupling and the light PCAC and heavy RGI quark masses, respectively. The entering renormalization constants ZP and Z and the improvement coefficient bm were calculated in [34].
to the subtracted bare heavy quark mass, m q,h , and its hopping parameter, κ h , in the O(a)improved theory. Here,  [35], and the renormalization factor Z and the improvement coefficient b m from [34]. The scale dependent renormalization constant Z P for the specific β-values in question was extracted in [34], following exactly the definition of [36]. Z P , b m and Z are also listed in Table B which represents the universal, regularization independent ratio of the RGI heavy quark mass to the running quark mass, m, in the SF scheme at the renormalization scale µ. h(L 1 /2) was evaluated by a reanalysis of available data on the non-perturbative quark mass renormalization in two-flavour QCD as published in [36].
Given the values 6, 7, 9, 11, 13, 15, 18, 21} (B.6) of the dimensionless RGI heavy quark mass in L 1 and resolutions L 1 /a = 20, 24, 32, 40, eqs. (B.3) -(B.5) can now straightforwardly be solved for the corresponding nine heavy valence quark hopping parameters κ h = κ h (z, g 0 ) that fix z to the numbers in eq. (B.6). 6 These hopping parameters and the associated z-values are collected in the two rightmost columns of  Let us still comment on the error budget arising from the above procedure of fixing z [34], which has to be accounted for in any secondary quantity analyzed as a function of z. From the uncertainties on Z A , Z P , Z, and b m quoted in the respective references [34,35] (see also  Table B.1), one obtains by the standard rules of Gaussian error propagation an accumulated relative error on z in the range 0.38% ≤ (∆z/z) ≤ 0.41% for all z-values and lattice resolutions in use here. The contribution from the universal continuum factor h(L 1 /2), eq. (B.5), represents with ∆h/h = 0.92% the dominating source of uncertainty in the total error budget of ∆z/z = 1.01%. Note, however, that the error on this universal factor h has to be propagated into the QCD observables Φ i only after their extrapolations to the continuum limit. It is not included in the errors in Table B.1. A similar tuning procedure has to be performed for the HQET simulations necessary for the matching. Since the matching proceeds through renormalized quantities and therefore in the continuum limit, there is no need to use the same lattice resolutions on both the HQET  Table B.2. From the two runs at L/a = 8 we estimate s = ∂ḡ 2 ∂z ḡ 2 =4.484 = 1.4(4), which we use to set a bound on the quark mass m l (which should be vanishing), such that its effect on the coupling is below the statistical error. The L/a = 12 * label refers to an additional run used to propagate the error onḡ 2 into the HQET observables. This is done by computing all HQET observables (at T = L and T = L/2) also at the bare parameters of 12 * in order to estimate the variation of our primary HQET observables with respect to a variation of the renormalized coupling. This procedure neglects the lattice spacing dependence of the variation, which is justified for an error computation. In practice, the uncertainty onḡ 2 (L) is the dominating piece of the errors of η 3 , η 4 shown in Fig. 2, while it can safely be neglected for all other observables within our present error budget.

Appendix C Simulation details
Our computations are in a natural way split into two parts, the generation of N f = 2 gauge field ensembles at the tuned parameters (L/a, β, κ sea ≡ κ l ), and the subsequent computation of all relevant correlation functions. The Sheikoleslami-Wohlert improvement coefficient c sw (g 0 ) is set to its non-perturbative estimate [19] for N f = 2. In order to have an improved action in the SF we include boundary counterterms to cancel boundary-induced lattice artefacts. The corresponding improvement coefficients, c t (g 0 ) and c t (g 0 ), are always set to their known 2-loop [37] and 1-loop [32] values, respectively. This guarantees that any boundaryto-boundary correlation function, such as f 1 , is O(a)-improved.
Ensemble generation: For simulating a doublet of mass degenerate, non-perturbatively improved dynamical Wilson fermions in the SF (θ sea = 0.5), we use the algorithmic implementation described in detail in [38] to produce the ensembles needed for the matching in the volume L 1 . Applying the step scaling technique in HQET to volume L 2 = 2L 1 also requires simulations at resolutions a/(2L 1 ) while keeping all other parameters fixed. Furthermore, our strategy involves simulations at fixed time extent T = L as well as T = L/2.
Since most of the simulations used for the HQET observables are fast and can usually be performed using several replica, we aimed for a total statistic of at least 8000 configurations in these ensembles. Only for our most expensive HQET ensembles with L 2 /a = 32 we did not reach this goal due to limited computing resources and thus restricted ourselves to have O(3000) configurations here. The QCD simulations have larger values of L/a (and thus smaller lattice spacings) compared to the HQET simulations. It is therefore more difficult to achieve high statistics for the QCD ensembles. Even more so since the gap in the spectrum of the Dirac operator, which allows to simulate at vanishing quark mass in the SF, decreases proportionally to the inverse time extent. Hence, for the production of ensembles used to measure QCD observables, our goal was just to reach a reasonable statistics, and thereby to obtain small and comparable errors in our final observables at the different resolutions. Thus, the ensemble size roughly increases from O(500) at L/a = 20 to O(1500) at L/a = 40.
Using the notation introduced in [38], we list the relevant algorithmic parameters and additional details in Table C.1. For QCD the bare parameters (L/a, β, κ l ) are those of Table B.1, while for HQET we use those in Table B.2 which are not marked by a star. The molecular dynamics (MD) is characterized by specifying the trajectory length, the integrator, and the step size(s). For the trajectory length we choose τ = 2 in MD units since we expect autocorrelation to be reduced [39] in that case. As integration scheme we always use multiple time scales with leap-frog integrator, also known as Sexton-Weingarten scheme [40]. The tunable algorithmic parameter ρ 0 is introduced as a mass-preconditioning of the Diracoperatorà la Hasenbusch [41]. The step sizes for the corresponding two pseudofermions are δτ 0 and δτ 1 , while for the gauge force we use δτ 0 /δτ g = 4 throughout. N (i) CG is the average number of conjugate-gradient iterations used to solve the symmetrically even-odd preconditioned Dirac equation during the trajectory, and P acc is the acceptance rate of the simulation. τ meas gives the MD time between configurations which have been stored on disk and used for measurements. In case we list more than one value of τ meas , we have performed independent simulations with different measurement frequencies which have been chosen to be of the typical size of observed integrated autocorrelation times τ int . Explicitly, we show the average plaquette and PCAC mass of the production runs together with their estimated autocorrelation time in Table C.2. We also list the results for the pseudoscalar boundary-toboundary correlation function f 1 , which typically is the quantity with the largest integrated autocorrelation time among the different SF correlation functions.
Measurements: since there is no explicit heavy quark mass contribution to correlation functions in HQET, we just need to specify κ l ≡ κ sea and the static quark action(s) in use to compute static-light correlation functions. The latter have been computed using the two static actions HYP1,2 described in [22]. For measurements in QCD we compute heavy-light observables in a partially quenched setup with κ l ≡ κ sea and heavy valence quark hopping parameters κ val,h ≡ κ (i) h , i = 1, . . . , 9. The latter slightly deviate from those in Table B.2 because at the time we fixed the values of z the updated g 2 0 -dependence of Z A [35] entering eq. (B. 3) was not yet at hand. This translates into a small mismatch ( 0.4%) of the z-values where we did the computation with respect to our target z-values. The QCD observables were interpolated to the target z-values to account for this.