{'fraction_non_alphanumeric': 0.03628746754779325, 'fraction_numerical': 0.013275902761387775, 'mean_word_length': 4.4157213612398145, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anomalous Dimension Matrix for Radiative and Rare Semileptonic B Decays up to Three Loops

We compute the complete O(alphas^2) anomalous dimension matrix relevant for the b to s gamma, b to s gluon and b to s l+l- transitions in the standard model and some of its extensions. For radiative decays we confirm the results of Misiak and Munz, and of Chetyrkin, Misiak and Munz. The O(alphas^2) mixing of four-quark into semileptonic operators is instead a new result and represents one of the last missing ingredients of the next-to-next-to-leading-order analysis of rare semileptonic B meson decays.


Introduction
Inclusive radiative B decays place very important constraints on models of physics beyond the Standard Model (SM). The present experimental world average for the branching fraction of the radiative B → X s γ decay is [1] BR(B → X s γ) exp = (3.34 ± 0.38) × 10 −4 , while the most recent SM prediction is [2,3] BR(B → X s γ) th = (3.70 ± 0.30) × 10 −4 .
The experimental error is rapidly approaching the level of accuracy of the theoretical prediction. The main limiting factor on the theoretical side resides in the perturbative QCD calculation and is related to the ambiguity in the definition of the charm quark mass in some two-loop diagrams containing the charm quark [2]. To improve significantly the present Next-to-Leading-Order (NLO) QCD calculation, one would need to include one more order in the strong coupling expansion, and compute at least the dominant NNLO effects: a very challenging enterprise, which seems to have already captured the imagination of some theorists [4]. The present calculation of the branching ratio of B → X s γ consists of several parts that are worth recalling. Perturbative QCD effects play an important role, due to the presence of large logarithms L = ln m b /M W , that can be resummed using the formalism of the operator product expansion and renormalization group techniques. The main components of the NLO calculation, which aims at resumming all the next-to-leading O(α n s L n−1 ) logarithms, have been established more than six years ago. They are i) the O(α s ) corrections to the relevant Wilson coefficients [5,6,7], ii) the O(α s ) matrix elements of the corresponding dimension-five and six operators [8], and iii) the O(α 2 s ) Anomalous Dimension Matrix (ADM) describing the mixing of physical dimension-five and six operators [9,10,11,12]. After the O(α s ) matrix elements of some suppressed operators have been computed last year [3], the NLO calculation is now formally complete. Higher order electroweak [13,14] and non-perturbative effects [15,16] amount to a few percent in the total rate and seem to be well under control. However, as the actual measurements adopt a lower cut on the photon energy -typically around 2 GeV -the total rate is not accessible experimentally. In fact, unless the photon energy cut lies well below 2 GeV, the measured rate depends significantly on the non-perturbative structure of the B meson, and in particular on the Fermi motion of the bottom quark inside the B meson. The ensuing theoretical uncertainty is significant [14] and is generally included in the reported experimental error, since the experimental measurements are extrapolated to conventionally defined total rates.
Nearly all the ingredients of the NLO QCD calculation involve a considerable degree of technical sophistication and have been performed independently by at least two groups, sometimes using different methods. However, the most complex part of the whole enterprise, the calculation of the two-loop dimension-five [11] and of the three-loop dimensionsix [12] O(α 2 s ) ADM, has never been checked by a different group. The main purpose of this work is to present an independent calculation of the O(α 2 s ) ADM governing the b → sγ and b → sg transitions.
The rare semileptonic decay B → X s ℓ + ℓ − represents, for new physics searches, a route complementary to the radiative ones. Experimentally, the exclusive modes have been recently observed for the first time [17], and we now also have a measurement of the inclusive branching fraction [18]. More progress is expected in the future from the B factories. Because of the presence of large logarithms already at zeroth order in α s , a precise calculation of the B → X s ℓ + ℓ − rate involves the resummation of formally nextto-next-to-leading O(α n s L n−2 ) logarithms. The Next-to-Next-to-Leading-Order (NNLO) QCD calculation of B → X s ℓ + ℓ − has required the computation of i) the O(α s ) corrections to the corresponding Wilson coefficients [7] and ii) the associated matrix elements at O(α s ) [19,20]. Moreover, it involves iii) the O(α 2 s ) ADM, but the operator basis must be enlarged to include the semileptonic operators characteristic of the b → sℓ + ℓ − transition.
The O(α 2 s ) mixing between four-quark and the relevant semileptonic operators is one of the last missing pieces of a full NNLO analysis of rare semileptonic B meson decays. Since we include the semileptonic operators in the O(α 2 s ) calculation, we are able to close this gap. Notice however, that a formally complete NNLO calculation of the B → X s ℓ + ℓ − rate would also require the knowledge of the O(α 2 s ) self-mixing of the four-quark operators [21], of the O(α s ) matrix elements of the QCD penguin operators, and of the renormalizationgroup-invariant two-loop matrix element of the vector-like semileptonic operator.
We perform the calculation off-shell in an arbitrary R ξ gauge which allows us to explicitly check the gauge-parameter independence of the mixing among physical operators. To distinguish between infrared (IR) and ultraviolet (UV) divergences we follow [22] and introduce a common mass M for all fields, expanding all loop integrals in inverse powers of M. This makes the calculation of the UV divergences possible even at three loops, as M becomes the only relevant internal scale and three-loop tadpole integrals with a single non-zero mass are known. On the other hand, this procedure requires to take into account insertions of non-physical operators, as well as of appropriate counterterms.
We have so far emphasized the inclusive modes, as they are amenable to a cleaner theoretical description. However, one should not underestimate the importance of the exclusive B meson decays like [23,24] and B → ρℓ + ℓ − . A thorough study of the exclusive channels can yield useful additional information in testing the flavor sector of the SM. These processes have received a lot of theoretical interest in recent years and their accurate description will also benefit from a firm understanding of higher order perturbative corrections.
The ADM we have computed can be used in analyses of new physics models, provided they do not introduce new operators with respect to the SM. This applies, for example, to the case of two Higgs doublet models [2,6,25,26], and to some supersymmetric scenarios with minimal flavor violation. See for instance [26,27]. On the other hand, in left-rightsymmetric models [26,28] and in the general minimal supersymmetric SM [29], additional operators with different chirality structures arise. In many cases one can exploit the chiral invariance of QCD and use the same ADM, but in general an extended basis is required.
Our paper is organized as follows: in Section 2 we recall the relevant effective Lagrangian and list all the dimension-five and six operators that will be needed in the calculation. In Section 3 we review the renormalization procedure and explain how the operator renormalization matrix can be extracted from the matrix elements at higher orders. The actual two-and three-loop calculation is described in Section 4, while the results for the ADM are presented in Section 5. We provide some intermediate results which can be useful in a number of related applications in an Appendix. They include the complete O(α s ) and the relevant entries of the O(α 2 s ) operator renormalization matrices.

