Non-perturbative tests of continuum HQET through small-volume two-flavour QCD

We study the heavy quark mass dependence of selected observables constructed from heavy-light meson correlation functions in small-volume two-flavour lattice QCD after taking the continuum limit. The light quark mass is tuned to zero, whereas the range of available heavy quark masses $m_h$ covers a region extending from around the charm to beyond the bottom quark mass scale. This allows entering the asymptotic mass-scaling regime as $1/m_h \to 0$ and performing well-controlled extrapolations to the infinite-mass limit. Our results are then compared to predictions obtained in the static limit of continuum Heavy Quark Effective Theory (HQET), in order to verify non-perturbatively that HQET is an effective theory of QCD. While in general we observe a nice agreement at the few-% level, we find it to be less convincing for the small-volume pseudoscalar decay constant when perturbative matching is involved.


Introduction
Heavy quark systems, notably B-and B s -mesons, provide a unique opportunity to perform stringent tests of the Standard Model and probe signals of new physics. To maximize the impact of experiments performed at the Large Hadron Collider, for example, it is imperative to control various aspects of the underlying theory. In particular, strong interaction effects must be understood at the quantitative level, including reliable estimates of systematic uncertainties.
Although in principle lattice QCD in a large physical volume L 3 allows for ab initio computations of hadronic matrix elements and energy levels, the presence of heavy and light quarks still renders computations very demanding. Current state-of-the-art lattice simulations of QCD with N f ≥ 2 dynamical quarks usually reach lattice spacings down to a ∼ 0.05 fm, while satisfying m π L 4 in order to keep the finite-size effects under control. Due to increasing computer power and algorithmic advances over the last decade, it has become possible to simulate close to the physical pion mass for the figures just given. However, additionally including relativistic b-quarks with their "heavy" physical mass of about 4 GeV requires very fine lattice spacings a 1/m h ∼ 0.05 fm in order to monitor its discretization effects in the spirit of Symanzik's local effective theory. Associated with this is the problem of large scale separations, m B /m π = O(100), to be accommodated in a single simulation, which remains out of reach in the near future. Therefore, one still needs to take a detour to an effective description for the heavy quark. Here we focus on the Heavy Quark Effective Theory (HQET) [1][2][3][4] which provides a natural framework to study heavy-light mesons through a systematic expansion in the inverse heavy quark mass, 1/m h . A non-perturbative implementation of HQET on the lattice [5], including the nextto-leading order in the 1/m h -expansion, has been tested and applied successfully for the quenched case in the past [6][7][8]. It requires to solve a set of matching relations between quantities in continuum QCD and lattice HQET in a small physical volume. In two-flavour QCD, the resulting non-perturbative set of HQET parameters [9] was recently used to extract phenomenologically relevant parameters such as the b-quark mass, the B-meson decay constant and hyperfine splittings from large-volume simulations [10][11][12]. At present there are efforts to extend the matching strategy to include the vector meson channel [13,14] in order to compute f B * and form factors of semi-leptonic B (s) → π (K) transitions.
In this paper we probe predictions of HQET by studying the asymptotic behaviour 1/m h → 0 of continuum-extrapolated lattice QCD observables, computed in a small volume (L ≈ 0.4 fm) with N f = 2 dynamical flavours of non-perturbatively O(a)-improved massless Wilson fermions. Extrapolations to the static limit are performed, which not only gives numerical evidence of the correctness of static order HQET but also allows assessing the size of the 1/m h -effects. As such, the present study constitutes a non-trivial non-perturbative test of HQET being an effective theory of QCD and extends earlier work in quenched QCD [15] to the physically more realistic situation with dynamical light quarks. In particular, we investigate a considerable set of observables at varied kinematics, which in the spirit of the general non-perturbative matching strategy mentioned above are usually employed to define suitable matching relations. Complementary to perturbative studies [13,14,16], this yields non-perturbative insights into their heavy quark mass asymptotics and provides criteria for which of them are to be preferred within the matching strategy of [5], such as impact of mass-dependent cutoff effects in the continuum limit extrapolations, numerical accuracy and magnitude of higher-order corrections in 1/m h .
In contrast to a non-perturbative matching of (lattice) HQET to continuum QCD, the perturbative approach relies on a perturbative evaluation of matching (resp. conversion) functions, often called Wilson coefficients. To properly recover the static limit in HQET as m h → ∞, also our QCD observables -non-perturbatively evaluated at finite quark mass -still have to be combined with matching functions of this kind. These are only perturbatively known, up to three-loop order in most cases. To disentangle in our tests the genuine non-perturbative properties of the theory encoded in the observables from perturbative effects induced by the conversion functions, we map out their m h -dependence towards the static limit for conversion functions of different perturbative orders. This comparison gives a rough idea on the systematic error that is involved, when the matching between HQET and QCD is performed perturbatively. As will be exposed by the example of the pseudoscalar decay constant below, we observe that -with perturbative matching at work -the agreement between the large-mass QCD asymptotics and the HQET prediction may not be as good as expected, even if a three-loop expression for the conversion function is used. We take this as an indication that in an effective theory framework for heavy quarks the matching should be done non-perturbatively, if one wants to have the systematic errors under control. In addition, we consider QCD observables, which have a non-trivial static limit but do not depend on any conversion function. Hence, their m h → ∞ limit is much less affected by systematic uncertainties such that they are among the cleanest observables in our non-perturbative tests of the effective theory approach to heavy-light physics.
Some preliminary results on a smaller subset of observables, data ensembles and statistics have already been reported in [17,18].

Observables
This paper follows up our previous work [19] on the definition of a line of constant physics in N f = 2 lattice simulations. It allows us to non-perturbatively study the quark mass dependence of relativistic QCD meson observables in a finite box of extent L 1 ≈ 0.4 fm [20]. We explore a wide range of quark masses that starts below the charm sector and goes beyond the bottom quark region. This is done in a partially quenched setup, i.e., the light quark mass is set to the approximately vanishing mass of a degenerate sea quark doublet and all other quarks are quenched. We generically refer to the latter as heavy quarks of mass m h . Equivalently, we will assign to them the dimensionless mass parameter z = L 1 M h from now on, where M h ≡ M denotes some fixed value of the renormalization group invariant (RGI) heavy quark mass. In [19] we have already shown that in such a small box the lattice spacing can be chosen small enough so that all heavy quarks up to a certain value can be simulated relativistically while keeping cutoff effects in the O(a) improved theory well under control. For any unexplained notation, the reader may consult [15,19].
Since our interest lies in relating predictions made by HQET non-perturbatively to the proper counterpart in QCD towards the limit 1/m h → 0, we furthermore take into account measurements of HQET observables that have been done in the framework of a general non-perturbative matching strategy of HQET and QCD in the very same volume L 1 , but to a much higher statistical accuracy. Additional details can be found in appendices B and C of reference [9]. Working in a finite (and small) volume, all matrix elements and energies become effective quantities which intrinsically depend on the scale L 1 . For notational brevity we often suppress this dependence in the following.
Our main observables are built from Schrödinger functional (SF) correlation functions [21,22] in a T × L 3 volume with T = L = L 1 fixed and periodic boundary conditions in space. The fermion fields are taken periodic only up to a phase, where we use θ ∈ {0, 0.5, 1}. In correlation functions, this periodicity angle θ amounts to a projection onto quark and antiquark momenta with components ±θ/L. In time direction, Dirichlet boundary conditions are imposed at x 0 = 0, T where source quark and antiquark fields are separately projected onto vanishing spatial momentum. We are interested in the pseudoscalar (PS) and vector (V) channel using finite-volume heavy-light QCD currents.
They are given by the time component of the axial vector and spatial components of the vector current, respectively: To be more precise, we use their O improved [23,24] and non-perturbatively renormalized [25] lattice versions. Furthermore, we need the pseudoscalar heavy-light current the renormalization factor of which, Z P (µ, g 2 0 ), has been determined non-perturbatively in the SF scheme at scale µ = 1/L 1 during the production runs reported in [9]. Since they have not been quoted in that reference we list them together with further details in table 4.

