B-physics from Nf=2 tmQCD: the Standard Model and beyond

We present a lattice QCD computation of the b-quark mass, the B and B_s decay constants, the B-mixing bag parameters for the full four-fermion operator basis as well as determinations for \xi and f_{Bq}\sqrt{B_i^{(q)}} extrapolated to the continuum limit and to the physical pion mass. We used N_f = 2 twisted mass Wilson fermions at four values of the lattice spacing with pion masses ranging from 280 to 500 MeV. Extrapolation in the heavy quark mass from the charm to the bottom quark region has been carried out on ratios of physical quantities computed at nearby quark masses, exploiting the fact that they have an exactly known infinite mass limit. Our results are m_b(m_b, \overline{\rm{MS}})=4.29(12) GeV, f_{Bs}=228(8) MeV, f_{B}=189(8) MeV and f_{Bs}/f_B=1.206(24). Moreover with our results for the bag-parameters we find \xi=1.225(31), B_1^{(s)}/B_1^{(d)}=1.01(2), f_{Bd}\sqrt{\hat{B}_{1}^{(d)}} = 216(10) MeV and f_{Bs}\sqrt{\hat{B}_{1}^{(s)}} = 262(10) MeV. We also computed the bag parameters for the complete basis of the four-fermion operators which are required in beyond the SM theories. By using these results for the bag parameters we are able to provide a refined Unitarity Triangle analysis in the presence of New Physics, improving the bounds coming from B_{(s)}-\bar B_{(s)} mixing.

values of the lattice spacing with pion masses ranging from 280 to 500 MeV. Extrapolation in the heavy quark mass from the charm to the bottom quark region has been carried out on ratios of physical quantities computed at nearby quark masses, exploiting the fact that they have an exactly known infinite mass limit. Our results are