The Effective Lagrangian
Let us briefly recall the formalism. We work in the framework of an effective low-energy theory with five active quarks, three active leptons, photons and gluons, obtained by integrating out heavy degrees of freedom characterized by a mass scale M ≥ M W . In the leading order of the operator product expansion the effective off-shell Lagrangian relevant for the b → sγ, b → sg and b → sℓ + ℓ − transition at a scale µ is given by Here the first term is the conventional QCD-QED Lagrangian for the light SM particles. In the second term V ij denotes the elements of the CKM matrix and C i (µ) are the Wilson coefficients of the corresponding operators Q i built out of the light fields. In our case it is useful to divide the local operators Q i entering the effective Lagrangian into five different classes: i) physical operators, ii) gauge-invariant operators that vanish by use of the QCD×QED Equations Of Motion (EOM), iii) gauge-variant EOM-vanishing operators, and iv) so-called evanescent operators that vanish algebraically in four dimensions. In principle, one could also encounter v) non-physical counterterms that can be written as a Becchi-Rouet-Stora-Tyutin (BRST) variation of some other operators, so-called BRSTexact operators. However, they turn out to be unnecessary in the case of the O(α 2 s ) mixing of the operators Q i considered below. See also [11].
Neglecting the mass of the strange quark and the electroweak penguin operators which first arise at O(α), the physical operators can be written as [30] where the sum over q and ℓ extends over all light quark and lepton fields, respectively. e (g) is the electromagnetic (strong) coupling constant, q L and q R are the chiral quark fields, F µν (G a µν ) is the electromagnetic (gluonic) field strength tensor, and T a are the color matrices, normalized so that Tr(T a T b ) = δ ab /2. Notice that since QCD is flavor blind it is not necessary for our purposes to consider the analogues of Q 1 and Q 2 involving the up instead of the charm quark.
The physical operators given above consist of the current-current operators Q 1 and Q 2 , the QCD penguin operators Q 3 -Q 6 , the magnetic moment type operators Q 7 and Q 8 , and the semileptonic operators Q 9 and Q 10 , relevant for the b → sℓ + ℓ − transition. We have defined Q 1 -Q 6 in such a way that problems connected with the treatment of γ 5 in n = 4 − 2ǫ dimensions do not arise [10]. Consequently, we are allowed to consistently use fully anticommuting γ 5 in dimensional regularization throughout the calculation.
The gauge-invariant EOM-vanishing operators can be chosen to be [6,7] where the sum over f runs over all light fermion fields, while D µ and ← D µ denotes the covariant derivative of the gauge group SU(3) C × U(1) Q acting on the fields to the right and left, respectively. Notice that the set of operators Q 1 -Q 16 closes off-shell under QCD renormalization, up to evanescent operators [30,7].
In order to remove the divergences of all possible one-particle irreducible (1PI) Green's functions with single insertion of Q 1 -Q 10 we also have to introduce the following gaugevariant EOM-vanishing operators where G a µ denotes the gluon field, and we have used the abbreviations G µ = G a µ T a and d abc = 2Tr({T a , T b }T c ).
It is important to remark that the EOM-vanishing operators introduced in Eqs. (5) and (6) arise as counterterms independently of what kind of IR regularization is adopted in the computation. However, if the regularization respects the underlying symmetry, and all the diagrams are calculated without expansion in the external momenta, non-physical operators have vanishing matrix elements [31]. In this case the EOM-vanishing operators given in Eqs. (5) and (6) play no role in the calculation of the mixing of physical operators. If the gauge symmetry is broken this is no longer the case, as diagrams with insertions of nonphysical operators will generally have non-vanishing projection on the physical operators. Since our IR regularization implies a massive gluon propagator, non-physical counterterms play a crucial role at intermediate stages of the calculation.
As far as the evanescent operators are concerned, another eight operators are needed in order to find the O(α 2 s ) mixing of the physical operators Q 1 -Q 10 . Following [12,10] they can be defined as

