$B^0$-$\bar{B}^0$ mixing: matching to HQET at NNLO

We compute perturbative matching coefficients to the Heavy Quark Effective Theory representation for the QCD effective local $\Delta B=2$ Hamiltonian that determines the mass difference in $B^0-\bar{B}^0$ system of states. We report on the results at NNLO in the strong coupling constant for matching coefficients of two physical operators in HQET. Our results provide the firm confirmation that the recent NLO sum rules analysis of the bag parameter $B_q$ is stable with regard of inclusion of higher order radiative corrections. As a by-product of our calculation we give a fully analytical solution for the one loop QCD-to-HQET matching problem: we present the explicit formulas for the renormalization of four quark operators of the full bases in both QCD and HQET and the expressions for matching coefficients in a closed form.


Introduction
The accuracy of theoretical predictions for the neutral meson mixing within the Standard Model (SM) has been steadily improving during last years [1,2,3]. The main reason for this progress is due to better numerical precision achieved for the numerical values for hadronic matrix elements in the lattice simulations (e.g. [4,5]). Nevertheless, recently some new results with a competitive level of accuracy appeared in the domain of analytical computation based on sum rules approach [6,7]. The latter calculation has become possible due to the technical advance of integral computation at three loops [8]. At present, the high precision of experimental material provides for good opportunities for searches of physics beyond SM [9], and these new theoretical predictions for the mixing within SM are important since the search for new physics depends heavily on the accurate knowledge of numerical values of low energy parameters [10,11]. The mixing of the neutral B mesons is dominated by the top-quark contribution and hence looks directly to the ultraviolet (UV) new physics having long distance effects under control. The concept of effective theories allows for getting better precision for the theoretical predictions due to separation of scales [12,13,14] that essentially improves on old electroweak results [15]. Clearly, the most crucial point of such an improvement is perturbation theory (PT) corrections at the lowest scales involved in the analysis. This motivates our research -we compute next-to-next-to leading order (NNLO) corrections to the matching of QCD to heavy quark effective theory (HQET) that is applicable at the scales of order few Λ QCD . Some partial results have been already published earlier [16].
The paper is organized as follows. In the next section we briefly introduce notation just for the paper to be self-contained (the details can be found in [16]) and give our main results. In Sect. 3 we describe the technique of the computation. In Sect. 4 we briefly discuss implications of our calculation for phenomenology. In summary section we give our conclusions.

Operators and Wilson coefficients
In the SM the B 0 -B 0 transition is described by a nonlocal expression of quark scattering at loop level. The most important corrections to the leading electroweak term are the contributions of strong interactions. They can be computed within QCD perturbatively as the relevant scale is of order the meson mass m B and is much larger than QCD infrared scale Λ QCD .
In this section we repeat some formulas from [16] to introduce notation. More details are given in [16].

QCD: below m W
A relevant scale for the description of B 0 -B 0 physics is around the B-meson mass m B that is kinematically saturated by the b-quark mass m b . After integrating out particles of the SM that are heavier than b-quark the effective Hamiltonian for the process of B 0 -B 0 mixing reads H = C(m W , m t , µ, α s (µ))Q(µ) (2.1) with Q(µ) =q L γ α bq L γ α b(µ) . (2. 2) The quantity Q(µ) is a local renormalized operator with ∆B = 2. One can choose a low normalization point µ ∼ m b for the operator Q while large logarithms of the scales ratio, ln(m W /m b ), in the Wilson coefficient C(m W , m t , µ, α s (µ)) can be regularly resummed with renormalization group (RG) techniques. At present the coefficient C(m W , m t , µ, α s (µ)) is known at the next-to-leading order (NLO) of the expansion in the strong coupling constant that gives the accuracy of a few percent [12,13,14]. The renormalization properties of the operator Q are a bit peculiar but well understood and intensively discussed in the literature. The point is that the calculation are usually and dominantly performed in dimensional regularization and one has to close up the algebra of Dirac gamma matrices which is infinitely dimensional in general d-dimensional space. This is an outstanding problem and has been much discussed the literature [17]. Dimensional reduction has been used for the computation of corrections to ∆F = 1 Hamiltonian [18]. A naive dimensional regularization prescription (NDR) is introduced in [19]. As a result of using dimensional regularization for the integral evaluation one requires a special treatment of Dirac structure and, eventually, the extension of the operator basis. For the four quark operator in question, i.e. Q(µ) in eq. (2.2), the procedure was discussed by Buras et al. [20]. The method has been further analyzed in [21,22]. Note that the anomalous dimensions of baryon or three quark operators have been computed earlier within a similar approach [23,24]. A clear presentation of the techniques is given in [25]. Thus, the effective Hamiltonian in eq. (2.1) should,in fact, contain additional evanescent operators and should explicitly read as where E is a general notation for a string of evanescent operators in QCD. Their matrix elements vanish but their presence in the basis influences the renormalization pattern and, therefore, evolution of physical operators Q(µ). The renormalized operator Q(µ) depends on the choice of the evanescent ones. By choosing E = E + a Q one obtains a different renormalized physical operator Q (µ). At the one loop level one gets the relation between the two Physical predictions stay independent of the choice of evanescent operators. We discuss how it is achieved later in the text.