Definitions
How to non-perturbatively set up a line of constant physics in the envisaged small volume L 1 ≈ 0.4 fm has already been reported in [19]. In that volume we have four different ensembles with non-perturbatively O improved dynamical Wilson fermions made of a doublet of massless quarks. The range of lattice spacings used is 0.01 fm a 0.02 fm. This ensures feasible continuum limit extrapolations of the QCD observables to be introduced below, which depend on the dimensionless RGI heavy quark mass fixed to 3, 3.3, 4, 6, 7, 9, 11, 13, 15, 18, 21} . (2.4) In contrast to our previous work, we added four additional z-values at the lower end to also cover the charm quark region. Details about the latter, which were needed to perform the additional measurements, are listed in table 4. In the computation of HQET observables we can naturally rely on larger lattice spacings (0.025 fm a 0.067 fm). We use the "HQET[L 1 ]" lattices with T = L = L 1 as specified in table C.1 and C.2 of ref. [9]. For completeness we list some additional details in table 5 that were not published there.
The HQET observables themselves are computed using two different static actions, referred to as HYP1 and HYP2 [26], in order to have an improved noise-to-signal behaviour compared to the original Eichten-Hill action. For the present purpose of testing HQET we want to express all observables in terms of matrix elements computed in a finite volume. The interested reader can find the corresponding notation in terms of the traditional SF correlation functions in appendix A and reference [15]. In passing we just note that only states, which are eigenstates of spatial momentum with eigenvalue zero, enter them.

Effective masses
The small-volume effective pseudoscalar and vector meson mass in terms of Hilbert space matrix elements read, adopting an operator notation for the quark bilinear composite fields, where |Ω(L) ≡ e −x 0 H |ϕ 0 (L) denotes the vacuum state given in terms of the Hamiltonian H and the SF intrinsic vacuum boundary state |ϕ 0 (L) at x 0 = 0. Having in mind a variation of the heavy quark mass for our non-perturbative tests later on, we denote a general heavy-light pseudoscalar state by |B(L) ≡ e −x 0 H |ϕ B (L) and a heavy-light vector state by |B * (L) ≡ e −x 0 H |ϕ B * (L) . Both are again given through the time evolution operator and the well-defined boundary states |ϕ X with quantum numbers in the respective channel, X = B, B * . Since all states naturally depend on the finite volume (or box) size L and are also considered as functions of the RGI heavy quark mass M (or z), we drop these dependencies for the moment to ease notation, as we do so for their additional dependence on the SF-specific periodicity angle θ of the fermion fields. Here and from now on, we fix x 0 = T /2 ≡ L/2 in the correlation functions employed to construct the observables above and those to be introduced in the subsections below. It is then worth to emphasize that the time evolution operator e −T H/2 suppresses high-energy states exponentially such that |Ω(L) and |B(L) , upon expanding them in terms of eigenstates of H, are dominated by contributions from states with energies of at most ∆E = O(1/L) above the ground state. Therefore, as HQET is expected to apply to correlation functions at large Euclidean time separations, it particularly describes their large-mass behaviour at large Whereas physical masses must be computed in large-volume simulations, the definition of our observables is such that they agree with the physical ones in the large-volume limit L → ∞.

Decay constants and ratios
Furthermore and analogously, we define the following QCD observables, suppressing again their dependence on M , L and θ: where Y PS and Y V are the finite-volume heavy-light pseudoscalar and vector decay constant, respectively. As L → ∞, they become proportional to the physical heavy-light pseudoscalar and vector meson decay constants. Accordingly, Y PS/V is the ratio of the two, which in large volume becomes proportional to f B /f B * , if the heavy quark is set to the b-quark. The ratio R spin is proportional to the spin splitting between the pseudoscalar and vector channel and, as predicted by HQET, has to vanish in the static limit owing to the heavy-quark spin symmetry. These observables also involve boundary-to-boundary SF correlation functions, see appendix A, which properly cancel the multiplicative renormalization factors of the boundary quark fields.
All quantities involving [•] R are understood to be renormalized non-perturbatively, while for others such factors either drop out or are not needed at all, such that alltogether they thus are finite and possess a well-defined continuum limit.

Quantities with different kinematics
The SF is especially useful, if one wants to probe physics with different kinematics. Here we do so by changing the fermionic phase angle θ as mentioned earlier. Whereas in the last section all observables were meant to be evaluated at the same values, i.e., θ 0 ∈ {0, 0.5, 1}, we now turn our attention to quantities that are made of two matrix elements of heavy-light composite fields referring to fermionic periodicity phases different from eachother. With the same notational conventions as before they read where in our actual calculations we consider the following pairs of phase angles: It is important to note that all multiplicative renormalization and improvement factors cancel in these ratios. 1 In this respect and because one expects cancellations of cutoff effects at every fixed value of z, these QCD observables (as well as their counterparts in HQET) may be seen as "gold plated" observables for the purpose of testing the asymptotic behaviour of heavy-light physics as z → ∞ for different kinematical setups.

Observables in HQET
After the continuum limit of the previously defined QCD observables has been taken, we aim for an extrapolation to the static limit, 1/z → 0, in order to compare their asymptotic behaviour to the one predicted by HQET in the continuum. According to the systematic heavy quark expansion, all quantities approach a well-defined value in that limit. Classically, the leading asymptotic behaviour of our effective masses (2.5), for instance, is linear in the heavy quark mass z, but it receives logarithmic modifications on the quantum level owing to the scale dependent renormalization of the effective theory to compare with, viz.
where C mass (z) denotes the conversion function that relates the heavy quark's pole mass to the RGI heavy quark mass M = z/L 1 .
1 They differ from 1 (or 0 for R1) only due to θ1 = θ2, both at finite z and in the static limit.
A generic conversion function C X (z) carries all the logarithmic dependence of a given quantity X to some order in perturbation theory such that only power corrections in 1/z remain in the effective theory. For more details, we refer to appendix B and appendix B of [15]. To avoid a remnant renormalization scheme dependence in the (static) effective theory, we favour to fully express the asymptotic behaviour in terms of renormalization group invariants (RGIs). For the quantities of section 2.3, i.e., X ∈ {PS, V, PS/P, PS/V, spin}, this means The ratios R PS/P , R PS/V and Y PS/V , along with the associated C X , approach 1 in the static limit of HQET, whereas the effective decay constants Y PS and Y V both approach the finite-volume RGI static-light decay constant X RGI (L) in this limit, as a consequence of the heavy-quark spin symmetry. In the two-flavour theory at hand [27], the renormalization scale was implicitly fixed by a value for the renormalized SF coupling ofḡ 2 (µ)| µ=L −1 max ≡ 4.61. However, since for the purpose of this study we are working at a slightly different physical volume L 1 L max , the universal part Z stat A,RGI /Z stat A (µ) of the total renormalization factor, which relates a matrix element of the static axial current renormalized at a scale µ, X R (µ), to the RGI one, had to be re-evaluated. The outcome for µ = L −1 1 based on the data of [27] is which allows us to directly compute the (renormalization scale and scheme independent) quantity X RGI instead of the renormalized (and thus scale dependent) static-light decay constant X R (µ) in finite volume. Some additional technical details on X RGI and its error budget are postponed to appendix C. By contrast, the RGI matrix element of the spin splitting operator, X spin RGI , is not known, because the corresponding RG running is not available for the dynamical flavour theory. Note that in principle it would be possible to extract it from our QCD data according to the given asymptotic behaviour, (2.16), but with the only perturbatively known conversion function C spin we do not expect the result to be particularly meaningful.
Our most stringent tests of HQET will finally arise from the QCD observables defined in section 2.4. In these quantities, the conversion functions cancel such that no logarithmic corrections are left to all orders in perturbation theory and non-perturbatively. Hence, all uncertainties from the matching between QCD and HQET are absent and one just faces the genuine power corrections in 1/z of the effective theory. This makes the comparison of continuum QCD in the limit 1/z → 0 and continuum HQET at static order entirely non-perturbative and very well controllable, without encountering any systematic errors induced by the perturbatively evaluated C X . Exploiting again the heavy-quark symmetry, the HQET counterparts of these QCD observables in the static limit are:

Results
In the following subsections we first discuss exemplary continuum extrapolations for some of our test observables before we turn our attention to the main results, i.e., their extrapolations as 1/z → 0 to the static limit of HQET. Additional details about our continuum extrapolations can be found in appendix C.

Representative continuum extrapolations
For the quantities of subsections 2.2 -2.4, which as properly renormalized QCD observables are now generically denoted by Ω QCD = Ω QCD (L, M, a), we perform extrapolations to the continuum limit (CL) using a global fit ansatz in order to have better control over massdependent lattice artefacts. Due to the latter, we exclude some points at coarsest lattice spacings from these global QCD continuum extrapolations. In general we aim at taking the continuum limit of a QCD observable according to the global fit ansatz Ω QCD (L, z, a) = Ω QCD (L, z) 1 + (a/L) 2 · ρ 0 + ρ 1 z + ρ 2 z 2 , (3.1) which accounts for terms proportional a/L × aM and (aM ) 2 . From earlier studies [15,28] it is known that an exclusion limit of aM > 0.7 has to be imposed on the data entering in the extrapolating fits, in order to avoid contaminations which are potentially dangerous for a reliable continuum limit. To assure stability in our CL results, different fit ansaetze have also been studied (e.g., allowing for a cubic term in a/L with mass-dependent coefficient or omitting some of the lightest masses, say z ≤ 4); these lead to consistent results. As the general scaling behaviour towards the continuum limit looks rather similar among the different observables, we only show two representative examples here. In the left panel of figure 1 we present the result for the effective pseudoscalar decay constant at θ 0 = 0.5 and in the right the outcome for a ratio of the same quantity evaluated at different kinematical parameters, namely at (θ 1 , θ 2 ) = (0.5, 1). As yet we have not taken into account the error stemming from the tuning of the heavy quark mass (2.4) at finite lattice spacing. The uncertainty of the continuum heavy quark mass M = z/L 1 is ∆M/M = ∆z/z = 1.01% [19]. We add its error quadratically, before performing any extrapolations to the static limit as they are presented in the following subsections. The derivative is estimated numerically from the data at hand. Its contribution to the total error budget is actually negligible, as can be inferred from the continuum mass dependence displayed in figure 1, for instance. In appendix C we list a representative selection of results at finite lattice spacing and its continuum limit. In some cases, continuum extrapolations of associated observables in HQET have also to be performed. Taking the continuum limit of observables in the static theory is much more straightforward, because there is obsviously no dependence on z and thus no massdependent cutoff effect that needs to be controlled in addition. However, they depend on the two static actions employed here (HYPi, i = 1, 2), and hence this suggests to adopt a joint continuum extrapolation of a given HQET observable according to As two explicit examples we show in figure 2 the continuum extrapolation of R stat X (θ 1 , θ 2 ) and R stat PS (θ 1 , θ 2 ). The numerical results are listed in table 2 and 3, respectively. Before turning our attention to the main results, we want to add some general remarks on their presentation and the analysis underlying them.

Results including perturbative conversion functions
As emphasized before, for certain (continuum-extrapolated) observables Ω QCD (L, M, 0) ≡ Ω QCD (L, z) we need to take into account logarithmic corrections when comparing HQET with QCD in the heavy quark mass limit. Hence, we divide by the corresponding conversion function C Ω (M/Λ MS ) and, in a few cases where convenient, cancel the leading mass dependence by an appropriate multiplication with some power of z, such as = −1 in the case of Ω = LΓ, for instance. An exactly -or even non-perturbatively -known conversion function by definition removes all logarithmic contributions, while the remainder can then be assumed to be organized as a "power series" in 1/z, resembling continuum HQET in the asymptotic regime of large z. We thus perform extrapolations to the static limit of HQET of quantities with perturbative conversion functions according to the fit ansatz where Ω [0] represents the static limit of the observable in question as extracted through this fit from the relativistic QCD data. (The choice 1/z 0.26 is motivated below.) In practice, we must rely on a perturbative evaluation of the conversion functions up to limited order that have uncertainties decreasing only logarithmically (see, e.g., ref. [29]). Since those have to be combined with our non-perturbative lattice data, we generically cannot disentangle logarithmic and power-like contributions exactly up to some definite value of 1/z. A safe statement that can be made, however, is that one asymptotically expects eq. (3.4) to yield a good description of the data in the sense that the l.h.s. of (3.4) approaches Ω [0] in the limit z → ∞. To what extent a quantity Ω [0] agrees with its HQET counterpart is subject of the discussion below.
To represent our data, we employed several unconstrained (and weighted) polynomial fits in 1/z along (3.4) by varying the range of points that enter the fit with the degree of the polynomial. Even though the full data set -down to masses of ∼ 70% of the charm quark mass -is well described by a quadratic interpolation within our statistical accuracy, it is not surprising that this variation, depending on the observable under consideration, has a visible effect on the result of the static extrapolation. More precisely, we observe the clear general trend (motivating the z-range choice in (3.4)) that a fit of degree n should include data in the interval [0, n · ∆z −1 ], with ∆z −1 ≈ 0.13, in order to keep the static limit unchanged (within errors) and thereby independent of the employed polynomial fit ansatz. This statement holds separately true for the data incorporating C Ω evaluated at two-or three-loop order in perturbation theory. However, as will become evident in the sequel, the dominantly significant, systematic effect on the static extrapolation results originates from the only perturbative knowledge of the conversion functions that are required to relate QCD observables at finite z to their HQET counterparts, before the static limit is taken.