Renormalization of the Effective Theory
Our aim is to study the renormalization properties of the physical operators Q 1 -Q 10 introduced in Eq. (4). Upon renormalization, the bare Wilson coefficients where the renormalization constants Z ij can be expanded in powers of α s as Following the standard MS scheme prescription, Z ij is given by pure 1/ǫ l poles, except when i = 25-32 and j = 25-32. In the latter case, the renormalization constant is finite, to make sure that the matrix elements of the evanescent operators vanish in four dimensions [32,33]. The calculation of an effective amplitude A eff , also involves the matrix element Q i ≡ F |Q i (µ)|I of the operator Q i between a initial state I and a final state F , which is renormalized by the usual coupling, mass and wave function renormalization factor characteristic of the operator, Q i → Z(Q i ). The renormalized effective amplitude is therefore given by where Z(Q i ) R denotes the matrix element of the operator Z(Q i ) after performing coupling, mass and wave function renormalization. Clearly, it is also possible to define the operator renormalization constantZ ij from the relation between unrenormalized and amputated Green's functions via Z( In this case, one simply has Z ij = Z −1 ij . In general Z(Q i ) will not be proportional to Q i . For example, in many of the EOMvanishing operator introduced in Eqs. (5) and (6) one has two different terms, only one of which has a factor of m b . Correspondingly, the m b renormalization of the operator is where Z m b denotes the mass renormalization constant of the bottom quark, and Q ′ i is the Figure 1: Some of the two-loop 1PI diagrams we had to calculate in order to find the O(α 2 s ) mixing of the complete set of operators Q 1 -Q 32 .
The product on the right-hand side of Eq. (10) must be finite by definition at any given order in α s . Therefore, requiring the cancellation of UV divergences we can extract Z (k) ij order by order. The result, up to third order, reads where the superscript (k) always stands for the k-th order contribution in α s .
If we leave aside the complication that in general Z(Q i ) will not be proportional to Q i , and write symbolically Z(Q i ) R = Z i Q i B the above relations can be rewritten in terms of bare quantities. Up to third order in α s we obtain The first line in Eqs. (13) recalls the familiar result that the one-loop renormalization matrix is given by the UV divergences of the one-loop matrix elements, after performing wave function and possibly coupling and mass renormalization. For example, in the case of the operators Q 1 -Q 6 , one has Z i = Z 2 q with Z q denoting the wave function renormalization constant of the quark fields, and Eqs. (13) take a particularly simple form, which upon expansion in α s reproduces the classical results derived more than ten years ago [32].
For a given set of operators and knowing the QCD renormalization constants, the solution of the above systems of linear equations requires the calculation of a sufficient number of Green's functions for different external fields with single insertions of the operators Q i . In our case, in order to determine the complete Z  Figure 2: Some of the three-loop 1PI diagrams we had to calculate in order to find the mixing of the four-quark operators into Q 7 -Q 10 at O(α 3 s ). Fig. 1. As we are interested in a subset of the three-loop ADM, we have actually calculated only the three-loop b → sγ and b → sg amplitudes involving insertions of Q 1 -Q 6 -see Fig. 2. We have calculated the complete off-shell amplitudes up to terms proportional to external momenta squared. By using the EOM it is therefore straightforward to extract the mixing into Q 7 -Q 10 . Notice that the results for the Z ij . The normalization of the physical operators adopted in Section 2 has been chosen [7] in such a way that the power of α s in Z ij is equal to the number of loops of the contributing diagrams. For instance, without the factor 1/g 2 in Q 7 -Q 10 , as in the standard normalization adopted in [2,3,20], both one-and two-loop diagrams contribute to the O(α s ) mixing matrix, because of the O(α s ) two-loop mixing of four-quark into magnetic operators. This choice simplifies both the implementation of the renormalization program and the resummation of large logarithms, since the redefinition enables one to proceed for b → sℓ + ℓ − in the same way as in the b → sγ and b → sg case.
In a mass independent renormalization scheme Z ij is µ-independent. This allows to check the renormalization of two-and three-loop matrix elements. The right-hand sides of the Eqs. (12) and (13) receive contributions from irreducible two-and three-loop diagrams as well as one-and two-loop counterterms. The µ-dependence is different in each case and governed by the n-loop factor (µ 2ǫ ) n . The UV structure of the k-th term is therefore given by where M (n,l) denotes the 1/ǫ l pole of the sum of all n-loop contributions. Expanding in powers of ǫ we find the following set of equations which have to be fulfilled to get a µ-independent Z ij up to three-loop order: This system of equations provides us with a powerful check of the renormalization of twoas well as three-loop diagrams. Notice that the locality of UV divergences also places some constraints on the renormalization matrix itself. We will turn to this point later on.