HQET: below m b
Since the scale m b is still QCD perturbative, m b Λ QCD , one can go further low in scales and remove the explicit dependence on m b from the matrix element or the mixing amplitude at low energy. This is achieved by using HQET [26,27,28].
The low scale for the operators involved in the mixing used to be necessary for the lattice computation, however at present there is sufficient power for lattice simulations directly at the scale m b . In a computational framework within analytical methods one matches the theory of QCD on to HQET where considerable technical simplifications occur for subsequent computation of sum rules. Thus the matching is an unavoidable part of the whole computation. Though there is an approach based on calculation at m b [29]. The results need update as the definition of the operator has been different.
The heavy quark expansion (HQE) of the operator Q(µ) goes The field h + annihilates the heavy quark in HQET (moving with the four velocity v), and h − creates the heavy antiquark (again moving with the four velocity v), which is a completely separate particle in HQET framework. In HQET one can define its own set of physical and evanescent operators (see, [16]). The general basis is and q is a light fermion which is usually a chiral one, q ≡ q L . A choice for a basis in HQET is an antisymmetrized product of transverse gamma's, Then γ n ⊥ is a notation for an antisymmetrized product of n transverse gamma matrices. We sometimes call the number n a rank of the product and, therefore, of the corresponding operator.
For our two loop calculation it is more convenient to change the operator basis in the physical sector from the standard The operators {O p , O l } do not mix under renormalization at the one-loop level [30]. The expression of operators {O l , O p } through the basis operators {O n , O n } is The matching pattern of the QCD operator Q(µ) at µ = m b then reads where dots denote the contribution of evanescent operators. We define a PT expansion of the matching coefficients C l,p (m b ) as where N is a number of colors for the SU (N ) gauge group. Note the difference with the {O l , O s } basis that is traditionally used at NLO in the literature (2.14) In the present paper we have computed the NNLO contributions to the coefficients {C l , C p } which is the main result of the paper.
The coefficient C Here I 0 is one of the master integrals of the computation that reads Note that the PT expansion in HQET goes over the n l -flavored coupling constant since there are just n l flavors in the low energy theory. In QCD we have in addition an active heavy quark b and, therefore, n l + 1 flavors.
For the number of colors N = 3 and the number of massless flavors n l = 4, the numerical values of the expansion coefficients are One sees that the convergence of PT series for quantity C p (m b ) in the renormalized coupling constant α s (m b ) with n l = 4 is marginal.
The coefficient C The coefficient of the n l structure is different from the one given in our early paper [16] for two reasons. First, the mixing of operators {O l , O s } has not been accounted for, and second, the expansion of Γ( ) at NLO has not been done to a necessary order in ε (up to O( ), in fact) that has caused a finite shift proportional to π 2 in the NNLO coefficient. So, the naive non-abelianization [33] works moderately well: the β 0 terms are about 2 times larger than those without β 0 .
The results for the matching coefficients correspond to the normalization point One sees that the NNLO corrections are rather large. Their calculation is crucial for estimating the reliability of the results at NLO in PT. We discuss some physical implications of the obtained results later.