Introduction
Physical processes in the B-sector are crucial to perform accurate tests of the Standard Model (SM) and search for possible signals of New Physics (NP). The experimental accuracy in flavour processes has recently been significantly increased by the B factories, and new measurements are being performed, mainly thanks to the remarkable performance of the dedicated experiment LHCb. On the theoretical side, lattice computations have entered during the past few years a precision era, in which the target per cent precision for some of the relevant hadronic parameters in Flavour Physics is becoming accessible. In particular, in the study of B-physics processes, there has been substantial progress thanks to alternative lattice methods and techniques aimed at treating the heavy quarks on the lattice with controlled systematic uncertainties. For a recent review of lattice results see Ref. [1]. Lattice methods are irreplaceable for the calculation of the so called golden plated processes since hadronisation effects are fully under control. They lead to accurate determinations of decay constants, form factors and bag-parameters. For example, the leptonic decays B → τ ν τ and B 0 (s) → µ + µ − receive precise input information from lattice determinations of the B-mesons decay constants, that are necessary for the experimental results to acquire their full physical interpretation. At present, the world average of the B-meson leptonic decays BR(B → τ ν τ ) = (1.14 ± 0.22)10 −4 [2,3,4], which is potentially sensitive to NP effects already at tree level, turns out to be in agreement with the SM prediction BR(B → τ ν τ ) SM = (0.81 ± 0.07)10 −4 [5,6], and also the recent measurements of the B 0 s → µ + µ − decay [7,8,9,10] have given a first, remarkable evidence of SM consistency [11,12,6]. The neutral B-meson mixings, which can only occur at the loop level in the SM, could be a privileged candidate process for detecting amplified NP effects, and indeed they play a crucial role in the Unitarity Triangle (UT) analysis (for recent results, see for example [13,14,15]). While in the SM the frequency of the oscillations, ∆M B (s) , receives the contribution from a single four-quark operator, the knowledge of the bag-parameters of the full four-fermion operator basis is required to get predictions in general NP extensions of the SM. Two of these B-parameters also enter the SM prediction of the lifetime difference ∆Γ s of neutral B s mesons, which has been recently measured rather precisely by LHCb [16].
In this paper we use gauge configurations with N f = 2 dynamical quarks at four values of the lattice spacing, generated by European Twisted Mass Collaboration (ETMC), to obtain the continuum limit results for a number of physical quantities that are relevant for B-Physics. These are the b-quark mass m b , the pseudoscalar decay constants f B and f Bs , and the bag-parameters of the full basis of ∆B = 2 four-fermion operators. In our previous paper [17], we provided a determination of the b-quark mass and of the B-meson decay constants obtained by studying the heavy quark on the lattice with the so called ratio method, proposed in Ref. [18]. The same strategy is also applied in the present study, and the lattice calculation presented here is based on the same set of gauge configurations used in our previous study. Nevertheless, several new results and improvements are presented in this paper. The main novelties, with respect to Ref. [17], are the following: -We have computed the full set of B-parameters for the ∆B = 2 four-fermion operators, which are relevant for B-meson mixings, within and beyond the SM, and for the theoretical predictions of the neutral B-meson lifetime differences ∆Γ (s) . For the full set of B-parameters, this is the first lattice calculation which takes into account the effect of dynamical quarks (preliminary results obtained with N f = 2 + 1 dynamical quarks have been presented in [19]). We have used these results to provide a refined Unitarity Triangle analysis improving the bounds coming from B (s) -meson mixing constraints on NP.
-We have computed 2-and 3-point correlation functions by employing optimized smearing techniques. Given the temporal extensions of our lattices, while the use of smearing interpolating operators is mandatory to extract a signal for the B-parameters from the 3-point correlation functions, the suppression of the excited states contribution helps in improving also the determination of the b-quark mass and the decay constants. For this reason, the results obtained in this paper for m b , f B and f Bs should be considered as an update of those given in Ref. [17].
-One of the main sources of uncertainty in the determination of the decay constant f B in Ref. [17] was due to the chiral extrapolation. In this study, we reduce this uncertainty by making use of the observation that the double ratio of decay constants (f Bs /f B )/(f K /f π ) exhibits a much smoother chiral behavior with respect to f Bs /f B itself, due to a large numerical cancellation of the corresponding chiral logarithms [20]. Therefore, also in this respect, the present determinations of f Bs /f B and of f B represent an improvement over the results of Ref. [17].
The plan of this paper is as follows. In Section 2, based on the results of this work for the ∆B = 2 bag parameters, we discuss the implications for NP of our updated Unitarity Triangle analysis. In Section 3 we give information about lattice simulation details and we describe the techniques that have been used in this work. In Sections 4 and 5 we present and apply our strategy, namely the ratio method, in order to get continuum limit determinations for the b-quark mass, pseudoscalar decay constants of the B and B s mesons, and the bag parameters for the full four-fermion operator basis, as well as other interesting quantities like ξ and f Bq B (q) i (i = 1, . . . , 5 and q = d/s). In Section 6 we summarise the final results and discuss our error budget. We also provide some comparison plots between our numbers and those obtained by other lattice collaborations.
For reader's convenience we immediately give here our final results. For each quantity the quoted error corresponds to the total uncertainty which is the sum in quadrature of the statistical and systematic error.
The result for the b-quark mass is given in the MS scheme at the scale of its own value, m b . We perform the running up to m b using either N f = 2 or N f = 4, we take the average over the two results and we consider their half difference as a systematic uncertainty 1  As a by-product of our work we have computed the decay constants for the D s and D mesons as well as their ratio. They read: f Ds = 250(7) MeV, f D = 208(7) MeV, f Ds /f D = 1.201 (21). (1.5) The most general form of the ∆F = 2 effective weak Hamiltonian is where in the so-called SUSY basis ( [21], [22]) the four-fermion operators O i andÕ i read For the neutral B-meson system h ≡ b and q ≡ d or s with α and β denoting color indices. Spin indices are implicitly contracted within square brackets. The Wilson coefficients C i andC i have an implicit renormalization scale dependence which is compensated by the scale dependence of the renormalization constants of the corresponding operators. Notice that the parity-even parts of the operators O i andÕ i are identical. Due to parity conservation in strong interactions, for the study of B 0 q − B 0 q oscillations it is sufficient to consider only the matrix elements B 0 q |O i |B 0 q , where by O i (i = 1, . . . , 5) we denote the parity-even components of the operators (1.7). We recall that in the SM only the matrix element of the operator O 1 is relevant.
The bag parameters, B i (i = 1, . . . , 5), provide the value of four-fermion matrix elements in units of the magnitude of their vacuum saturation approximation. More explicitly, they are defined by the equations [23,24] with In Table 1 we collect our results for the bag parameters (see Eqs (1.9) and and (1.10)) of the full operator basis (i.e. the parity even componenents of the operators in Eq. (1.7)) in the RI/MOM scheme at the scale of the b-quark mass (Eq. (1.1)). In Table 2 we gather results for the bag parameters expressed in the MS scheme of Ref. [25] at the scale of the b-quark mass. Also, in Table 3 5). Bag parameters are expressed in the MS scheme of Ref. [25] at the scale of the b-quark mass.
Finally, in Table 4 we collect our results for the quantities f Bq B (q) i where q = d, s and i = 1, . . . , 5. Again the bag parameters are expressed in the MS scheme of Ref. [25] at the scale of the b-quark mass. For convenience we also give our results for the SM relevant quantities in which the bag parameters are expressed in the RGI scheme and we have employed in the running N f = 5 and Λ [27]. We get The RGI values of the bag parameters corresponding to the SM four-fermion operators read 2 has been obtained with N f = 2 + 1 dynamical quarks in [28], while preliminary results for the full basis with N f = 2 + 1 have been presented in [19]. ∆F = 2 processes provide some of the most stringent constraints on NP generalizations of the SM. Several phenomenological analyses of ∆F = 2 processes have been performed in the last years, both for specific models and in model-independent frameworks [29,30,31,32,33,34,35,36,37,38,39]. A generalization of the Unitarity Triangle (UT) analysis, which allows for NP effects by including the most significant flavor constraints on NP available at the time was performed in Ref. [29]. The result was a simultaneous determination of the CKM parameters and the size of NP contributions to ∆F = 2 processes in the neutral kaon and B (s) meson sectors. The NP generalization of the UT analysis consists in including in the theoretical parametrization of the various observables the matrix elements of operators which, though absent in the SM, may appear in some of its extensions.
In a previous paper [39] we have presented the first (N f = 2) unquenched, continuum limit, lattice QCD results for the matrix elements of the operators describing neutral kaon oscillations in extensions of the SM. In the same paper we have updated the UT analysis allowing for possible NP effects, improving the bounds coming from K 0 −K 0 mixing constraints.
In a similar way, we present here the N f = 2 lattice QCD results for the bag parameters of the full basis of ∆B = 2 four-fermion operators and we use them in updating the UT analysis beyond the SM. The new ingredients entering the analysis are collected in Table 2. For all the other input data we use the numbers quoted in Ref. [6] in the Winter 2013 analysis.
In the NP-oriented analysis, the relations among experimental observables and the CKM matrix elements are extended by taking into consideration the most general form of the ∆F = 2 effective weak Hamiltonian, given in Eq. (1.6). In the present analysis we focus on ∆B = 2 processes. The effective weak Hamiltonian is parameterized by Wilson coefficients of the form where F i is the (generally complex) relevant NP flavor coupling, L i is a (loop) factor which depends on the interactions that generate C i (Λ), and Λ is the scale of NP, i.e. the typical mass of new particles mediating ∆B = 2 transitions. For a generic strongly interacting theory with an unconstrained flavor structure, one expects F i ∼ L i ∼ 1, so that the phenomenologically allowed range for each of the Wilson coefficients can be immediately translated into a lower bound on Λ. Specific assumptions on the flavor structure of NP correspond to special choices of the F i functions. Following Ref. [29], in deriving the lower bounds on the NP scale Λ, we assume L i = 1, that corresponds to strongly-interacting and/or tree-level coupled NP. Two other interesting possibilities are given by loop-mediated NP contributions proportional to either α 2 s or α 2 W . The first case corresponds for example to gluino exchange in the minimal supersymmetric SM. The second case applies to all models with SM-like loopmediated weak interactions. To obtain the lower bound on Λ entailed by loop-mediated contributions, one simply has to multiply the bounds we quote in the following by α s (Λ) ∼ 0.1 or α W ∼ 0.03.
The results for the upper bounds on the |C B d i | and |C Bs i | coefficients and the corresponding lower bounds on the NP scale Λ are collected in Tables 5 and 6, where they are compared to the previous results of Ref. [29]. The superscript B d or B s is to recall that we are reporting the bounds coming from the B d -and B s -meson sectors we are here analyzing. The constraints on the Wilson coefficients of the non-standard operators and, consequently, on the NP scale turn out to be significantly more stringent than in Ref. [29], in particular for the B s sector. Both experimental and theoretical inputs have been updated with respect to Ref. [29] (see Ref. [6]). We notice, in particular, that the input values used in Ref. [29] for B (d/s) i were obtained in Ref. [24] in the quenched approximation, at rather large pion masses and at only one lattice spacing (a ∼ 0.1 fm).
We observe that the analysis is performed (as in [29]) by switching on one coefficient at the time in each sector, thus excluding the possibility of having accidental cancellations among the contributions of different operators. Therefore, the reader should keep in mind that the bounds may be weakened if, instead, some accidental cancellation occurs.
In Figs. 1 and 2 we show the comparison between the lower bounds on the NP scale obtained for the case of a generic strongly interacting NP with generic flavor structure by the constraints on the |C B d i | and |C Bs i | coefficients coming from the present generalized UT analysis, and the previous results of Ref. [29].
Comparing with the results of the UT-analysis in Ref [39], we notice that (at least for generic NP models with unconstrained flavour structure) the bounds on the NP scale coming from K 0 -K 0 matrix elements turn out to be the most stringent ones.  Table 6: 95% upper bounds for the |C Bs i | coefficients and the corresponding lower bounds on the NP scale, Λ, for a generic strongly interacting NP with generic flavor structure (L i = F i = 1). In the lower panel the results of [29] are displayed for comparison.  The lower bounds on the NP scale, provided by the constraints on |C Bs i | (i = 1, . . . , 5) for generic NP flavor structure, are shown as brown bars. For comparison, we plot the bounds of Ref. [29] as yellow bars.