The Calculation of the Operator Mixing
In the renormalization of QCD and QED at higher orders the standard method of extracting the UV divergence structure of a Feynman integral is to perform the calculation with massless propagators. However, if one uses massless propagators to compute three-point or higher Green's functions one might generate spurious IR infinities which, in dimensional regularization, cannot be distinguished from the UV divergences one seeks. There exist several methods [34] to overcome this problem, but they are generally quite involved and not suitable to the automated evaluation of a large number of diagrams. In the approach of [22] the so-called IR rearrangement is performed by introducing an artificial mass. For the calculation of the renormalization constants this means that we can safely apply Taylor expansion in the external momenta after introducing a non-zero auxiliary mass M for each internal propagator, including those of the massless vector particles. The auxiliary mass regulates all IR divergences and the renormalization constants can be extracted from the UV divergences of massive, one-scale tadpole diagrams, that are known up to the four-loop level [35]. Following references [11,22], the starting point of our procedure is the exact decomposition of a propagator: Here k is a linear combination of the integration momenta, p stands for a linear combination of the external momenta, and m denotes the mass of the propagating particle. If we assume that the dimensionality of the operators in our effective theory is bounded from above, and we apply recursively the above decomposition a sufficient number of times, we will reach the point where the overall degree of divergence of a certain diagram would become negative if any of its propagators were replaced by the last term in the decomposition. We are then allowed to drop the last term in the propagator decomposition, as it does not affect the UV divergent part of the Green's function after subtraction of all subdivergences. This algorithm can be also simplified by the following observation [22]. The terms containing powers of the auxiliary mass squared in the numerators contribute only to UV divergences that are proportional to those powers of M 2 . The latter are local after the subtraction of all subdivergences, and must precisely cancel similar terms originating from integrals with no auxiliary mass in the numerators. Since the decomposition of Eq. (16) is exact, no dependence on M 2 can remain after performing the whole calculation. This observation allows one to avoid calculating integrals that contain an artificial mass in the numerator. Instead of calculating them, one can just replace them by local counterterms proportional to M 2 which cancel the corresponding subdivergences in the integrals with no M 2 in the propagator numerators. Nevertheless, the final result for the UV divergent parts of the Green's functions are precisely the same as if the full propagators were used.
The counterterms proportional to M 2 in general do not preserve the symmetry of the underlying theory, specifically they do not have to be gauge-invariant. Fortunately, the number of these counterterms is usually rather small, because their dimension must be two units less than the maximal dimension of the operators belonging to the effective theory. For instance, in QCD only a single possible gauge-variant operator exists that fulfills the above requirement. It looks like a gluon mass counterterm, and cancels gauge-variant pieces of integrals with no M 2 in the numerators. To ensure that our renormalization procedure with the fictitious gluon and photon mass is valid, we have checked explicitly the full MS renormalization of QCD and QED up to the three-loop level, finding perfect agreement with the results given in the literature [36]. In our case, beside the term in Eq. (17), we also have M 2 counterterms of dimension-three and four, some of which explicitly break gauge invariance: where A µ denotes the photon field. As already mentioned in Section 2, another side effect of our IR regularization is that we have to consider insertions of non-physical effective operators in our calculation. Let us explain this point in more detail. Non-physical counterterms generally arise in QCD calculations, but the projections of their matrix elements on physical operators vanish unless the underlying symmetry is broken at some stage. Due to the exact nature of the decomposition Eq. (16), the UV poles of the diagrams obtained by our method are correct after the subtraction of all subdivergences. However, the UV poles related to subdivergences and their subtraction terms both depend on the finite parts of certain lower loop diagrams, which in our approach are not necessarily correct and do not comply with the usual Slavnov-Taylor identities. For instance, the introduction of the IR regulator invalidates the argument that guarantees vanishing on-shell matrix elements for the non-physical operators. One therefore expects non-negligible contributions to the counterterms from all possible operators with appropriate dimension. Consequently, all EOM-vanishing operators, gauge-invariant or not, and in general even BRST-exact operators must be included in the operator basis. The "incorrect" subdivergences are present in both counterterm and irreducible diagrams, but they cancel in their sum, provided the calculation is carried out in exactly the same way. The operator renormalization constants calculated in this way are correct for all the operators in the complete basis.
The large number of diagrams which occurs at higher orders makes it necessary to generate the diagrams automatically. For the evaluation of the ADM presented here all diagrams have been generated by the Mathematica [37] package FeynArts [38], which provides the possibility to implement the Feynman rules for different Lagrangians in a simple way. We have adapted it to include the effective vertices induced by the operators Q 1 -Q 32 . We have processed the FeynArts output using two independent programs. In one case the output is converted into a format recognizable by the language Form [39]. The group theory for each graph as well as the projection onto all possible form factors is performed before the integrals are evaluated. The very computation of the integrals is done with the program package MATAD [40], which is able to deal with vacuum diagrams at one-, two-and three-loop level where several of the internal lines may have a common mass. The calculation of the tadpole integrals in MATAD is based on the so-called integration-byparts technique [41]. The second program is entirely a Mathematica code, which for the three-loop integrals uses the algorithm described in detail in [22].