Asymptotic behaviour of effective meson masses
First, we study the large quark mass behaviour of the effective meson masses (2.5) according to its leading asymptotics in this regime that is given by eq. (2.10). The raw data and the continuum-extrapolated values for representative θ-values are collected in tables 6 and 7 of appendix C. Following eq. (3.4), we combine the effective meson masses with their associated conversion function C mass and remove the leading heavy quark mass dependence. The remaining finite piece of LΓ X (L, M, 0) zC mass (M/Λ MS ) then should approach 1 in the static (i.e., 1/z → 0) limit for all X = PS, V, P. Figure 3 confirms this expectation, where for each individual observable and θ 0 ∈ {0, 0.5, 1} we depict the remaining (i.e., subleading) asymptotic behaviour of our continuum data points for the effective meson masses using C mass evaluated at two-and three-loop order, cf. eq. (B.7). For ease of presentation, only the statistical errors of LΓ X (z, θ 0 ) are included in the figures. With the fit ansatz (3.4) and all errors taken into account, the results in the static limit, Ω [0] X , agree for all X = PS, V, P with 1 within errors indeed. In order to better judge the validity range of the asymptotic 1/z-expansion reflected by the extrapolating fits to the static limit -as well as the size of possible particular systematic effects at the physical scale of the b-and c-quark -here and in subsequent figures we add vertical error bands corresponding to z b = L 1 M b ≈ 13.25 [10] and z c = L 1 M c ≈ 3.04 [30].
Another non-trivial test consists in studying ratios of masses, in which the leading asymptotics drops out completely. Such a case is displayed in the bottom-right panel of figure 3. This is a first explicit example with the heavy-quark spin symmetry at work, according to which HQET predicts lim 1/z→0 [Γ PS /Γ V ] = 1. The deviation from one at the b-quark scale (z = z b ), dominated by the spin-splitting term, is about 2%. With increasing 1/z, higher-order terms appear to become relevant and contribute with opposite signs, leading to a deviation of about 20% at the charm quark mass scale. At least for Γ PS /Γ V , this supports the general expectation that HQET does not provide an all too good description for charm physics any more.

Asymptotic behaviour of ratios of heavy-light currents
Owing to the definitions in subsections 2.2 and 2.3, the finite-volume continuum QCD observable Y PS/V = Y PS /Y V approaches in the large-volume limit a combination of ratios of pseudoscalar to vector heavy-light meson decay constants and masses, which becomes of  phenomenological relevance at the b-quark mass scale: Thus it is interesting to also inspect the 1/z-dependence of Y PS/V . As before, after extrapolating to the continuum limit and accounting for the proper full (logarithmic) mass dependence via attaching the function C PS/V (z) to two-, three-and even four-loop accuracy, cf. eq. (B.6), we obtain Y PS/V (z, θ 0 ) and the corresponding extrapolations to the static limit of HQET reproduced in the left panel of figure 4. Again, we find that for every θ 0 the static HQET prediction (= 1) is reached within errors, with an almost linear approach as 1/z → 0 for z > 10; its slope grows with the flavour-twisted momentum θ 0 .
As discussed in appendix B, the conversion formula for the ratio of decay constants f B /f B * , C PS/V , is even known to four-loop accuracy [31], because in the entering difference of anomalous dimensions the unknown four-loop anomalous dimensions of the currents themselves drop out. It was already noted in [29,31] (and is also reflected in the middleright panel of figure 10 in appendix B) that the perturbative expansion of C PS/V exhibits a bad behaviour, since the perturbative coefficients grow further with the loop order such that the concept of asymptotic convergence of the perturbative series appears to be meaningful only for rather small couplings or masses far above the mass of the b-quark. At the scale of the b-quark mass, for instance, every known perturbative order approximately contributes by an equal amount. 2 The worrying behaviour of perturbation theory for C PS/V may also be read off from the data points in figure 4, where in fact no signs of an asymptotic convergence with the loop order in C PS/V from about the b-quark mass scale on of any of the curves Figure 5. Comparison of the static extrapolations of Y PS (θ 0 , z) and Y V (θ 0 , z) to the nonperturbative HQET results in the continuum (indicated as the data points in the middle bewteen the panels). In the left panels we use conversion functions C X , X = PS, V, evaluated at two-loop, while in the right panels they are evaluated at three-loop order of perturbation theory. Additionally, the lower panels show (weighted) quadratic fits to all data points, while the upper panels only include data points in the expected applicability domain of HQET (1/z < 0.2). A linear extrapolation over the further restricted range 1/z 0.13 would lead to compatible results.
through the points can be stated. This is even more so for in the right panel of figure 4, as a direct effective finite-volume estimates of f B /f B * involving C PS/V , too, where Γ PS /Γ V cancels the ratio of meson masses in Y PS/V . The growing deviations of the data points from the polynomial fits in 1/z towards the charm quark scale for this quantity are inherited from the corresponding behaviour of the higher-order terms in Γ PS /Γ V (see previous paragraph).