Description of computation
QCD operators can be expanded in 1/m in terms of HQET operators: where Q(µ) is the column of renormalized QCD operators, O(µ) is the column of renormalized HQET operators, and C(µ) is the matrix of matching coefficients. For example, the leading matching coefficients C for heavy-light quark currents were calculated in [34,33,35,36].
Here we follow the same procedure but for the 4-quark operators. First we calculate the bare matching coefficients where Z(µ) andZ(µ) are the matrices of renormalization constants in QCD and HQET (we omit 1/m corrections). This is done by equating the on-shell matrix element of Q 0 where Z os Q is the on-shell heavy-quark field renormalization constant, Z os q is the on-shell light-quark renormalization constant, and Γ 0 is the vertex function of Q 0 , to C 0 times the corresponding on-shell matrix element of O 0 . If all light quarks (including the charmed quark c) are considered massless, all loop corrections to the HQET quantitiesZ os Q ,Z os q ,Γ 0 vanish in dimensional regularization because they contain no scale.
The 2-loop results for Z os Q [37] and Z os q [35] are known. We calculate Γ 0 up to 2 loops using the reduce package recursor [38], similarly to [33]. The basis of antisymmetrized products of γ ⊥ allows us to project onto individual HQET operators easily. The calculation is done in an arbitrary covariant gauge; we check that the on-shell matrix element (3.3) is gauge invariant (Z os Q and Z os q are gauge invariant up to 2 loops). If we keep only the (gauge-invariant) subset of diagrams in which only b quark and the light quark belonging to the same color-singlet current interact, we successfully reproduce the 2-loop matching coefficients [33,35] of the heavy-light currents with the Dirac structures 1, / v, γ ⊥ . The quantities Γ 0 , Z os Q , Z os q are calculated via the n f -flavor bare coupling g (n f ) 0 ; we re-express them via the renormalized α (n f ) s (µ). The renormalization constant matrix Z(µ) is also expressed via α (n f ) s (µ). However, the renormalization constant matrixZ(µ) is expressed via α (n l ) s (µ), n l = n f − 1. We have to use the decoupling relation between α (n f ) s (µ) and α (n l ) s (µ), see, e. g., an introductory review [39]. It is most convenient to perform matching at µ = m, the on-shell mass of the b quark. Then with T F = T (R) = 1/2 for fermions in fundamental representation. We need this O(ε) term because the 1-loop QCD matrix element contains poles in 1/ε. The computation has been done in leading logs in [40,41] where LO anomalous dimension for O l has been found. It happens to be equal to that of the simple product of two heavy light currents. In higher orders this factorization does not necessarily persist. The standard result for the coefficients C l,p at NLO has been obtained in [31,30,32].
In our paper [16] we have calculated the NNLO contribution of the leading order in n l only. This allows one to estimate the full results for the 2-loop matching coefficients using the approximate method of naive non-abelianization [33]. Presently available techniques allow for the analytical calculation of any number of fermionic bubbles (e. g. Chapher 8 in [28], also available as [42]) that can be converted into the β 0 -dominance estimates of the matching coefficients at any order of perturbative expansion. While the estimate is technically feasible, the quantitative validity of the approximation for a phenomenological analysis is not immediately clear (e. g. see some discussion in [43]).
General description of matching calculations given above is well known. In our case of matching four-quark operators of QCD onto four-quark operators in HQET there is a subtlety of using dimensional regularization caused by the presence of spurious operators that formally vanish in four-dimensional space. Therefore one organizes the basis of operators in both QCD and HQET in physical and evanescent sectors. In QCD the operators are symbolically Q = {Q, E} (do not confuse the concrete operator Q and general notation for the whole set in QCD). In HQET the operators are O = {O l , O p , e i }. The matching is by necessity a relation between whole basis sets of the operators in QCD and in HQET and it reads in general matrix form Here CO denotes a matrix multiplication of the matrix of matching coefficients C and the string (one dimensional matrix) of operators of the basis. The renormalization pattern of the operators in eq. (3.5) (including evanescent ones) is with Z,Z being the renormalization matrices and Q B , O B are bare images of renormalized operators in QCD and HQET. One obtains for the matching in eq. (3.5) The matching coefficients can be found by taking the matrix elements on shell in HQET from the both sides of this relation. And for HQET one has because all loops are scaleless and vanish in dimensional regularization. This is independent of whether the operator is physical one or evanescent one. Let O| j be an operation of projecting on a particular state P j (operator) and one chooses a complete system of P j for the HQET basis such that O B i | j = δ ij . For the basis of four quark operators this operation has a simple realisation: one takes traces in both Dirac strings. Since the bare operator of the basis has a structure of a direct product O B n = A n ⊗ B n , symbolically, The matching coefficients C become This is a working formula. The quantity (Q B p | j ) is basically a bare matching coefficient of the bare operator Q B p from the QCD basis to the bare operator O B j in an HQET basis. The bare matching coefficient depends on one scale m b and is represented by loop integrals on shell. At NNLO they are two loop integrals.
At one loop level the renormalization matrix in QCD reads where z EQ is obtained from the requirement that matrix elements of the renormalized is related to an anomalous dimension of the physical operator Q. It is independent of details of the basis. The quantity z QE = −T F describes the admixture of the evanescent operator E to the physical one. The quantity z EE is anomalous dimension of evanescent operators. Note that there are more than one evanescent operator in the basis, however for our computation we need only one independent combination of those, and z EE is irrelevant altogether. At the two-loop level we need only one additional entry to the renormalization constant Z which is z The value of z QQ is reconstructed from a two-loop anomalous dimension of the physical operator Q. The anomalous dimension of the operator Q(µ) is and n f = n l + 1 is the number of flavors. For the perturbative expansion of the anomalous dimension of the form the coefficient γ 1 in the basis described reads [20] In QCD one needs only one evanescent operator for the whole computation of matching in two-loops, just that one which admixes to Q at NLO even if it can be composed of several basis operators Q n .
In HQET the renormalization matrix for relevant operators looks similar to that one in QCDZ but now the quantity z OO is a 2 × 2 matrix in the subspace of physical operators {O l , O p }.
In HQET two evanescent operators are necessary for renormalization of physical operators, one for O l and one for O p . For computing the matching coefficients the whole set of evanescent operators is relevant or at least the one what Q can match onto at NLO (in fact, there is one operator that is multiplied by the poles of the matching coefficient and the other with only finite parts). It shows that, indeed, in general the whole operator basis should be matched onto.
As for bare coefficients we need two-loop values for C(Q → O l , O p ) and one-loop values for C(Q → j) for any operator O l , O p , e i and C(E → l, p). One more ingredient is Z Q on shell at NNLO from [37]. The calculation contains rather delicate cancellations of infinities (poles in ) and fixing the finite parts according to sophisticated conventions. To give just an example, the renormalization matrix Z is expanded, by convention, over the renormalized coupling α (n l +1) s (µ) while the renormalization matrixZ is expanded, by convention, over the renormalized coupling α (n l ) s (µ). In the course of computation of matching coefficients the poles in cancel.
We have computed the necessary quantities. The quantity z QQ can be extracted from earlier calculations through anomalous dimensions of Q at two loops. The matrix z (2) OO is related to the anomalous dimension matrix of the physical pair (l, p). One of the entry has been considered in [44] where the result for the anomalous dimension of O l has been given. The whole basis of operators was not explicitly specified in [44]. We do not compute the corresponding anomalous dimension independently. It can be extracted from our calculations though. In fact, one can extract only the difference γ 1 −γ 1 that reads Assuming the value for γ 1 from eq. (3.16) with n f = n l +1 we have extracted the anomalous dimensionγ 1 of the HQET operator O l The coefficient n l agrees with our paper [16]. The whole expression disagree with the results quoted in ref. [44]. Numerically we find for n l = 4 and QCD with N = 3. The irony of life is that the two quantitiesγ (1) andγ (1) | G are rather close numerically though the analytical expressions for them are quite different. As for the calculation at one loop level the whole program of renormalization of bases of four quark operators in QCD and HQET and matching onto one another has been explicitly realized in a closed form. The corresponding formulas are given in the Appendix.