Lattice setup and simulation details
The ETM Collaboration has generated N f = 2 gauge configuration ensembles at four values of the inverse bare gauge coupling, β, and at a number of light sea quark masses. The values of the simulated lattice spacings lie in the interval [0.05, 0.1] fm. Dynamical quark simulations have been performed using the tree-level improved Symanzik gauge action [40] and the Wilson twisted mass action [41] tuned to maximal twist [42]. Bare quark mass parameters, corresponding to a degenerate bare mass value of the u/d quark, are chosen so as to have the light pseudoscalar mesons ("pions") in the range 280 ≤ m PS ≤ 500 MeV. A list of the simulated charged pseudoscalar meson masses is given in [48]. Discussion about the computation of the neutral pseudoscalar meson mass using twisted mass fermions has been presented in [44,45]. More details on the action and our N f = 2 gauge ensembles can be found in Refs. [44,45,43]. We stress that the use of maximally twisted fermionic action offers the advantage of automatic O(a) improvement for all the interesting physical observables computed on the lattice [42].
In the present work we treat the strange and the charm quarks as quenched. We have computed 2-and 3-point correlation functions using valence quark masses whose range is extended from the light sea quark mass up to 2.5-3 times the charm quark mass. Simulation details are given in Table 7, where µ ℓ , µ s and µ h indicate the bare light, strange-like and heavy (i.e. charm-like and heavier) valence quark masses respectively.  We have set light valence quark mass values equal to the light sea ones, aµ ℓ = aµ sea . Renormalised quark masses, µ R , are obtained by the bare ones using the renormalisation constant (RC) Z µ = Z −1 P , µ R = µ/Z P [41,46]. The values for Z P at the three coarsest lattice spacings have been computed in [47] using RI-MOM techniques. Following the same method we have also computed Z P at the finest value of the lattice spacing corresponding to β = 4.20. All Z P expressed in the MS scheme at 2 GeV are gathered in Appendix C of [39]. Here we use the corresponding Z P values at the scale of 3 GeV, which is in the region of momenta directly accessible in the RI-MOM calculation of Ref. [47]. In Table 8 we collect the values of Z P at each value of β as well as the corresponding lattice spacing values. The latter have been computed employing SU(2) (NLO) ChPT formulae to fit in a combined way our data for the pion mass and decay constant by using as an input the experimental value of the pion decay constant [48]. Moreover in [48], we have computed the values We have computed 2-and 3-point correlation functions by employing smearing techniques on a set of 100-240 independent gauge configurations for each ensemble and evaluated statistical errors using the bootstrap method. Smeared interpolating operators become mandatory in the cases where relativistic heavy quarks are involved. Smearing proves to be beneficial in reducing the coupling of the interpolating field with the excited states, thus increasing its projection onto the lowest energy eigenstate. The usual drawback, i.e. increase of the gauge noise due to fluctuations of the links entering in the smeared fields, is controlled by replacing thin gauge links with APE smeared ones [49]. With this technical improvement we can extract heavy-light meson masses and matrix elements at relatively small temporal separations while keeping noise-to-signal ratio under control. We employed Gaussian smearing [50,51] for heavy-light meson interpolating fields at the source and/or the sink. The smeared field is of the form: where Φ L is a standard local source and ∇ APE is the lattice covariant derivative with APE smeared gauge links characterised by the parameters α APE = 0.5 and N APE = 20. We have taken κ G = 4 and N G = 30. We have noticed that in practice we get better overlap with the ground state when the source, rather than the sink, is smeared. Thus 2-point Smeared-Local (SL) correlation functions yield more improved plateaux for the lowest energy mass state than Local-Smeared (LS) or Smeared-Smeared (SS) ones. Even stronger overlap with the ground state is achieved with the use of an optimised source constructed as follows: where we have introduced the tunable parameter w. In practice, the use of these sources does not involve other inversions than those of the local and smeared sources. We constructed correlators that have Φ opt as a source and either Φ L or Φ S in the sink. We verified that in general the optimised correlators fulfil the expectations of providing an earlier Euclidean time projection on the ground state than the (SL) correlators. In Fig. 3(a) we show an example of an improved ground state plateau using smeared source for the quark masses (aµ ℓ , aµ h ) = (0.0080, 0.5246) at β = 3.80. Fig. 3(b) shows the improvement that we achieve when we use the optimised source of Eq. (3.2) with an appropriately tuned w parameter. In the figure the effect of tuning w at the level of the first decimal place is illustrated 3 . In our application however, we tried an even better tuning, e.g. up to the second digit. Our general conclusion is that employing the optimal source (3.2) leads to significant improvement that results in earlier time plateaux (i.e. at shorter time separations) for the effective pseudoscalar meson mass. Further details on the implementation of the method for computing pseudoscalar meson masses, decay constants and bag parameters are given in Appendix A. (a) improved plateau using a smeared field at the source and a local field at the sink (SL) -green crosses -compared to the correlators (LL) -magenta upward triangles -or (SS) -blue downward triangles; (b) further improvement is obtained with optimised field at the source and local field at the sink after tuning the parameter w. Blue points represent the values of M eff for different values of the parameter w; red full circles correspond to a case close to the optimal improvement; green crosses are for the w=1 case which corresponds to the (SL) correlator case shown in panel (a).