Asymptotic behaviour of effective decay constants
We now turn to an example, where the serious concerns about the usefulness of perturbation theory for the evaluation of the conversion functions raised in the foregoing discussion yet becomes evident in a mismatch between the large-mass asymptotics on the QCD side and a non-trivial, non-perturbative HQET prediction itself. These are the effective finitevolume pseudoscalar and vector meson decay constants Y PS and Y V , which according to static limit in QCD static HQET static limit in QCD  (29) 1.299 (14) 1.2400 (41) 1.2630(45) Table 1. Selected results of static extrapolations using the three-loop C X , X = PS, V, and the associated non-perturbative result computed in static effective theory; all numbers refer to the continuum limit.
subsection 2.5 have to obey the predictions in the asymptotic regime of 1/z → 0. Recall that X RGI is the renormalization group invariant matrix element of the static axial current, eq. (2.17), and its occurrence as the static limit of both QCD observables is a consequence of the degeneracy of pseudoscalar and vector channels at static order of HQET owing to the heavy-quark spin symmetry. As outlined around (2.17) and in appendix C, the renormalization factors entering in X RGI were determined non-perturbatively in the Schrödinger functional renormalization scheme [27] so that X RGI is numerically available without perturbative uncertainties at an overall precision of about 1%, see table 1.
The comparison between the z-dependence of the small-volume pseudoscalar and vector meson decay constants, together with its static extrapolations, and the non-perturbative HQET results, after prior continuum limit extrapolations of all individual pieces involved, is presented in figure 5. The various panels distinguish between two-and three-loop perturbative evaluations of the conversion functions C PS , C V , respectively, as well as between extrapolations quadratic in 1/z including data with 1/z < 0.2 only and over the whole range. We also remark that by linear fits over the further restricted range 1/z 0.13 we arrive at compatible extrapolations. All errors from our numerical simulations and extrapolations were taken into account as explained earlier, but we obviously can not do so for any systematic error from the conversion functions because of their perturbative nature. While the degeneracy of pseudoscalar and vector channels at static order is nicely reproduced by the unconstrained fits via the coincidence of the respective 1/z = 0 limits of Y PS /C PS and Y V /C V , the agreement of the extrapolation of the relativistic QCD results with the associated predictions at static order of HQET does not look very convincing. As can be inferred from table 1 and figure 5, the results obtained in the static effective theory (black data points in the center of the figure) differ systematically for each value of θ 0 from the results of the static extrapolations using unconstrained quadratic fits that represent the data very well. Although these differences tend to decrease when going from the two-loop to the three-loop evaluation of the C X , X = PS, V, the disagreement still remains at the 1−2σ level of the statistical errors. One thus may speculate whether this is just an unfortunate statistical effect, but given the previously discussed doubts on the reliability of perturbation when matching the quark-mass dependent QCD results to HQET via the conversion functions, it may also very well be attributed to the only perturbative approximation of C PS and C V .
In order to provide further support for the conjecture that the use of perturbation theory for the C X is actually responsible for the observed disagreement between the results of the 1/z → 0 extrapolations of the QCD data and the genuinely non-perturbative static HQET predictions, let us go back to the very definition of X RGI in eq. (2.17) and appendix C. It is an example for a renormalization group invariant, which is independent of schemes and scales, allowing a clean factorization of observables into a non-perturbative matrix element of some composite field operator and a multiplicative matching (resp. conversion) function that possesses a perturbative expansion. In the situation at hand, we always refer to the axial vector current A 0 in the lowest-order (i.e., static) effective theory, where its multiplicative renormalization factor is not protected against a scale dependence by a suitable axial Ward identity, as it is the case for the axial current in QCD. Along the lines of ref. [29] one can express the lowest-order HQET approximation of the QCD matrix element of A 0 , in slight adaption of our notation introduced before, as such that at leading order in the inverse heavy quark mass, 1/m h , the conversion function C PS defines a RGI-mass scaling function that contains the full (logarithmic) mass dependence, whereas the non-perturbative matrix element in the static effective theory, X RGI , becomes a pure mass independent number. Here, the renormalization scale µ of the static current is identified with µ = m * , where m * is implicitly defined by the solution of m * = m(m * ), m being the renormalized (running) heavy quark mass. The mass dependence is induced via the renormalized coupling g * ≡ḡ(m * ) that is understood as a function of the ratio of renormalization group invariants, g * = g * (M/Λ), and can be determined in perturbation theory for any value of M/Λ from an integral expression for Λ/M , analogous to eq. (3.8), but in terms of the beta-function and the quark mass anomalous dimension. Moreover, γ match in (3.8) denotes the anomalous dimension in a renormalization scheme for the static axial current, A stat 0 , called the matching scheme, which is defined by the condition that -in this scheme and at µ = m * -matrix elements of A stat 0 are equal to the QCD ones up to O(1/m h ) (cf. appendix B, and refs. [15,29,32] for more details). 3 The crucial observation at this point is now that, in the transition to renormalization group invariants leading to (3.8), the perturbative running enters in C PS (M/Λ) = Y PS (m * )/X RGI through the only perturbative knowledge of the beta-function and the anomalous dimensions of the quark mass and A stat 0 . Hence, it does not come as a too big surprise that the large-mass extrapolation of Y PS (m * )/C PS (M/Λ), which combines fully non-perturbative QCD results with a perturbative matching function, fails to meet the expected static order HQET limit X RGI = X R (µ = L −1 1 )/ X R (µ = L −1 1 )/X RGI , a fully non-perturbative number, both factors of which are precisely known for our setup via the data from the non-perturbative computation in [27], see eq. (2.17) and appendix C. On the other hand, replacing in the calculation of X RGI the non-perturbative value X RGI /X R (µ = L −1 1 ) = Z stat A,RGI /Z stat A (µ = L −1 1 ) = 0.875(7) of (2.17) by the corresponding value at µ = L −1 1 obtained from perturbative running (using the available three-loop beta-function and two-loop anomalous dimension of A stat 0 in the SF renormalization scheme), the black HQET data points in the center of figure 5 receive a downward shift of about 5% to finally coincide with the polynomial extrapolation results of Y X /C X , X = PS, V, in the static limit. This finding clearly suggests that the quark mass dependence of the conversion function C PS , which employs inputs from the perturbative matching between HQET and QCD only, is then also only enough to reproduce the perturbative prediction for the matrix elements in the static effective theory. However, it is not able to reproduce the non-perturbative result, since it does not comprise the full non-perturbative mass dependence that is required for a fully consistent matching between HQET and QCD.
All in all we therefore conclude that with only perturbative knowledge of the conversion functions one is not automatically guaranteed to extrapolate relativistic QCD data at finite values of the (heavy) quark mass to the correct static limit of HQET as z → ∞ (resp. M → ∞). In fact, C PS (and C V ) at three-loop accuracy do not seem to be qualified for use in conjunction with non-perturbative QCD results on the small-volume meson decay constants to extrapolate their heavy quark mass dependence to the genuinely non-perturbative prediction of the effective theory. The distinct mismatch between the QCD extrapolations and the expected HQET results in figure 5 clearly illustrates this. A possible reason for the use of three-loop conversion functions such as C PS not to correctly recover the expected static limit can be traced back to the not that well behaved perturbative expansion of the underlying anomalous dimension γ match in the aforementioned matching scheme, which enters in C PS in addition to the well behaved perturbative series of the beta-function and the quark mass anomalous dimension. Although the overall mass dependence of, say, C PS in figure 10 of appendix B evaluated for different perturbative orders looks quite innocent, a more careful estimation of the coefficients of γ match for the different known orders (up to three loops) still gives rise to worries about neglecting higher-order terms at values of the renormalized coupling around the b-quark mass scale or even below [29].
Our observations on the small-volume decay constants should also be taken as a warning that the method of extracting, e.g., the B-meson decay constant via an interpolation of large-volume lattice data in the charm region and HQET data in the static limit (as sometimes adopted when applying lattice QCD to B-physics phenomenology) can easily yield misleading results, as long as only perturbation theory is employed in the matching step of relating the HQET numbers to the quark mass-dependent ones in QCD.
In general, of course, these statements on the influence of the conversion functions on the validity of HQET as an effective theory of QCD may depend on the individual observable in question and should rather be investigated on a case-by-case basis as we do it in the present study. For instance, in the discussion of the meson masses and the heavy-light . Asymptotic behaviour of the spin splitting observable R spin (z, θ 0 ) together with a quadratic unconstrained (and weighted) fit including data at 1/z ≤ 1/4. In the right panels we compare the static extrapolation in the asymptotic scaling region for a linear (n = 1), quadratic (n = 2) and cubic (n = 3) polynomial fit ansatz.
current ratios above we have already seen that their asymptotic 1/z → 0 behaviour together with the perturbative conversion functions meets the corresponding HQET predictions very well. Even more convincing tests of the effective theory will follow shortly in subsection 3.3, when we come to consider observables in which perturbative factors such as C PS , C V drop out completely.

Asymptotic behaviour of the spin splitting
As can be inferred from figure 6, our spin splitting observable R spin (z, θ 0 ) shows the expected asymptotic behaviour towards the 1/z → 0 limit, where it has to vanish due to the heavy-quark spin symmetry. Opposed to the case of the static axial current, data for the corresponding RGI matrix element, X spin RGI , are not available for the two-flavour theory. Hence, we refrain from studying the static limit of our data in the form of applying eq. (3.4) via (2.16) to R spin /C spin , as we did for the axial and vector meson decay constants in the previous paragraph. Rather, we studied different static extrapolations using a linear (n = 1, z ≥ 13), quadratic (n = 2, z ≥ 4) and cubic (n = 3, z ≥ 3) fit ansatz for the static extrapolation. They are presented for better visibility in the asymptotic region only (right panel of figure 6). The left panel shows the n = 2 case with all available data points; the fit ansaetze are able to describe the data very well, and its behaviour confirms the HQET expectation.

Results without perturbative conversion functions
We now turn our attention to quantities that do not depend on any conversion functions and as such are free of any influence of perturbative uncertainties; they thus are expected to exhibit an unambiguous static extrapolation compatible with their heavy quark mass   Figure 7. Extrapolations of Y PS/PS (z, θ 1 , θ 2 ) and Y V/V (z, θ 1 , θ 2 ) to the static limit for all three combinations of (θ 1 , θ 2 ). The results are listed in table 2. The right panel shows a scaled excerpt of the heavy quark mass region. Black data points indicate the continuum results of the corresponding quantity at static order of HQET. Its continuum extrapolation is displayed in figure 2. For comparison, the lower panel displays an extrapolation with just a linear function in 1/z for 1/z ≤ 0.2 that leads to an equally well confirmation of the HQET expectation.
static limit in QCD static HQET static limit in QCD static HQET  expansion. In most cases we find that a fit function such as  models their z-dependence in the continuum limit very well over the whole region of available heavy quark mass values.