Implications for phenomenology of the mixing
The splitting between the mass eigenstates of B 0 -B 0 system is determined, in the SM, by the non-diagonal matrix element of the effective Hamiltonian (2.1) or, more specifically, if the concrete calculations are made within dimensional regularization, by the expression (2.3). After the coefficient function have been determined in QCD perturbation theory the task is the computation of hadronic matrix elements of four-quark operators, Q, E. While for evanescent operators E this task can be trivially accomplished since by construction B 0 |E(µ)|B 0 = 0, an accurate determination of numerical values for physical operators is a genuine QCD-bound-states problem. The nonperturbative techniques required for the computation can be QCD sum rules and direct lattice simulation. The sum rules results for the matrix element B 0 |Q(m b )|B 0 have been presented earlier [45,29,46,47] within different approximation schemes for the Green functions used in the OPE. The obtained precision has been rather limited by modern standards though. Till recently the lattice analysis with fully relativistic heavy quarks was impossible due to insufficient computational power, therefore the matching to HQET (2.5) has been a necessity. The NLO results for matching coefficients to both HQET and to the lattice representation for the operators have been obtained more than a quarter of a century ago. The technical breakthrough with computation of three loop HQET integrals in ref [8] allowed for a NLO analysis of the mixing matrix element using sum rules in HQET. In our recent paper [6] we have computed the bag parameter B d for B 0 d −B 0 d mixing at the next-to-leading order of perturbative expansion for matching coefficient and for the Green functions entering the sum rules analysis. To evaluate the matrix element of the mixing we use a vertex (three-point) correlation function [48]. The analysis uses the splitting of the whole Green function necessary for the calculation within OPE and for the sum rule approach [48,46,47] into factorizable and non-factorizable parts. It happens that the non-factorizable part starts only at NLO of perturbation theory and turns out to be small. These features allow for getting a numerical result of high precision for the bag parameter.
The techniques have been also used for other four quark operators in [7] 1 . The computation of ref. [7] is pinning down the important uncertainty for lifetime differences of the heavy mesons with both b and c quarks. It is rather a complete analysis but it is limited to only NLO of the perturbation theory for HQET Green functions. Future experimental data may require even more accurate theoretical predictions. For obtaining still better theoretical accuracy, the NNLO perturbative corrections to matching coefficients can be useful. The first step in this direction has been done in our recent paper [16] where the NNLO corrections proportional to the n l factor have been computed that allowed us to perform an approximate evaluation of the coefficients within the naive non-abelianization (β 0 dominance) approach [33]. In the present paper we have computed the full NNLO results for matching coefficients. They read numerically where a s = α (4) s (m b )/(4π). Our calculation of matching coefficients is an important step in the program of NNLO description of mixing within analytical methods of computation. From our results we see that NLO approximation for matching coefficient may not be reliable since the NNLO corrections are large. On the other hand, we also see that the values of corrections to the coefficients by themselves do not lead to immediate physical conclusions. One has to add corrections to corresponding Green functions which determine the OPE for sum rule analysis. Thus, the NLO correction to C l requires NLO correction to the correlator of the operator O l for consistent NLO analysis (see ref. [16] for more details). In case of the O p operator the corresponding correlator of the O p operator can be taken at tree-level approximation for obtainig the NLO result for the QCD matrix element. With the NNLO correction included to C p , one needs the NLO correction to the correlators with insertion of O p ; this turns out to be feasible [8] and it is even available for some cases [7]. In case of the NNLO corrections to C l one has to compute the NNLO corrections to the correlator (4.4) as well that leads to four-loop integrals, which are currently beyond known technology.
Note that with at NNLO there will appear new evanescent operators which we do not specify explicitly. However, one can try to choose the basis of evanescent operators such that C (2) l becomes smaller or even vanishing (that will change the numerical value of C (2) p accordingly) and one can hope to obtain a chance to avoid the necessity of the computation of NNLO corrections to the correlator (4.4).
The other possibility to get the most of our NNLO calculation of matching coefficients even without having the NNLO results for the correlator (4.4) is to perform a direct comparison between physical quantities. While most observables in QCD contain nonperturbative contributions, there are, in fact, some observables for which one can show that one may construct a perturbative relation between the two observables. In such a case the perturbative expansion acquires an immediate physical meaning and statements about the size of coefficients as well as on the convergence of the perturbative series become meaningful. To this end, while for the individual observables large coefficients may appear (depending on the definition of the operator matrix elements), in a relation between these observables the large coefficients my cancel, once the matrix elements are defined in the same way in both observables.
As we pointed out above, the bag factor B B turns out to be small and hence to a good approximation we expect a perturbative relation between ∆M d and the B-meson decay constant f B . In fact, it is known that the matching coefficients of the axial current are also large [33], and we may expect that in a direct comparison of the axial vector current with ∆M d a well-behaved perturbation theory results.
The matching coefficient of the axial current to HQET interpolating operator is [33,35] with a rather sizable coefficient at NNLO. The HQET matching (2.5) has a great deal of arbitraryness in distributing the contributions between coefficients and operators. One concrete choice of fixing the MS-scheme for the definition of the operators is a current standard. If we call Q(m b ) to be our physical quantity determining δm (again up to the freedom of redefinition the matrix elements and the coefficients in QCD but we set this aside now) then the expansion in HQET is not unique in a variety of aspects: choosing a physical pair of the operators, changing renormalization scheme, etc. Also clearly, the freedom of the definition of evanescent operators is not only the choice of a physical pair but deviation from minimality. By adding a physical operator to an evanescent one with a coefficient vanishing at = 0 we obtain different renormalized physical operators, see (2.4) and therefore a different matching coefficient.
While physics, i.e. predictions for observables, does not depend on this reshuffling of the operator bases, the independence restores only after adding matrix elements computed within the same scheme up to same accuracy. The magnitude of the correction to a particular coefficient is not an invariant characteristics of PT expansion and one should really collect Wilson coefficients and matrix elements together which is difficult in QCD since there is no any quantitative scheme for computing hadronic matrix elements. Lattice simulation may be an exception, but then the perturbative short distance part of the matrix element cannot be easily identified.
The full NNLO analysis of the mixing seems to be feasible, and the most intriguing part is the possibility to find a basis of evanescent operators in which the NNLO correction to the matching coefficient of the operator O l turns out to be small. This assertion deserves to be validated in future work.
Note also that the two-loop anomalous dimensions of the operators {O l , O p } are not very important quantitatively. The point is that the difference of scales is not large and summation of logarithms of the scale differences is not crucial numerically. Indeed, matching is done at m b with m pole b ∼ 4.8 GeV (e.g. [49,50]) while the Green functions necessary for sum rules analysis are computed at the scale of the order of w 0 ∼ 1 GeV. The leading logs of the form (α s (w 0 ) ln(w 0 /m b )) n , with n > 0 can be summed with the leading anomalous dimensions of the operators {O l , O p } which are known while the subleading logs of the form α s (w 0 )(α s (w 0 ) ln(w 0 /m b )) n with n ≥ 1 are not large and can be retained in an expanded form.
Here we check that our results for matching coefficients allow for a reasonable perturbative expansion of physical parameters. Recall that the matrix element in QCD is represented by the expression and therefore B q (µ) gives a relation between physical parameters ∆m and f 2 B (we factor out the redefinition freedom of the operator Q(µ) in QCD considering it fully under control in perturbation theory). The perturbation theory series for the relation between ∆m and f 2 B is unambiguous and can depend only on N c and α s (with all possible reservations that they are not completely PT quantities). The realization of this idea is given in eq. (2.18) of ref. [16] in the form Here B l,s are bag parameters for the HQET operators O l,s (see, ref. [16] for more detail). The perturbative corrections to the parameter B l has the form [6] represents an (yet unknown) NNLO correction to B l from HQET sum rules (compare to the parameter X used in [6] for estimating higher order contributions). We have extracted the correction to the parameter B s from ref. [7] with the result Substituting our results for the coefficients C l , C p and the expression for C J , we finally obtain the relation between physical observables in the form The quantity x (2) l emerges from NNLO corrections to HQET sum rules and it is expected to be in the range of NLO correction which is just of order unity (compare assumptions on the parameter X in [6], |X| < 20). The perturbative expansion in eq. (4.12) has reasonably small coefficients as we had expected. We see that all large coefficients in C l , C p , C J mutually cancel each other and numbers of order unity survive in the final expression rending the PT expansion to be reliable.
The conjecture about the above pattern of perturbative expansion is quite feasible and can be explicitly checked once x (2) l is determined. Note that the idea of re-expressing the physical quantities through one another with a resulting PT factor is rather old and has been widely used. It is often applied as well to observables which are not quite fully perturbatively related (see, e.g. [51,14].
Note that the lattice results are based on NLO analysis for the matching coefficients between continuum and lattice representations of the operators. While that matching is different, our computation shows that PT corrections at NNLO to the matching coefficients can be important at this level of precision.
If there is no particularly large NNLO contribution to the sum rule determination of x (2) l , our numbers show that the precision of the matrix element at the level of few percents even in the presence of NNLO matching can be obtained. The large corrections are mainly hidden in C 2 J . In other words, the single out of B parameter is very efficient physically since it has a reasonable perturbative expansion a posteriori and a bunch of perturbative correction to the matrix element simply reproduces the correct value of f B .
Note in passing that with the NNLO accuracy of the leading term one may need to account also for nonleading terms of HQE [52] while small power corrections have been accounted in [47].
In general, up-to-date lattice results turn out to be more precise than sum rule estimates for classical quantities like decay constants [53,54]. Only due to a special structure of the observables occurring in the B 0 −B 0 mixing one can still use the sum rules to obtain predictions that are still competitive with lattice computations.
Extension of our results to the case of B 0 sB 0 s has been discussed in [16]. As the mass m s is not large it can be taken into account in expansion in m s /w 0 (an interesting example is given in [55]). We have shown that the strange quark mass appears in non-factorizable quantities only at NLO level while the leading order contributions are hidden in factorizable parameters like f Bs [53]. Numerically the parameter (m s /w 0 )α s (w 0 ) [56] corresponds to NNLO level but achieving this accuracy requires only one loop calculations. With formulas given in the Appendix it is rather a straightforward computation which we are going to present in future publications.

Summary
We have calculated NNLO corrections to matching coefficients necessary for the analysis of mixing in B 0 −B 0 system, i.e. for the calculation of bag parameters in sum rules at three-loop level in HQET [6,7]. The NNLO corrections happen to be large, however, to rather large extent they cancel the large NNLO corrections in the matching coefficients for the axial current that determines the B-meson leptonic decay constant, f B . The relation between experimentally measured quantities ∆m and f 2 B turns out to be rather well behaved as a perturbative series up to NNLO. This observation gives a strong ground to our estimate of the uncertainties for the QCD bag parameter B q along the lines of ref. [6,16].
We also discuss possible ways of getting invariant physical predictions from our results independent of the introduction of evanescent operators in the physical sector due to renormalization. We have constructed and present the completely analytical framework for analyzing the renormalization and matching of four quark operators at one loop level.

Acknowledgment
A.G. is grateful to Siegen University for hospitality; his work has been partially supported by the Russian Ministry of Education and Science. This work is supported by the DFG Research Unit FOR 1873 "Quark Flavour Physics and Effective Theories".
For brevity we'll write it as M 0 = T 1 Γ ⊗ Γ where The 1-loop matrix element is We are interested only in the UV 1/ε divergent terms. Therefore we may treat all quarks as massless and set all external momenta to 0. Then we need some IR regulator, say, replacing all massless denominators by the ones with some non-zero mass, or a hard IR cutoff in euclidean momentum space. Such a regulator will be implied, not written explicitly. The MS quark field renormalization constant is where ξ is the gauge fixing parameter. Averaging over directions of the loop momentum k in the integrands, we easily obtain The matrix element (A.4) is gauge invariant: where only the UV 1/ε divergences are kept in the 1-loop terms. For the operator the only difference is that the color factors: Now we specifically consider the operators is the antisymmetrized product of n γ matrices. We have [21] γ µ Γ n γ µ = (−1) n (d − 2n)Γ n , Using these relations twice, we can re-write the matrix elements (A.6), (A.8) as where <O> in the right-hand side are the Born-level matrix elements. We are interested only in the UV 1/ε divergent terms. Therefore we may treat light quarks as massless, set their external momenta to 0 and external residual momenta of HQET (anti-) quarks to 0. An IR cutoff is implied. The MS HQET field renormalization constant is Averaging the integrands over k directions (in particular, using (k · v) −2 = −(d − 2)(k 2 ) −1 [33,57]), we easily obtain The matrix element (B.2) is gauge invariant: Setting γ µ = / vv µ + γ ν ⊥ and using the (d − 1)-dimensional versions of (A.11), we can rewrite the matrix element (B.4) as −<Õ n+2 > + (d + n − 2)<Õ n+1 > + n(3d − 2n − 4)<Õ n > − n(d − n)(d + 2n − 4)<Õ n−1 > − n(n − 1)(d − n + 1)(2d − n − 2)<Õ n−2 > + n(n − 1)(n − 2)(d − n + 1)(d − n + 2)<Õ n−3 > ,