The Anomalous Dimension Matrix
The anomalous dimensions γ ij defined by can be expressed in terms of the entries of the renormalization matrix Z ij as follows In a mass independent renormalization scheme the only µ-dependence of Z ij resides in the coupling constant. In consequence, we might rewrite Eq. (20) as where β(ǫ, α s ) is related to the β function via The finite parts of Eq. (21) in the limit of ǫ going to zero give the anomalous dimensions. Expanding the anomalous dimensions and the β function in powers of α s aŝ we find in accordance with [22] up to third order in α s : On the other hand the pole parts of Eq. (21) must vanish. From this condition one obtains relations between single, double and triple 1/ǫ poles of the Z ij , which constitute a useful check of the calculation. In agreement with [22] we find (1,1) , Having summarized the general formalism and our method, we will now present our results for five active quark flavors. For completeness we start with the regularization-and renormalization-scheme independent matrixγ (0) , which is given bŷ While the matrixγ (0) is renormalization-scheme independent,γ (1) andγ (2) are not. In the MS scheme supplemented by the definition of evanescent operators given in Eq. (7) we obtain The question marks stand for unknown entries that we have not calculated. As far as the remaining entries are concerned, our results for the one-and two-loop mixing of Q 1 -Q 6 agree with those of [10], and therefore also with previous results [9] that were obtained in a different operator basis [42]. We also confirm other well-established results, like the two-loop mixing of Q 1 -Q 6 into Q 7 -Q 10 [43]. See [7] for the conversion to our basis. For the two-loop mixing of Q 7 and Q 8 we confirm the results of [11], and for the entries γ and γ (2) 8i , that is the three-loop mixing of four-quark operators into magnetic ones, those of [12]. Neither of the latter results have been checked before, with the exception of the self-mixing of Q 7 [44].
We recall that the ADM reported in [12] refers to so-called effective coefficients. The relation with our results will be discussed below. Finally, the entries γ (2) 9i which contain terms proportional to the Riemann zeta function ζ 3 , are entirely new. The ADM entries involving non-physical operators can be readily obtained from the renormalization matrices provided in the Appendix. Finally, let us mention here once again that for a NNLO analysis of rare semileptonic B decays one would also need to know the self-mixing of Q 1 -Q 6 , which will be discussed in a separate communication [21].
It might be useful to recall explicitly the relation between the ADM in our basis and the ADM in an operator basis where Q 7 -Q 10 are not rescaled by 1/g 2 . The latter is frequently used for phenomenological applications [2,3,12]. The Wilson coefficients in that basis, C i (µ), are given by while the coefficients in the expansion in powers of α s of the corresponding anomalous dimensions γ ij , take the following form: ij , for i = 1-6 , and j = 7-10 , Notice that the perturbative expansion of the anomalous dimensions γ ij , that govern the evolution of C i (µ) already starts at zeroth order in α s , whereas γ ij contains no such terms.
In the case of radiative B decays it is sometimes convenient to use effective Wilson coefficients C eff i (µ) [45], defined in such a way that the leading order b → sγ and b → sg matrix elements are proportional to C eff 7 (µ) and C eff 8 (µ), respectively. In particular, dropping the semileptonic operators which are irrelevant for radiative B decays, the effective Wilson coefficients in the operator basis of [12] are given by In the MS scheme with fully anticommuting γ 5 one has y (7) = (0, 0, − 1 3 , − 4 9 , − 20 3 , − 80 9 ) and y (8) = (0, 0, 1, − 1 6 , 20, − 10 3 ). In consequence, the first two coefficients in the perturbative expansion of the effective ADMγ eff , that governs the evolution of C eff i (µ) read which are in perfect agreement with [12].