Computation of the b-quark mass and decay constants f B and f Bs
The determination of the b-quark mass, the decay constants f B and f Bs as well as for their ratio, f Bs /f B , is carried out by adopting the so-called ratio method presented in Refs. [18,17]. We refer the reader to these papers for a detailed presentation of the method. We will here discuss how our improved 2-point correlation functions lead to reduced systematic uncertainties. We start by recalling the main steps of the ratio method for the computation of the b-quark mass. The HQET suggests the asymptotic behavior, where M hℓ is the heavy-light pseudoscalar meson mass and µ pole h is the heavy quark pole mass. The key observation is that the static limit of appropriate ratios of the quantity in the (lhs) of Eq. (4.1) taken at nearby values of the quark pole mass is equal to unity. This knowledge can be exploited in order to compute the b-quark mass by interpolating between relativistic data obtained in the charm quark mass region and somewhat above it, and the infinite heavy quark mass limit. To this aim it is convenient to consider a sequence of heavy quark masses (µ h ) which have a fixed ratio, λ, between any two successive values: µ At each value of the lattice spacing we then build the following ratios, where the function ρ(µ (n) h , µ), that is known up to N 3 LO in perturbation theory [52,53,54,55,56], relates the MS renormalised quark mass (at the scale of µ = 3 GeV) to the pole mass: µ pole h = ρ(µ h , µ) µ h (µ). By construction, ratios of pseudoscalar meson masses at successive values of the heavy quark mass are expected to show small discretisation errors even for rather large values of µ h . At each value of µ (n) h we can thus perform a well controlled combined chiral and continuum fit on the ratios of Eq. (4.2) to extract the quantity y(µ h ) ≡ y(µ h , λ; µ u/d , a = 0). As an example of the quality of the fit we report in Fig. 4(a) the linear fit in µ ℓ of the data for y(µ (n) h ) at the largest value of the heavy quark mass. Relying on the well-known matching of heavy-light meson mass evaluated in QCD onto HQET, we have defined the ratio y(µ h ) in such a way that its dependence on µ h can be described by the fit ansatz that implements the constraint lim µ h →∞ y(µ h ) = 1. The fit parameters could be, in general, functions of log(µ h ). However in the range of the currently explored heavy quark mass values this logarithmic dependence can be safely neglected. Hence we approximate η 1,2 to constants whose value will be determined by the fit to the available ratio data. Data (with ρ from NLL order perturbative formulae) and fit are shown in Fig. 4(b). Finally the value for the b-quark mass can be computed from the chain equation .
Here on the one hand λ, K and µ (1) h are such that M hu/d (µ  of its role to Eq. (4.4) we will call it the triggering point. The quality of the linear fit in µ ℓ is shown in Fig. 5(a). 4 For the present analysis we use (µ (1) h , λ) = (1.05 GeV, 1.1784), for which Eq. (4.4) is satisfied for K = K b = 9. Many other choices of (µ (1) h , λ, K b ) could equally well be used. The b-quark mass result reads (4.5) In section 6 we discuss the error budget attached to this result. We have also verified that using M hs instead of M hℓ leads to fully compatible results for the b-quark mass, the difference being at the per mille level. We apply an analogous strategy to compute the pseudoscalar decay constant of the meson B s , f Bs and the ratio of the decay constants f Bs /f B . The appropriate HQET asymptotic conditions in the first two cases are where by 'constant' we denote some finite non-zero value. Based on QCD to HQET matching of heavy-light CL -phys. point  h ) against the renormalised light quark mass µ ℓ . Both fit ansätze are linear in µ ℓ and in a 2 . Empty black circle is our result at the physical u/d quark mass point in the continuum limit for both cases.
meson decay constant (f hℓ ) and mass (M hℓ ) we define the ratios The factor C stat A (µ * , µ h ), known up to N 2 LO in PT [58], provides the matching between the decay constant in QCD for a heavy quark mass µ h and its static-light counterpart in HQET (the arbitrary renormalization scale µ * of HQET cancels in the ratio above). From Eqs (4.8) and (4.9) we also form the double ratio (4.10) By construction the ratios z d , z s and ζ have an exactly known static limit equal to unity and show a smooth chiral and continuum combined behavior. As in the case of the y-ratios, this is a consequence of the fact that z d , z s and ζ are simply ratios of quantities calculated at nearby values of the heavy quark mass for which much of the discretisation errors cancel. Figs 6(a) and 6(b) are two examples illustrating the quality of the combined chiral and continuum fits for z s (µ (7) h ) and ζ s (µ (7) h ) respectively, at the largest heavy quark mass values used in the decay constant analysis. 5 6 In Figs 7(a) and 7(b) we show the dependence of z s (µ h ) and ζ(µ h ) on the 5 Note that at the largest value of the heavy quark mass, µ (8) h , which has been used in the b-quark mass analysis, our estimates of the pseudoscalar meson decay constants proved to be rather noisy. Hence, in the decay constants' analysis we decided to use data corresponding up to the next largest heavy quark mass value, µ h . 6 The rather high χ 2 /dof value in the fit of ζs(µ (7) h ) data ratio is not representative of the quality of the ratio fits performed in β = 4.20 h ) (Eq. (4.10)) against µ ℓ are shown in panels (a) and (b) respectively. In both cases ratios for the largest value of the heavy quark mass are reported (n = 7). Empty black circle is our result at the physical u/d quark mass point in the continuum limit.
inverse heavy quark mass, respectively. The fit ansätze we have used are polynomial fit functions in the inverse heavy quark mass analogous to the one specified in Eq. (4.3). For the case of the double ratio ζ(µ h ) we have also tried a linear fit in 1/µ h always implementing the static condition lim µ h →∞ ζ(µ h ) = 1.
Determinations of f Bs and f Bs /f B are obtained by means of the equations . (4.12) The (lhs) of the above equations are taken from the fits of Figs 7(a) and (b) respectively. Setting µ (cf. Eq. (4.5)) and having determined the values of f hs (µ h )] from a combined chiral and continuum fit, we finally get our results for f Bs and [f Bs /f B ] respectively. As expected the combined chiral and continuum fit for the quantity f hs (µ (1) h ) is smooth and shows tolerably small cutoff effects, as well as a very weak dependence on the light quark mass (see Fig. 5(b)). the present work. The fits performed in the present analysis are dominated by systematic uncertainties. Correlation matrices are not taken into account since they turn out to be affected by large uncertainties, and the χ 2 definition incorporates the contribution of priors. Furthermore, we would like to stress that in order to control systematic effects in the ratio analysis for all the physical quantities studied inn this work, we have repeated the whole procedure excluding the heaviest quark mass. We have treated the difference between the analyses as a systematic uncertainty in the final result. See also the discussion in Section 6. by means of which any use of the heavy quark pole mass is avoided. In analogy to Eq. (4.9) we can define the ratioz . (4.14) We then determine the value of f Bs by means of the equatioñ where we set M hs (µ Note that, as it can be seen from Fig. 8(a), the triggering point calculation of the quantity f hs √ M hs presents very small discretisation effects, and, though it is an accidental fact, it contributes to an accurate computation of the continuum limit. The fit ansatz used in fitting the data of Fig. 8(c) against the inverse heavy quark mass is of the same form as the one presented in Eq. (4.3). We anticipate here (see also Section 6) our finding that determinations of f Bs computed either via z s orz s ratios are fully compatible differing by less than 0.5%. We also note that ratios defined in terms of M hs (or M hℓ ), rather than the heavy quark pole mass, could be used for all the matrix elements discussed in the present paper.   We also need to estimate the triggering point ratio [f hs (µ h )] i.e. the value it takes after extrapolation to the continuum limit and the physical light quark mass. To this aim we make use of the useful observation [20,59] that by forming the double ratio of the decay constants [(f hs /f hℓ )/(f sℓ /f ℓℓ )] one can exploit the possibility for a large cancellation of the chiral logarithmic terms. One can get the triggering point value by multiplying the above expression by an appropriate estimate of the ratio of the K to π decay constants, (f K /f π ). For notation simplicity we define the quantity and we plot it against µ ℓ , see Fig. 9. We have used two fit ansätze. The first fit ansatz is linear in µ ℓ , the second one is suggested by the combined use of the SU(2) ChPT and HMChPT. They read In the fit based on HMChPT, we take for the parameterĝ the valueĝ = 0.61 (7) [60] obtained from the experimental measurement of the g D * Dπ coupling. We choose this value, instead of the HQET predictionĝ = 0.44 (8) [61], because we fit data that are close to the charm mass region and in order to conservatively include in the average the maximum spread resulting from the different ways of performing the chiral extrapolation of our data.
As it can be noticed from Fig. 9 discretisation effects are small. The two estimates for the triggering point ratio at the physical light quark mass are compatible within two standard deviations. We take their average value as our final result and we consider their half difference as a systematic uncertainty. In this computation we have used the result f K /f π = 1.193(5) by FLAG [62] coming from an average over lattice determinations using N f = 2 + 1 dynamical quark simulations. This value is completely uncorrelated with relevant determinations by ETMC. The latest PDG result for the same ratio could be an alternative choice. This differs by one standard deviation from the above one [2]. In order to to get a (conservative) estimate of this particular systematic uncertainty we consider the spread between the above value, f K /f π = 1.193 (5), and the one determined from N f = 2 dynamical quark simulations, f K /f π = 1.210(18) [62]. Finally, we sum in quadrature the two systematic uncertainties, that is the one coming from the fit ansatz choice and the other from f K /f π . Our result at the triggering point reads f hs f hu/d µ (1) h = 1.201 (7)(20), (4.19) where the first error is statistical and the second denotes the systematic uncertainty we have just discussed.

Computation of Bag parameters and ξ
For the evaluation of the four-fermion matrix elements on the lattice we use a mixed fermionic action setup where we adopt different regularizations for sea and valence quarks as proposed in Ref. [46]. The Mtm-LQCD action of the light quark flavor doublet that is used to generate unquenched gauge configurations can be written in the so-called "physical basis" in the form [42] S Mtm sea = a 4 The field ψ denotes a mass degenerate up and down doublet with bare (twisted) mass µ sea . The parameter M cr is the critical mass that one has to fix non-perturbatively at its optimal value [43,44] to guarantee the O(a)improvement of physical observables and get rid of all the unwanted leading chirally enhanced cutoff effects. In the gauge sector the tree-level improved action proposed in Ref. [40] has been used.
For valence quarks we use the OS regularization [63]. The full valence action is given by the sum of the contributions of each individual valence flavour q f and reads [46] S OS val = a 4 x,fq where the index f labels the valence flavors and M cr is the same critical mass parameter which appears in Eq. (5.1). We denote by r f and µ f the values of the Wilson parameter and the twisted quark mass of each valence flavor. This particular setup offers the advantage that one can compute matrix elements that are at the same time O(a)-improved and free of wrong chirality mixing effects [64]. These two properties have already proved to be very beneficial in the study of neutral Kaon meson oscillations [39,65,66,67]. For a detailed discussion about the choice of the action and its implementation for the calculation of the matrix elements we refer to Section 4 and Appendix A of Ref. [39].
The only significant difference concerning the four-fermion matrix elements computation of the present work with respect to Ref. [39] is that here we are dealing with heavy quarks (i.e. charm-like and heavier ones) rather than with the strange quark. With heavy quarks the signal to noise ratio deteriorates quickly as the time separation increases and moreover large time separations are necessary for the projection onto the ground state. We overcome both problems by employing smeared interpolating operators for the meson sources and in this way we are able to reduce the source time separation, T sep . The latter turns out to be less than half of the lattice   Table 7.

Ratio method for the bag parameters and ξ
We apply the ratio method to compute the bag parameters relevant to the neutral B-meson system. We introduce ratios that approach unity in the static limit. In HQET, the static limit of each of the five bag parameters is a constant. Hence, following a procedure analogous to that used for the b-quark mass and decay constants determination, we build ratios for each bag parameter computed at nearby values of the heavy quark mass that we define to be 7 ratios above. It incorporates the QCD evolution from the scale µ to the reference scale identified by the heavy quark mass and the HQET evolution from the same heavy quark mass scale to some arbitrary scale µ ⋆ (the dependence on which cancels in all ratios above), as well as the matching between the two theories. At NLL order the W matrix has a (3 × 3 ⊕ 2 × 2) block diagonal form. Note that at TL or LL order the SM bag parameters B . For more details on the QCD-HQET matching of the B-parameters see Appendix B. Moreover in Appendix C we show that using the ratio method (with QCD to HQET matching at TL order) one can numerically verify with good precision the relationship which connects, via the equations of motion, the operators O It has been noticed that, due to strong cancellations of systematic effects, it is convenient for the unitarity triangle analysis, to compute the SU(3)-breaking parameter ξ, whose knowledge fixes the ratio of the CKM matrix elements, |V td |/|V ts |. According to the philosophy of this paper we define the ratio (see Eqs .

(5.10)
Bag parameters defined at µ (1) h are determined as usual by a combined chiral and continuum fit that shows small discretisation effects. In Figs 11(a) and 11(b) we present the fit for the cases of B i cases HMChPT predicts a logarithmic chiral behaviour [69]. We have thus tried besides a linear fit ansatz, fit functions of the following type:  We now pass to the discussion of the dependence on the inverse heavy quark mass of the ω-ratios, evaluated in the continuum limit at the light quark physical point. As before we try polynomial fit ansatz of the form, using alternatively TL, LL, NLL order expressions for the matrix W . The relevant fits (for the case of W at NLL order) are illustrated in the five panels of Fig. 14  We follow the same procedure to determine the numerical value of the ratio B 1 and the parameter ξ. Our final results will be presented in the next section. In Fig. 15(a) we show the combined chiral and continuum fit for the bag parameter ratios at µ (1) h for four values of the lattice spacing. We have made use of two fit ansätze. The first one suggested by the HMChPT reads  h ) (Eq. (5.10)) against µ ℓ are shown in the panels (a) and (b), respectively. In both cases ratios for the largest value of the heavy quark mass are reported (n = 7). Linear fit ansatz in µ ℓ and fit to a constant are shown. Colored lines show the linear fit to the data. Empty black circle and full grey circle are the results, respectively, at the physical u/d quark mass point in the continuum limit. withĝ = 0.61(7) [60]. The second is a linear fit with no logarithmic terms. It can be seen from Fig. 15(a) that both fit ansätze lead to compatible results within half standard deviation. 8 We average over the two results and we take their half difference as a systematic error: We study the dependence of ξ(µ (1) h ) on µ ℓ by making use of the quantity R f defined in Eq. (4.16) in order to better control and test the impact of the logarithmic terms from our data. The quantity that we fit against the light quark mass is defined by The fit is illustrated in Fig. 15(b). We have tried a fit ansatz following SU(2) ChPT and HMChPT combining the formulae given in Eqs. (4.18) and (5.15). We have also tried a fit assuming only a linear dependence on µ ℓ . Following the same reasoning as in the case of the ratio of the decay constants we take our final result as the average over the two values coming from the respective fit ansätze. We also consider similar arguments to get the estimate of our systematic uncertainty. Hence our result at the triggering point reads ξ(µ  Linear fit ansatz in µ ℓ and fit to a constant are shown. Empty black circle and full grey circle are the results, respectively, at the physical u/d quark mass point in the continuum limit.

Summary of results and discussion
In this section we collect the results for the b-quark mass, the pseudoscalar decay constants of the B and B s mesons, the bag parameters for the full four-fermion operator basis that quantify the QCD effects in the B The first error (see the first entry for m b in the error budget Table 9) has been computed using the bootstrap method. It includes the statistical uncertainties on the pseudoscalar meson mass at the triggering point (typically less than 0.5%) and on the ratios (less than 0.1 − 0.2%), the uncertainty due to the quark mass renormalisation constant, Z P , in employing renormalised light and heavy quark masses, as well as the statistical uncertainties of the lattice scale. Notice that the pure statistical error of the pseudoscalar mass values computed from the relevant 2-point correlator functions lies typically at a sub-percent level. The second error includes: (i) systematic uncertainties of the lattice scale; (ii) The estimate of discretisation effects that has been obtained by repeating the whole analysis without including the data corresponding to the coarsest value of the lattice spacing; (iii) systematic uncertainties due to fit choices by (a) fitting the ratios y against the inverse heavy quark mass adding an extra cubic dependence in the fit ansatz (see Eq. (4.3)) and (b) excluding from the fit all the data from the ratio corresponding to the heaviest quark mass. Each of these checks in the fitting procedure leads to a small shift of 0.2-0.3% from the central value. (iv) We consider he systematic error due to the Λ QCD value as well as due to the half difference between results obtained employing in the running either N f = 2 or N f = 4 dynamical flavours.
In Table 9 we present a detailed description of the error budget giving the percentage error for each of the systematic sources of uncertainty. Our total error is computed by summing in quadrature all the above errors and it amounts to less than 3%. By averaging the two results above and including their half-difference as an additional systematic uncertainty we finally obtain Our current result for the b-quark mass is fully compatible with our previous determination in Ref. [17]. where the first error is due to statistical uncertainties estimated using the bootstrap method. The statistical error for the matrix element from which decay constants are computed is at a sub-percent level. The statistical uncertainties at the triggering point and those of the ratios are about 0.5% and 0.2%, respectively. It is also included the statistical uncertainty of the lattice scale. The second error refers to systematic uncertainties all added in quadrature. They are discussed in the following. (i) For dimensionful quantities the systematic uncertainty of the lattice scale (about 2%) is taken into account. (ii) We estimate the systematic uncertainty due to residual discretisation effects by repeating the whole analysis without including data from the coarsest value of the lattice spacing. (iii) The f Bs value has been obtained using the NLL approximation of the factors C stat A and ρ, (see Eq. (4.9)). 9 Had we used LL or TL approximation the total shift to the central value would be about 2 MeV which is well covered by the quoted errors. This is also included in the systematic error budget. A shift to the final value at the level of per mille is noticed if we exclude from the fitting procedure of ratios against the inverse heavy quark mass data corresponding to the heaviest quark mass. (iv) The second error in the decay constant ratio (6.5) as well as the third error in f B is due to the systematic uncertainty arising from the chiral fit at the triggering point -see the discussion in Section 4 and the result (4.19). In Table 9 we present in detail the full error budget. The alternative ratio method computation of f Bs based on the use of the HQET asymptotic behaviour expressed in Eq. (4.13) leads to an estimate which differs from the above of Eq. (6.4) by only 0.5%, but its total error is a bit larger due to statistical and systematic uncertainties coming from the meson mass, M hs . Our total uncertainty is given by the sum in quadrature of the above errors. For f Bs , f Bs /f B and f B the total error is 3.4%, 2.0% and 4.0%, respectively. In Table 9 a detailed description of the error budget is presented. We also notice that our present results are compatible within one standard deviation with the older ones of Ref. [17] which have been computed using local sources for the various propagator inversions.
As a by-product of our work, since the heavy quark mass value at the triggering point has been set equal to the charm quark mass, we have computed the decay constants for the D s and D mesons as well as their ratio. where the first error is statistical and the second one is systematic.
The results for the B-meson decay constants in Eqs. (6.4)-(6.6) are obtained in the N f = 2 theory and do not account for the dynamical sea quark effects of the strange and the heavier quarks. In Figs. 18 (a) to (d) we compare these results and the one for the b-quark mass of Eq. (6.3) with those obtained by other lattice studies using either N f = 2, N f = 2 + 1 or N f = 2 + 1 + 1 simulations. The good agreement observed among the various results in the plots provides a first indication that the systematic effects due to the quenching of the strange and charm quarks are smaller than present uncertainties from other sources.
More quantitatively, one can compare our results for the decay constants with the averages quoted by the FLAG working group from N f = 2 + 1 calculations [74]:  Table 9: 1 and ξ. Each of the numbers in the first row includes the uncertainties coming from the statistical errors of the correlators, the uncertainty of the quark mass renormalisation, the error of the combined chiral and continuum fits and the statistical uncertainty from the scale setting. In the second row we estimate the uncertainty due to the lattice scale systematics. In the third row we give an estimate of the systematic uncertainty related to discretisation effects. In the fourth row we display the systematic uncertainty connected to the fit of the ratios against 1/µ h . In the fifth row the systematic uncertainty at the triggering point fit is shown. The sixth row is for the systematic error coming from the uncertainty in the value of Λ QCD and the half difference between results obtained employing either N f = 2 or N f = 4 in the evolution from µ = 3 Gev to the b-mass point ∼ 4.3 GeV. The last row gives the total error. than other uncertainties. For the heavier charm quark, the quenching effect is expected to be even smaller in size and the comparison between the N f = 2 + 1 + 1 results of Ref. [75] with the N f = 2 and N f = 2 + 1 determinations presented in Fig. 18 does not show, indeed, any systematic deviation.
In the case of the charm quark, a parametric estimate of the quenching effect can be also derived using a perturbative argument. Since the sea quark contributions enter at the loop level, involving the exchange of at least two gluons, and are quadratically suppressed in the decoupling limit by the inverse quark mass, one expects them to be of order α s (m c ) 2 (Λ QCD /m c ) 2 for the charm quark, corresponding to less than 1%. Even though the accuracy of perturbation theory at the charm scale is limited, this estimate is of the same order of magnitude of the one indicated by the previous numerical comparison between lattice results. It supports the conclusion that for the decay constants this error is negligible with respect to other uncertainties.
The estimate of the partial quenching effects in the determination of the quark masses, like the b-quark mass in the present study, is more subtle, since it also depends on the details of the renormalization procedure. In our calculation, we have renormalized the quark mass with the non-perturbative RI-MOM method at the scale µ = 3 GeV and we have chosen to quote our final result in the MS scheme at the scale m b . Thus, a dependence on the number of dynamical flavors also appears in the matching between the initial and the final scheme and scale. This effect is quantified by the difference between the results in Eq. (6.1) and (6.2) and it is at the level of 1%. In this case, being not fully negligible with respect to other uncertainties, it has been taken into account in the systematic error quoted in Eq. (6.3). Table 10 we gather our results for the bag parameters of the full four-fermion operator basis. Results are given in the MS scheme of Ref. [25] at the scale of the b-quark mass m b of Eq. (6.3). Refs. [75,77,77,78,79]; (c) Refs. [75,77,78]; (d) Refs. [80,81,79]. (The results of Ref. [79] are still preliminary. Note also that the results in panel (a) indicated as HPQCD'11 and HPQCD'12, as well as those indicated as HPQCD'12 in panel (b), have both been obtained using N f = 2 + 1 (MILC) gauge ensembles but employ different valence quark regularisations.)   1, . . . , 5), renormalized in the MS scheme of Ref. [25] at the scale of the b-quark mass. See the text for discussion about the quoted errors.

Bag parameters: In
We consider our final values to be the average over the results obtained by using N f = 2 and N f = 4 in the evolution from the QCD renormalization scale (µ = 3 GeV) and the b-mass scale (∼ 4.3 GeV), while their difference, which is less than 1% (in the worst case), is taken as an additional systematic error. In the results of Table 10 the first error corresponds to the sum of the statistical errors (typically about 1%) on the correlators, the RCs uncertainty (which is the largest among the others and amounts to 2 − 3% depending on the case) and the fit error. The second error includes: (i) the uncertainties due to the possible choices of fit ansatz in 1/µ h and the maximum spread in the results induced by using the TL, LL or NLL order formulae for W in the ratios of Eqs (5.4) and (5.5). The latter uncertainty is for most of the cases rather small compared to the first error i.e. about 0.1% − 1%, with the exception of B  i ); (iii) the systematic uncertainty of residual discretisation effects, estimated by repeating the whole analysis without including data from the coarsest value of the lattice spacing. our systematic uncertainty coming from the value (iv) the systematic error due to the Λ QCD value as well as due to the half difference between results obtained employing in the running either N f = 2 or N f = 4 dynamical flavours. In Tables 11 and 12 we give the detailed error budget for B i , respectively. We notice that the results for B (d/s) i in Ref. [24], obtained in the quenched approximation and using Wilson-clover quarks are in the same ballpark as the present ones. A detailed comparison, for example on the impact of quenching, is not very meaningful because older quenched results have been obtained at rather large pion masses and at only one lattice spacing (a ∼ 0.1 fm), while the unquenched ones of this work have been extrapolated to the continuum limit and to the physical pion mass. The first error in both results is of statistical nature and it is due to the fitting procedures employed for the triggering point, ∼ 1% (in both quantities) and for the corresponding ratios (less than 1% for each one of them). The second error is our estimated uncertainty coming from employing different types of ansatz in fitting the ratios (less than ∼ 0.5%) as a function of the inverse heavy quark mass, as well as the uncertainty coming source of uncertainty (in %)  Table 11: Error budget (in % ) for B d i (i = 1, . . . , 5). First row includes the statistical uncertainties of correlators, of the fits' extrapolation and interpolation as well as the uncertainties of the RCs. The second row includes the systematic uncertainties due to discretisation effects. In the third row we give the estimates of the systematic errors by employing different choices of fit ansätze in fitting the ratios as a function of µ ℓ and the systematics related to the fit of the continuum value of the ratios against 1/µ h . We have also included the systematic uncertainty estimated by repeating our analysis without using data from the heaviest quark mass pair. In the fourth we show the systematic uncertainty from the fitting procedure at the triggering point. In the fifth row we display our systematic uncertainty coming from the value of Λ QCD and the half difference between results obtained employing in the running either N f = 2 or N f = 4 dynamical flavours. The last row gives the total error. source of uncertainty (in %) B   Table 11.
from using different orders for the QCD and HQET running and matching (less than 0.5% in both cases). The third error in ξ corresponds to the uncertainty due to possible fit ansatz choices at the triggering point discussed in Section 5, (see Eq. (5.18)). Our final uncertainty is taken as the sum in quadrature of the various uncertainties and it amounts to 2.0% for B (s) 1 and 2.5% for ξ. The error budget for both quantities is summarized in Table 9.
In Table 13 we collect our results for the quantities f Bq B (q) i with q = d, s and i = 1, . . . , 5. Bag parameters are expressed in the MS scheme of Ref. [25] at the scale of the b-quark mass. The error for each of these quantities is determined by combining the errors previously discussed on decay constants and bagparameters. For convenience we also give our results for the SM relevant quantities in which the bag parameters are expressed in the RGI scheme. In this case for the running of the coupling constant we take N f = 5 and   5). Bag parameters are expressed in the MS scheme of Ref. [25] at the scale of the b-quark mass. [27]. We find where the first error is statistical and the second one is systematic.
In Fig. 19 we compare our results for ξ, f Bd B (d) 1 and f Bs B (s) 1 , with bag-parameters expressed in the RGI scheme, with the ones obtained by other lattice collaborations.
In conclusion, in this work using smeared interpolating operators in computing the 2-and 3-point correlation functions: (i) we update our older continuum limit results [17] for the b-quark mass and the pseudoscalar decay constants for the B and B s pseudoscalar mesons for which local interpolating operators had been used in computing the 2-point correlation functions. Our final errors now are smaller though they are principally determined by uncertainties in the renormalisation constant for the quark mass and the lattice scale. Besides that, the fitting procedure of our ratios against the inverse heavy quark mass towards the exactly known static limit has gained more accuracy with respect to our previous work thanks to the higher precision in calculating the ratios; (ii) we have presented results, extrapolated in the continuum limit and to the physical light quark mass point, about the bag-parameters for the full four-fermion operator basis that control the neutral B-meson oscillations, as well as results for the parameter ξ and the quantities f Bd B Finally, preliminary results for the b-quark mass employing the ratio method on N f = 2 + 1 + 1 dynamical quark configurations generated by ETMC have already been presented [84]. ETMC will present in the near future a complete analysis of the physical observables studied in the present paper but employing N f = 2+1+1 gauge ensembles.  [82,19]; (c) Refs. [82,19].

A Optimised smeared interpolating operators
Optimised two-point correlators can be built using optimized fields Φ opt as suggested in Eq. (3.2) which we report here for the reader's convenience Without loss of generality we consider in practice correlators with an optimized field at the source and a local field at the sink, which read 10 The choice of the normalization time x 0 , where by construction C opt-L 2 (x 0 ) = 1, is in principle arbitrary and only affects the optimal value of w and its behaviour as a function of the mass parameters. The x 0 values used at each β in this study are gathered in Table 14. They have been chosen so as to yield a reasonably smooth behaviour of the optimal w-value as a function of the heavy quark mass µ h , which in turn eases its numerical search. We recall that here Φ L and Φ S always carry the quantum numbers of a pseudoscalar density.
The optimal values of w are those for which earliest Euclidean time projection on the ground state can be achieved. We now describe the procedure we have followed. We consider the correlator given in ps . In Table 14 we collect the Euclidean time intervals (∆x 0 ) (j) that we used at each β. Given these time intervals, at each β and each µ h we vary w with a sufficiently fine resolution (typically around 0.05), then for all w-values we compute the maximal spread among the pseudoscalar mass estimates, M , is taken at some larger Euclidean time (but not too large so as to avoid introducing too much statistical noise).
Using maximally twisted mass lattice fermions, the decay constant of a pseudoscalar (non-singlet) meson made out of two valence quarks with masses µ 1 and µ 2 , regularized with Wilson parameters r 1 = −r 2 , is evaluated via the (Ward-identity based) formula ( [41,42,59]: 10 We should note that it would also be possible to construct an optimised correlator formed by a local source Φ L and a sink given by the operator defined in Eq. (3.2). 11 For instance, at β = 3.80 there are five time intervals: (∆x0) (   where we have set a = 1 and P =q 1 γ 5 q 2 is the pseudoscalar density. Using optimised interpolating fields, f ps is extracted from the following formula (in lattice units) where the correlator C opt−opt 2 is given by In Fig. 20 we illustrate an example on the pseudoscalar decay constant computation using optimised and smeared interpolating operators.
This computation implies the employment of optimised P † A 0 and P † P 2-point correlation functions, for i = 1 and i = 2, 3, 4, 5, respectively. At β = 3.80 we have compared results for the (bare) bag parameters B i (i = 1, . . . , 5) obtained through Eq. (A.8) with the ones coming through the use of smeared interpolating operators. We have found out that within our current statistical errors it is hardly noticed any difference on the plateau values of the bag parameters. In the two panels of Fig. 21 we illustrate two examples supporting the above numerical observation for the cases of B i with i = 1, 2, respectively. A similar behaviour has been observed also for the bag parameters B 3 , B 4 and B 5 . In the two figures we have also included data corresponding to the case where local sources have been employed. In this last case, as it might be expected, the plateau quality results problematic if heavy quark masses are to be employed as we do in the present study. Based on the above observations we have thus opted for using smeared interpolating fields in the computation of the bag parameters throughout this work. Finally, we note that statistical errors on pseudoscalar meson masses, pseudoscalar decay constants as well as on four-fermion operator matrix elements have been evaluated using the jackknife procedure. With 16 jackknife bins for each configuration ensemble we have verified that autocorrelations are well under control. In order to take into proper account the cross correlations we compute statistical errors on the fit results using the bootstrap method and for that we employ 1000 bootstrap samples produced by independent gauge configuration ensembles.

B QCD-HQET matching of B-parameters
According to the HQET scaling laws, each B-parameter scales with the inverse heavy quark mass µ h as a constant up to some perturbative logarithmic corrections of order ∼ 1/ log (µ h /Λ QCD ). These corrections are expected to be tiny in the range of heavy quark masses we are dealing with. This point is checked, and the corresponding systematic uncertainty is quantified, by matching at different (TL, LL, NLL) perturbative orders the QCD B-parameters to their counterparts in HQET, thereby yielding B-parameters with a welldefined static limit and correspondingly smaller ( O(α(µ h )), O(α(µ h ) 2 ), O(α(µ h ) 3 ) ) logarithmic corrections to their leading 1/µ h -behaviour. It turns out (see the error budget discussion in Section 6) that the impact of logarithmic corrections to the power-like 1/µ h -behaviour on the chain equations (5.9) and (5.10), through which B-parameters are determined, is small (O(1%)) compared to other uncertainties in the calculation. It is only to this well controlled extent that a perturbative QCD-to-HQET matching of the quantities of interest, carried out by using the formulae given in the present Appendix, enters in our ratio method computation.
In this Appendix, we consider the QCD four-fermion operators in the SUSY basis of Eq. (1.7). The corresponding HQET operators have exactly the same form but with the relativistic heavy quark field h replaced by the (infinitely) heavy quark field of HQET. In the HQET due to the heavy quark spin symmetries, the operator O 3 is related to O 1 and O 2 by the relation O 3 = −O 2 − 1/2 O 1 . Nevertheless, we find it convenient to work with the redundant basis of five operators, which includes O 3 , in order to deal with squared 5x5 evolution and matching matrices in HQET as in QCD.
The relation between QCD B-parameters evaluated at the heavy quark mass µ h and their counterparts in HQET can be expressed as:B where B i (µ h ; µ) denotes the B-parameters in QCD renormalized at the scale µ andB i (µ h ; µ * ) are the HQET B-parameters renormalized at the scale µ * . The latter satisfy the heavy quark scaling laws with logarithmic corrections dictated by their renormalization group evolution, which are therefore easy to take into account when applying the ratio method. The matrix W (µ * , µ h , µ) can be decomposed as follows: The matrix U (µ h , µ) encodes the full QCD evolution from the scale µ to the scale µ h for the five ∆B = 2 B-parameters.Ũ (µ * , µ h ) is the corresponding evolution matrix in HQET, from µ h to µ * . Finally, C (µ h ) provides the matching from HQET to QCD at the common scale µ h .
At next-to-leading-log (NLL), the explicit expression of W reads where the superscript T stands for "transposed". The corresponding leading-log (LL) expression of W is obtained by setting J B =J B = c • β 0 is the leading order coefficient of the QCD beta function, β 0 = 11 − 2N f /3 B is the scheme-independent one-loop anomalous dimension matrix (ADM) of the QCD B-parameters. It is obtained by combining the one-loop ADM of the corresponding four-fermion operator γ (0) with the one-loop ADM of the pseudoscalar density γ (0) P as: • J B is the scheme-dependent two-loop ADM of the B-parameters in QCD. It is obtained from the ADMs of the four-fermion operator J and of the pseudoscalar density J P through where in the MS scheme with N f = 2 active flavors one obtains J P = 8134/2523. The expression of J in the SUSY basis and in the MS scheme of Ref. [25], with N f = 2, is •J B is the scheme-dependent two-loop ADM of the B-parameters in HQET. It is obtained from the ADMs of the four-fermion operatorJ and of the axial densityJ A , both in HQET, as whereJ A = 1037/2523 − 28π 2 /261 in the MS scheme with N f = 2. The HQET two-loop ADM of the four-fermion operators has been computed in Ref. [87]. In the SUSY basis, it reads • c (1) B is the LL order coefficient of the matching matrix between HQET and QCD. By combining the matching relations for the four-fermion operators, the axial and the pseudoscalar densities, one finds A δ i1 δ j1 − 2c In this appendix it will be shown that using the ratio method and data taken from relativistic quark simulations one can verify the validity of Eq. (C.1) in the infinite heavy quark mass limit. Following the discussion in Ref. [39], see Eq. We will show numerically, using just TL formulae for connecting matrix elements in QCD to their HQET counterparts, that the following equation is satisfied in the static limit with good precision. Following the familiar procedure that we have used throughout this paper we define suitable ratios of R In Figs 22(a) and (b) we display the combined chiral-continuum fit against the light quark mass for R 2 and R 3 , respectively, at the triggering point (µ (1) h ). In both fits we have used a linear fit ansatz in µ ℓ as well as in a 2 .
For each pair of nearby heavy quark masses we get the continuum limit result at the physical point of the ratios defined in Eqs (C.5) and (C.6) through a combined chiral-continuum extrapolation. In  3 (µ h ) (right panel) against µ ℓ for the largest value of the heavy quark mass considered in this work (n = 7). The full black line corresponds to a linear fit ansatz in µ ℓ and a 2 , while the dashed one corresponds to the continuum limit curve in the case of a linear fit in a 2 without dependence on µ ℓ (constant fit). In both panels colored lines correspond to the linear fit ansatz in µ ℓ and a 2 .
we illustrate this extrapolation for r , respectively at the largest value of heavy quark mass used in this work (n = 7).
In Fig. 24 we show the dependence of r respectively, as the number of N s increases. We observe that for N s > 25 the asymptotic value compatible with zero within one standard deviation is reached.