Ratios of currents
As prototype observables, we first consider and probe eq. (2.18c), viz.
where -again owing to the heavy-quark spin symmetry -the static extrapolation of the ratio of effective vector current matrix elements computed at different kinematical parameters, Y V/V (z, θ 1 , θ 2 ), is expected to agree in the limit 1/z → 0 with the associated ratio in the pseudoscalar channel, Y PS/PS (z, θ 1 , θ 2 ), and to approach the common, nontrivial leading-order HQET prediction R stat X (θ 1 , θ 2 ), which denotes the corresponding ratio of matrix elements computed by replacing the relativistic fields by the static ones.
Performing an unconstrained extrapolating fit of all data points (z ≥ 2), according to the fit ansatz (3.9), gives an asymptotic behaviour as depicted in figure 7 for the three available θ-combinations. The black circles (slightly moved to the left towards negative 1/z for ease of presentation) represent the results for the ratio of static-order HQET matrix elements, R stat X (θ 1 , θ 2 ). As usual throughout this work, all data points were extrapolated to the continuum limit first; in particular, the continuum extrapolation of R stat X (θ 1 , θ 2 ) is presented in the left panel of figure 2 and its numerical values can be read off together with the results of the static extrapolations from table 2.
In contrast to the decay constant discussed in the previous subsection, the comparison in figure 7 illustrates a very good agreement between the static extrapolations of the results on Y V/V , Y PS/PS and the numbers for R stat X , computed to even much higher statistical accuracy directly in the static approximation. Moreover, from the fits in the lower panel of the figure one reads off an equally excellent agreement with the HQET expectation from an extrapolation with just a linear function in 1/z for z ≥ 5 (1/z ≤ 0.2). These findings not only demonstrate the correctness of the effective theory itself, but also strongly advocate such observables, which are not affected by any perturbative imperfections deriving from conversion functions, as most promising candidates for observables for use in a non-perturbative strategy for matching HQET to QCD in the spirit of [5]. Two additional observations worth to be mentioned are: a) the error for all three combinations of (θ 1 , θ 2 ) grows with increasing θ 2 2 − θ 2 1 in a similar way for all three observables; b) also the 1/zterms, i.e., the slope at 1/z = 0, grows with that difference. These features also hold true for the quantities considered next and in principle can serve as further helpful criteria for a sensible choice of matching observables.

Ratios of boundary-to-boundary matrix elements
At next, we look at the static extrapolation of the quantities entering the prediction (2.18a). Compared to observables studied before, the QCD data points obtained in the vector (R k ) and pseudoscalar channel (R f ) lie quite close to each other already at finite quark mass. From figure 8 (left panel) one concludes that superficially the quadratic fit ansatz (3.9) very well represents the data points, which approach a (due to spin symmetry) common HQET limit for the three θ-combinations. Note that the static HQET results for R k and R f , independently computed in the continuum limit from different simulations, have not been constrained to be equal, but their agreement is just excellent though. However, only the extrapolation for the combination (θ 1 , θ 2 ) = (0, 0.5) leads to the correct result in the static limit, R k = R f , represented in the figure by the leftmost black data points. The results for this particular static extrapolation and the corresponding HQET results are given in table 2.
This partial mismatch of the results from extrapolations of the QCD data over the whole available z-range down to z = 2 (1/z = 0.5) with the HQET predictions should be taken as a warning that the validity of the HQET-inspired 1/m h -expansion of heavy-light QCD obeservables in the large-mass regime can not be tacitly trusted down to the charm quark scale or below. In fact, if we restrict the fits to include data points for 1/z < 0.26 only, as done in the previous section, we obtain a static extrapolation as shown in the right panel of figure 8, where the static numbers can be seen to be covered by the error of the extrapolation results. An even further restriction of the fit interval to 1/z < 0.1 would easily allow for a feasible linear interpolating fit including the HQET result as a constraint. Here, extending those fits to the lower z's that were excluded yields only mild deviations of about 1σ between the fit functions and the data points, though, but this may also be more pronounced for other observables. Turning the argument around, an interpolating fit somewhat arbitrarily extended, at fixed polynomial degree, to include data below the charm region (z = 2), which appears to give a good description of the data, can lead to an underestimation of the error of the static extrapolation and thereby pretend to miss the HQET prediction.  1/z Figure 9. Static extrapolations of R PS/PS (z, θ 1 , θ 2 ), R V/V (z, θ 1 , θ 2 ) and R P/P (z, θ 1 , θ 2 ) for all three combinations of (θ 1 , θ 2 ) in comparison to their HQET counterparts after taking the continuum limit in the static effective theory (these data points are slightly shifted for better readability). The results are listed in table 2. The right panel shows a scaled excerpt of the heavy quark mass region individually for the three available (θ 1 , θ 2 )-combinations such that error bars of the HQET results become visible.
static limit in QCD static HQET  Table 3. Results for static extrapolations of R PS/PS , R V/V and R P/P in continuum QCD as displayed in figure 9 and its corresponding non-perturbative continuum extrapolation results in the static effective theory.

Ratios of boundary-to-bulk matrix elements
Finally, we check the static extrapolation of eq. (2.18b). The HQET results for R stat PS (θ 1 , θ 2 ) follow from the continuum extrapolation presented in figure 2. Their values are given in table 3, together with the results that stem from a static extrapolation of the (continuum) QCD observables R PS/PS , R V/V and R P/P as defined in eq. (2.9). Once more, the data sets themselves are very well represented by a quadratic polynomial fit ansatz over the whole range of available data points, whereas the static HQET result (common to all three QCD observables, again due to spin symmetry) is only met for the θ-combination (θ 1 , θ 2 ) = (0, 0.5) as in the foregoing case of ratios of boundary-to-boundary matrix elements. Accordingly, the same discussion (and warning) carries over literally here, except for R P/P that very well extrapolates to the HQET result for all θ-combinations studied. Since the errors on the data points for R P/P stay roughly the same when going from the scale of the b-quark mass down to z = 2, the slight mismatch in the other cases may likely be attributed partly to statistical effects.