Summary
We have recalculated the complete O(α 2 s ) mixing relevant for the NLO analysis of radiative B decays in the SM and some of its extensions, confirming all previous results. In addition, we have also calculated the O(α 2 s ) three-loop mixing between four-quark and semileptonic operators relevant for the NNLO analysis of rare semileptonic b → sℓ + ℓ − transitions. This is a new result and provides one of the last missing pieces in the NNLO calculation for B → X s ℓ + ℓ − . We plan to complete the calculation of the NNLO ADM required for this process and to study its phenomenological implications in a separate communication.
The calculation involves the UV divergences of diagrams up to three loops. Our results have been subject to several cross-checks, following from: i) the locality of the UV divergences, ii) the independence of the ADM from the external states used in the calculation, iii) the completeness of our operator basis, iv) the gauge-parameter independence of the mixing among physical operators, and v) the absence of mixing of non-physical into physical operators. We have also reproduced the full MS renormalization of QCD and QED up to the three-loop level.

A.1 One-and Two-Loop QCD Renormalization Constants
The renormalization of the conventional QCD Lagrangian containing a massive bottom quark proceeds as usual. First, we introduce the renormalized fields and variables via where the subscript B denotes bare quantities and u a are the ghost fields. The gaugeparameter ξ is kept unrenormalized. This is legitimate, because the non-renormalization of the gauge-parameter is guaranteed by the usual Slavnov-Taylor identity, which is unaffected by the IR regularization adopted for the Yang-Mills theory. On the other hand, our IR rearrangement requires the introduction of the gauge-variant subtraction in Eq. (17), which can be interpreted as a counterterm for a fictitious gluon mass M.
Using the notation introduced in Eq. (9), the MS renormalization constants at one-loop order take the following form where C A = 3 and C F = 4/3 are the quadratic Casimir operators of SU(3) C . As usual N f stands for the number of active quark flavors. Our result for Z (1,1) M agrees with the expression given in [22]. At the two-loop level the poles of the MS renormalization constants are given by and Except for Z M , which have never been given explicitly, our renormalization constants agree with the results in the literature [46], if one bears in mind that the original papers contain some typing errors. We have also calculated the three-loop renormalization constants [36], but we do not report them here, as they are not needed in our calculation.