Conclusions
We have studied the asymptotic large-mass behaviour of heavy-light meson observables in lattice QCD with two massless dynamical sea quarks, in a small volume of linear extent L 1 ≈ 0.4 fm and for heavy quark mass values within a range from beyond the b-to below the c-quark scale, in order to confront them with their HQET predictions.
Having taken the continuum limit in all parts of our entirely non-perturbative calculations on both the QCD and the HQET side and subsequently performed unconstrained static extrapolations along the limit z = L 1 M h → ∞, in most cases we generically find (within the numerical precision on our results with errors at the few-% level) a very satisfactory -sometimes an even excellent -agreement between the large-mass asymptotics of the QCD observables and their expected leading-order HQET limits.
Moreover, the overall quality of the polynomial fits in 1/z to the heavy quark mass dependence of the (continuum) heavy-light QCD observables convincingly demonstrates that the theory is very well described by simple 1/m h -corrections to the static limit of the effective theory. In particular, for 1/z 0.1 our extrapolating fits to the HQET predictions can consistently be modeled by functions linear in 1/z. We are thus led to conclude that the effective theory is very well tested and that the regime with 1/z 0.1, which is a key ingredient to the finite-volume matching part of the ALPHA Collaboration's B-physics programme based on HQET non-perturbatively renormalized and matched to QCD at O(1/m h ) [5][6][7][8][9][10][11][12], lies very well within the applicability domain of HQET. This is, for instance, also in line with the onset of the linear behaviour reported in the tree-level study [13] of the 1/z-dependence of the HQET parameters contributing to the full set of heavy-light flavour currents in HQET at O(1/m h ). Alltogether, these findings are very reassuring, since they imply that in the finite-volume setup of the aforementioned nonperturbative matching strategy between HQET and QCD, higher-order corrections beyond the O(1/m h ) ones already included can be expected to be suppressed by a factor of about 10.
A prominent exception to this favourable outcome of our tests of HQET consists in the large-mass asymptotics of the small-volume heavy-light axial vector and vector meson decay constants, which fails to meet the leading-order HQET prediction in the static limit. But, as argued in section 3.2, rather than interpreting this as an inherent shortcoming of HQET being not an appropriate and predictive effective theory of QCD, one should take this as an advice that it is generally safer not to use conversion functions such as C PS , C V (which inevitably enter in a consistent comparison of QCD and HQET decay constants) from perturbation theory, even if they are determined at three-loop accuracy, i.e., with relative errors of orderḡ 6 (m * ) at some intrinsic mass scale m * ; e.g., their perturbative convergence appears to be relatively poor still at the b-quark scale. Otherwise, the combination of nonperturbative QCD data with the mass dependence via perturbative matching functions (encoding the running in HQET) to recover the correct HQET limit can easily lead to inconsistencies between QCD and non-perturbatively computed matrix elements in the effective theory. For the decay constant, we encounter a systematic effect of up to ∼ 5% from relying on perturbative running in the effective theory. Therefore, we consider this as a warning when, e.g., the physical (i.e., large-volume) B-meson decay constant is being extracted involving interpolations between QCD data below the b-scale and the static limit, and advocate to perform the matching entirely non-perturbatively.
As already noted earlier in [29], a more quantitative understanding of this deficiency can be gained from the relative error that results from a truncation of the perturbative matching expression at l-loop order, viz.
I.e., since this perturbative uncertainty only decreases logarithmically as m * becomes large, at some point it starts to dominate over the power correction that one needs to include at next-to-leading order in the HQET expansion for precision physics at the b-quark mass scale. This underlines once more that with an only perturbative conversion function, a consistent next-to-leading order expansion with errors decreasing as 1/m 2 h can not be achieved. Another interesting aspect of our work is that the continuum extrapolations of the HQET observables, such as those presented in figure 2 of section 3.1 for two discretizations of the static quark, provide strong numerical evidence that an universal continuum limit of the static effective theory exists and that the non-perturbative renormalizability of HQET along the finite-volume matching strategy of ref. [5] can be established indeed.
Finally, let us emphasize that exploring the size of the higher-order corrections in our QCD observables to test HQET non-perturbatively may also readily serve as a guide so single out preferred choices among observables, suitable for a specific HQET-QCD matching problem in question, that have only small O(1/m 2 h ) contributions. In addition to that, any flexibility in having different matching equations made of different observables to determine the same (set of) HQET parameters enables further useful checks in actual computations, because the final results should be independent of any specific but sensible choice of matching equations and kinematical parameters (such as T /L, x 0 and the θ's) entering them, up to small O(1/m 2 h ) corrections.
Forschungsgemeinschaft. N. G. is funded by the Leverhulme Trust, research grant RPG-2014-118. Finally, P. F. thanks the organizers and participants of the workshop "Highprecision QCD at low energy" at the "Centro de Ciencias de Benasque Pedro Pascual" for stimulating discussions and comments.

A Definitions
Here we provide the definitions of the observables of this study in terms of the traditional notation for Schrödinger functional (SF) correlation functions. We do not repeat all details here; for explicit expression for SF heavy-light meson correlators in lattice QCD as well as in the static limit of lattice HQET, the reader may, e.g., consult [5,13,33,34].

QCD observables
For the effective masses we use SF correlation functions f X with bulk insertions of O(a) improved heavy-light currents (cf. section 2) in the pseudoscalar and vector channel and define Effective decay constants associated with these correlators as well as suitable ratios thereof require renormalization and are given by The renormalization constants Z X , X = A, V, P, are known non-perturbatively in twoflavour QCD; as in our earlier work [19], Z A and Z V have been taken from [25], while Z P,RGI is available through the scale dependent renormalization factor Z P computed in refs. [9,35]. The improvement coefficients b X (multiplying the bare subtracted heavy quark mass) are known to one-loop order of perturbation theory and can be found in [24]. The spin splitting observable, however, is constructed from a ratio of SF boundaryto-boundary correlators with pseudoscalar and vector channel composite fields, free of any improvement coefficient and renormalization factor: As already explained in the main text, the foregoing finite-volume observables have been computed for one fermionic phase angle out of θ 0 ∈ {0, 0.5, 1}. To enlarge the variety in probing QCD and HQET at different kinematics, it remains to specify the observables that depend on two such angles. Here we build them from SF correlation functions with different θ's, i.e., where in practice we have chosen to extract them from our simulations for the non-trivial combinations (θ 1 , θ 2 ) ∈ {(0, 0.5), (0.5, 1), (0, 1)}. Again, renormalization factors drop out in these ratios.

HQET observables
Among the QCD observables defined above depending on a single θ-angle, only the effective decay constants in eq. (A.2) possess a (common) non-trivial static limit, which is given by up to renormalization factors (cf. subsection 2.5) to be discussed in appendix C, while the ratios in (A.3) approach 1 as a consequence of the heavy-quark spin symmetry that entails pseudoscalar and vector channels to coincide. The ratios in eq. (A.5) depending two θangles, on the other hand, approach in the static limit one of the following observables at static order of HQET: The corresponding static-light correlation functions f stat A and f stat 1 entering in these expressions have first been defined in [34].

B Conversion functions
In this section we summarize the expressions of the perturbative conversion functions C X , X ∈ {PS, V, spin, PS/P, PS/V, mass}, that have been employed to compare our nonperturbatively renormalized observables in QCD to their counterparts in HQET. They are computed and parameterized in the so-called matching scheme, which has been introduced in [32] and specified in appendix B of [15] (see also [29] for another detailed discussion). Our formulae for all relevant operators below refer to the two-flavour theory and are based  Figure 10. HQET-QCD conversion functions C X for heavy-light currents as described in the text. The solid and dashed curves correspond to the three-loop and two-loop order anomalous dimensions in the matching scheme, respectively. In the case of C PS/V (middle panel on the right), the four-loop expression, which is known from [31], is also included as dash-dotted curve. Comparing the curves for different loop orders suggests that perturbation theory converges rather slowly. Vertical dotted lines indicate the fixed values of z = L 1 M used in our study.
on their anomalous dimensions known up to three-loop order in continuum perturbation theory [36][37][38][39][40][41][42][43][44][45][46][47][48], to be combined with the matching coefficients between QCD and the effective theory up to two loops [2,40,[49][50][51], see appendix B of [15]. In order to judge the impact of the order of the perturbative expansion on this comparison between QCD and HQET, introduced by the conversion functions C X as far as they enter the observables under study, we evaluate the C X including the two-loop and three-loop anomalous dimensions (together with the respective matching coefficients in one-and two-loop accuracy) separately. In both cases we use the four-loop beta-function for the coupling [52][53][54][55][56][57], while the conversion function for effective masses involves the quark mass anomalous dimension τ at four-loop [57][58][59].
Moreover, thanks to the three-loop calculation of the matching coefficients for the heavy-light currents available from [31], the anomalous dimensions of ratios of currents become effectively known to four-loop order in the matching scheme, because the unknown four-loop anomalous dimensions of the currents themselves cancel out in the difference of anomalous dimensions contributing to these ratios. As an example, we therefore also include the conversion function C PS/V of the ratio of axial vector (A 0 ) to vector (V k ) currents at four-loop accuracy in our study.
Since so far only C PS for the N f = 2 theory was given earlier in [27], we here list the parametrization of all conversion functions C X entering this study that were obtained along the lines detailed in appendix B of [15]. They are expressed as smooth functions in terms of the variable with the product L 1 Λ MS = 0.629 (36) in the two-flavour theory taken from [20] and M denoting the renormalization group invariant (RGI) heavy quark mass as in the main text. The functions decompose into a pre-factor which encodes the leading logarithmic asymptotics (if any) as x → 0, multiplied by a polynomial in x which reflects the perturbative order of the underlying anomalous dimension in conjunction with the associated matching coefficients.
The functional dependence of the C X on Λ MS /M is shown in figure 10. The solid curves represent the matching functions to (in most cases) highest available perturbative accuracy, i.e., corresponding to three-loop anomalous dimensions of the heavy-light currents involved and the four-loop beta-function. The dashed curves are those obtained with twoloop anomalous dimensions only. For illustration, we plot as vertical dotted lines the values of z which have been fixed in order to non-perturbatively compute our test observables. In the exemplary case of the ratio of pseudoscalar to vector currents, C PS/V , which we are able to consider even up to the maximally available four-loop precision (dash-dotted curve in the middle-right panel of figure 10), one observes an increase in the size of correction when going from two-to three-and three-to four-loop order. This is in accordance with the conclusion of the authors of [31] that the perturbative series for this ratio of currents "converges very slowly at best". We finally remark that, since the previous expressions derive from continuum perturbation theory, the functions C X must be properly attached as factors to the HQET-QCD test observables in question after the lattice results on them have been extrapolated to the continuum limit, cf. section 2.5.

C Further details and tables
Tree-level improvement Even though our simulations are performed at rather fine lattice spacings, one ultimately encounters mass-dependent cutoff effects, which essentially grow with the heavy quark mass, i.e., to some power of z × a/L = aM in the non-perturbatively O improved theory. To attenuate these effects, we also apply perturbative improvement to some of our lattice observables under consideration. Here we restrict ourselves to tree-level improvement (TLI), which for a generic test observable Ω amounts to the replacement Ω(ḡ 2 , a/L, z) → Ω (0) (ḡ 2 , a/L, z) = Ω(ḡ 2 , a/L, z) and thereby removes all the O ( a L ) n effects to produce classically perfect observables. For additional details on the actual extraction of the improvement terms δ (0) Ω , see appendix D of [6]. Besides its obvious dependence on z and a/L, δ (0) Ω also depends on the kinematical setup of the Schrödinger functional, inherited from the observable itself. These are the choice of boundary gauge field, the ratio T /L and the fermionic periodicity angles involved in the definition of Ω. Except for Ω ∈ {R PS/PS , R V/V , R P/P , R 1 , R f , R k }, where TLI has not been accounted for, and observables such as R spin , where it vanishes exactly, the presented results always correspond to the tree-level improved quantities.
We furthermore remark that before any continuum limit in QCD and HQET is taken numerically, we use the data from appendix B of [9] to correct for small deviations in our observables from the renormalized (i.e., constant physics) trajectory at (L 1 m l ,ḡ 2 (L 1 )) = (0, 4.484). The associated errors are propagated quadratically.

Special case: Axial current renormalization
We add some remarks about the continuum extrapolation of the RGI static decay constant as defined in eq. (2.17), because we need to fully specify the renormalization scheme applied.
The total renormalization factor to obtain the RGI matrix element of the static axial current, X RGI , from the bare matrix element, X bare , decomposes into As indicated here, the universal part Z stat A,RGI /Z stat A (µ), which links a renormalized matrix element of the static axial current at a renormalization scale µ to the RGI one, is defined within the massless SF scheme for vanishing boundary field and (θ, ρ = T /L) = (0.5, 1) nonperturbatively and in the continuum limit. For the particular scale µ = L −1 1 corresponding to the physical volume employed in the present work, Z stat A,RGI /Z stat A (µ) has been re-computed from the two-flavour data of ref. [27] to yield the estimate quoted in (2.17). (Note that in [27], the universal ratio was denoted as Φ RGI /Φ(µ)). In the definition of the scale dependent renormalization factor itself, Z stat A (g 0 , aµ), the tree-level normalization factor X(0, L/a) enters at the respective values of (θ, ρ) such that Z stat A approaches one in the limit g 0 → 0, cf. [32]. In fact, as we always have ρ = T /L = 1 in the present work, too, we can take the continuum limit according to with the particular values θ 0 ∈ {0, 0.5, 1} and θ = 0.5. 4 The second equality in (C.3) picks up the notation of eq. (2.17), and the results are listed in table 1. Finally, let us point out that its total error is dominated by the error of the universal continuum factor for the non-perturbative running, Z stat A,RGI Z stat A (µ) (as quoted in (2.17)), and that this part of the error is only to be accounted for after X RGI has been extrapolated to the continuum limit. Tables   In table 4 we collect addition details concerning the measurements of heavy-light QCD observables in the present publication, thereby extending the parameter set of table 3 in [9], also referred to as set "QCD[L 1 ]". Table 5 repeats the bare parameters relevant for the measurements in the static theory, defining ensemble set "HQET[L 1 ]".  Table 5. Bare parameters for the five ensembles used to compute HQET observables in L 1 . The gauge field ensembles have been produced with vanishing SF boundary fields and a kinematical setup specified by (T /L, θ sea ) = (1, 0.5).
Tables 6 -12 list the results at finite lattice spacing and the corresponding continuum limits of some selected observables. Listing all results in full detail would be beyond the scope of this paper 5 , since the qualitative and quantitative behaviour can also be well inferred from the plots shown in the main text. As an example for the θ-dependence of an observable, we reproduce LΓ PS for all values of z ∈ {2, 2.7, 3, 3.3, 4, 6, 7, 9, 11, 13, 15, 18, 21} and θ = θ 0 ∈ {0, 0.5, 1} in table 6. For all other observables we only list values for θ 0 = 0.5 or the combination (θ 1 , θ 2 ) = (0.5, 1). Note that values in squared brackets have not been taken into account in the continuum extrapolations as detailed in section 3.1 and that the phase angle θ ≡ θ sea = 0.5 has been used for the doublet of sea quarks in the production runs to generate the underlying two-flavour gauge field configuration ensembles.
V (x 0 = T /2) and LΓ     Table 9. Tree-level improved ratios F PS/V and R PS/V together with the spin splitting R spin for all available values of z, a/L at θ 0 = 0.5. Continuum limit (CL) has been taken according to eq. (3.1). Table 12. Observables Y PS/PS and Y V/V for all available values of z, a/L at (θ 1 , θ 2 ) = (0.5, 1). Continuum limit (CL) has been taken according to eq. (3.1). Both observables are expected to agree in the static limit.