A.2 The Complete Operator Renormalization Matrix
The general structure of the operator renormalization matrix iŝ where P = 1-10 denotes the physical operators, N = 11-24 the EOM-vanishing operators, and E = 25-32 the evanescent operators. Throughout this section we set N f = 5. The mixing of non-physical into physical operators must vanish at all orders in α s . This is in fact only a requirement on the ADM, but we have seen in Eq. (24) that the one-loop renormalization matrixẐ (1,1) is proportional toγ (0) , and therefore at one-loop it implies the vanishing ofẐ (1,1) NP andẐ (1,1) EP . Sinceγ (0) for the physical operators can be found in Eq. (26), it is sufficient to give here only the non-physical parts ofẐ (1,0) andẐ (1,1) . By definition, the only non-vanishing parts ofẐ (1,0) areẐ The 6 × 4 block in the upper left corner ofẐ (1,0) EP agrees with the expression for the upper 6 × 4 block ofĉ given in Eq. (46) of [10].
The one-loop mixing of physical into non-physical operators is described byẐ The 4 × 6 block in the upper left corner ofẐ (1,1) PE agrees with the expression for the 4 × 6 block in the upper left corner ofb given in Eq. (45) of [10].
At one-loop we have moreover the mixing among EOM-vanishing operators, given bŷ and the mixing among evanescent operators, which readŝ Similarly to what happens in the case ofẐ PE not all entries ofẐ (2,1) PE are needed to find the O(α 3 s ) ADM of physical operators considered in this article. Needless to say, the mixing of physical into EOM-vanishing operators, described byẐ (2,1) PN , is not required to determine the mixing of physical operators at the three-loop level. However, some entries are important to verify the O(α 2 s ) mixing of magnetic into non-physical operators which has been discussed in part in [11].
As far as the mixing among EOM-vanishing operators is concerned, we have calculated only the first two rows of the corresponding matrixẐ where the question marks stand for entries that we have not calculated.
In the case of the mixing of evanescent into other operators, we have calculated only the first four rows of the corresponding matricesẐ Here once again question marks denote entries that we have not computed. Clearly, the mixing of evanescent into other operators does not affect the O(α 3 s ) mixing of physical operators at all, and thus is given here only for completeness. At the three-loop level we have calculated only a small subset of entries ofẐ (3,1) which are summarized below. AgainẐ (3,2) andẐ (3,3) can in principle be obtained using Eqs. (25). The single poles we have calculated read