General non-leptonic $\Delta F=1$ WET at the NLO in QCD

We reconsider the complete set of four-quark operators in the Weak Effective Theory (WET) for non-leptonic $\Delta F=1$ decays that govern $s\to d$ and $b\to d, s$ transitions in the Standard Model (SM) and beyond, at the Next-to-Leading Order (NLO) in QCD. We discuss cases with different numbers $N_f$ of active flavours, intermediate threshold corrections, as well as the issue of transformations between operator bases beyond leading order to facilitate the matching to high-energy completions or the Standard Model Effective Field Theory (SMEFT) at the electroweak scale. As a first step towards a SMEFT NLO analysis of $K\to\pi\pi$ and non-leptonic $B$-meson decays, we calculate the relevant WET Wilson coefficients including two-loop contributions to their renormalization group running, and express them in terms of the Wilson coefficients in a particular operator basis for which the one-loop matching to SMEFT is already known.


Introduction
Non-leptonic ∆F = 1 decays mediated by down-type b → d, s and s → d transitions offer some of the most important probes of the quark-mixing mechanism in the Standard Model (SM). In particular, they play an important role in the Cabibbo-Kobayashi-Maskawa (CKM) [1,2] matrix determination, and enable a detailed understanding of CP violation in the quark sector. The most prominent examples are the "gold-plated" decays B d → J ψ K ( * ) and B s → J ψ φ that are driven by the b → scc tree-level transitions. Interference of their amplitudes with time-dependent ∆F = 2 meson mixing amplitudes gives the most precise determination of the CKM angles known as β d and β s , respectively. Another example are the measurements of the CKM angle γ in B → DK decays mediated by tree-level b → cud transitions. Probably one of the most sensitive observables for CP-violating effects beyond the SM in quark-flavour physics is the ratio ε ′ ε in K → ππ decays where a large hierarchy among magnitudes of the CKM elements makes the SM contribution strongly suppressed. A recent review of these topics and the present status of non-leptonic meson decays can be found in Ref. [3]. See also Refs. [4][5][6][7].
The main obstacle to fully profiting from the potential of non-leptonic decays in tests of the SM and searches for New Physics (NP) is our imprecise knowledge of nonperturbative hadronic bound-state effects. In the SM, the left-handed nature of weak interactions leads to a specific form of the corresponding Weak Effective Theory (WET) 1 that allows us to develop specially tailored strategies for elimination of uncertainties that stem from such bound-state effects in the aforementioned non-leptonic decays. However, contributions beyond the SM (BSM) can upset such strategies. Firstly, they may introduce additional CP-violating phases. Secondly, they may give rise to additional operators in the WET, thereby increasing the number of unknown observable-specific hadronic matrix elements of operators. The latter are governed by strong dynamics, which requires first-principle nonperturbative methods for their evaluation.
The available accuracy of first-principle methods in calculations of hadronic matrix elements is strongly process-dependent. Several new lattice results for K → ππ became available over the last five years [8][9][10]. Further insight can be gained from large-N c QCD calculations [11,12]. A good description of the ∆I = 1 2 rule [10,11] and improved predictions for ε ′ ε have been obtained [10,13]. Two-body decays of B mesons can be treated with QCD-factorization methods [14][15][16] that are based on an expansion in inverse powers of the heavy-quark mass m b . The leading-power contributions depend on a few universal hadronic matrix elements. They usually provide a good qualitative description of decay rates. However, predictions for CP-violating observables often suffer from unknown subleading-power corrections, for which a proper factorization has not yet been established. In consequence, precision analyses can hardly be performed. Certain observables in B physics can be treated with the Heavy Quark Expansion (HQE), as, for example, the decay-width difference ∆Γ d,s of neutral B mesons [4,17]. The corresponding hadronic matrix elements are accessible to lattice methods [18], and high precision is expected in the future [19,20]. On the other hand, a rigorous description of bound-state effects in non-leptonic charm decays is currently not available [21].
A necessary prerequisite is the knowledge of the WET Wilson coefficients (effective couplings) at the relevant low-energy scale µ had ≲ 5 GeV at which the non-leptonic processes of interest take place. They must be determined with the help of Renormalization-Group (RG) improved perturbation theory. The first necessary step is a determination of the initial Wilson coefficients through the matching to a UV completion at the electroweak scale µ ew ∼ 100 GeV. Secondly, Anomalous Dimension Matrices (ADMs) have to be calculated. They originate from QCD and QED interactions associated with the residual SU(3) c ⊗U(1) em gauge symmetry. The RG evolution of the Wilson coefficients from µ ew to µ had is governed by the ADMs. When crossing the quark-decoupling thresholds, the corresponding corrections need to be included in the Next-to-Leading Order (NLO) RG evolution. This concerns, in particular, the ∆S = 1 decays for which subsequent decoupling of the bottom and the charm quarks implies working with N f = 4, 3 active flavours. 2 Eventually, in the last step, the aforementioned hadronic matrix elements appearing in specific observables need to be evaluated with nonperturbative methods.
Except for the latter step, calculations of the QCD and QED effects can systematically be performed within perturbation theory, enabling a precise calculation of the Wilson coefficients at µ had . Cancellations of intermediate renormalization-scheme dependences are well understood. The remaining final scheme dependence gets cancelled in predictions for observables by the corresponding scheme dependence of hadronic matrix elements. In fact, based on the universal character of the WET, 3 specific observables can be predicted adapting the most general WET operator basis, while keeping the initial matching conditions of the Wilson coefficients at µ ew arbitrary. These so-called master formulae for various observables can then be used in any UV completion by simply substituting the corresponding initial matching conditions. All radiative corrections below µ ew are properly included, and cancellation of scheme dependences is guaranteed, provided the initial matching conditions are determined in the same scheme as the ADMs. 4 The perturbative part of this program has been carried out for the operator basis generated by the SM to very high orders: the initial Wilson coefficients are known up to the Next-to-Next-to-Leading Order (NNLO) in QCD [24], with partial NLO electroweak (EW) corrections [25]. The ADMs are also calculated to the required precision, including the full NLO QCD and QED [26][27][28], and partial NNLO QCD contributions [29]. Further, also the NLO QCD and QED threshold corrections for N f = 5, 4 are available [30], as well as the NNLO QCD ones [31,32]. As far as the observable-specific hadronic matrix elements are concerned, efforts towards their determination have mainly been focused at the ones that arise in the SM.
On the other hand, no systematic study exists to similar precision for the most general ∆F = 1 WET basis of non-leptonic operators beyond the SM. First steps towards the inclusion of NLO perturbative corrections for non-leptonic ∆F = 1, 2 processes in the presence of new operators have been taken in the literature. The ADMs for the complete set of nonleptonic ∆F = 1, 2 processes were calculated at the NLO in QCD in Ref. [33], extending an earlier analysis of Ref. [34]. The specific choice [33] of the operator basis, called here BMU basis, reduced the number of operators that mix into the QCD-and QED-penguin operators as much as possible, thereby simplifying the structure of the RG evolution for N f = 3, 4, 5. The BMU basis contains the conventional SM basis [26][27][28] as a subset. 2 Currently most lattice calculations of K → ππ matrix elements are done with N f = 3 active quark flavours. In the future, treating the charm quark as dynamical should allow to raise µ had sufficiently above 1 GeV, into a more perturbative regime, and work with N f = 4 only. 3 As long as no BSM light degrees of freedom matter for the considered observables. 4 Examples of master formulae can be found in the case of SMEFT and WET for ∆F =2 observables in Ref. [22], and for g − 2 in Ref. [23].
In recent years, the idea of a considerable mass gap between the scale of new physics Λ and the electroweak scale µ ew ≪ Λ gave rise to broader interest in the Standard Model Effective Field Theory (SMEFT) [35,36]. It assumes the SM field content and the SM gauge symmetry group SU(3) c ⊗ SU(2) L ⊗ U(1) Y . In this context, another WET operator basis has been introduced [37] for the tree-level matching of SMEFT onto WET, called throughout the JMS basis. 5 This basis has recently been used in determination of the complete one-loop matching 6 of SMEFT onto WET [44]. Actual renormalizable UV completions are understood to be matched onto a stack of effective theories, the last two being subsequently the SMEFT and the WET.
In this work, we provide a connection of the JMS basis to the BMU basis beyond the leading order, which requires a correct treatment of scheme dependences. This combines the advantages of the JMS basis in the matching to SMEFT with the simplified RG evolution and known matrix elements in the BMU basis. In Section 2 and Appendix A, we recall definitions of the non-leptonic operators in the JMS and BMU bases for N f = 5. A transformation of the Wilson coefficients beyond Leading Order (LO) from the JMS to the BMU basis is provided in Section 3. Solutions of the RG equations in the BMU basis together with a discussion of the N f = 4, 3 cases are given in Section 4. In Section 5, we elucidate how the results presented here should be incorporated into a general SMEFT NLO analysis of non-leptonic decays. In particular, we indicate how various non-physical scheme dependences cancel between each other. In Section 6, we provide a numerical analysis for N f = 5 that is relevant for nonleptonic B decays, as well as for N f = 3 that matters for the BSM master formula for ε ′ ε in K → ππ at the NLO in QCD. A summary of our results is given in Section 7. The relevant ADMs and threshold correction matrices are collected in the Appendices B and C, respectively.

Non-leptonic operators
We consider in this work non-leptonic ∆F = 1 transitions of down-type quarks d i → d j with i ≠ j and i, j ∈ {d, s, b}. We focus on ∆C = 0 processes that proceed through transitions with only the down-type quarks or both the down-and up-type quarks: If we did not require ∆C = 0, another type of transition would need to be considered in addition, namely A generalization to ∆F = 1 up-type transitions c → u is straightforward. Below the electroweak scale µ ew ∼ 100 GeV, these processes are described within the WET. It allows us to systematically parameterize effects beyond the SM by including all possible operators consistent with the local gauge symmetry. Flavour-changing interactions in the WET Lagrangian that matter for our purpose are contained in where the sum runs over all the relevant four-quark operators Q i . We are going to restrict our attention to four-quark operators, leaving aside the chromo-and electro-magnetic dipole ones. The corresponding Wilson coefficients C i depend on the renormalization scale µ. They are obtained at the electroweak scale µ = µ ew in the matching from either the SMEFT or some of its possible extensions. As far as the RG evolution to lower scales is concerned, we prefer to perform it in the BMU basis, see Appendix A.2. In this basis, the operator mixing is considerably simplified. Moreover, the ADMs that govern the RG equations are known already at the NLO in QCD for the complete set of BSM operators in the BMU basis. Last but not least, calculations of the NLO corrections to various observables in the literature have been done in this basis. A transformation from the JMS to the BMU basis at the NLO in QCD is given in Section 3. The non-leptonic operators under consideration are listed in Appendix A. They can be grouped into six sectors according to their Dirac structures. The sector names and the corresponding numbers of operators in each sector are: VLL(8), VLR (16), SRR (16), VRR(8), VRL (16), SLL (16). Operators from the latter three sectors (called chirality-flipped ) are pairwise related to the former ("non-flipped") ones via interchange of the chirality projectors P L ↔ P R , where P L,R ≡ (1 ∓ γ 5 ) 2. The overall number of operators is 2 × 40 = 80 for N f = 5 and ∆C = 0. Since QCD and QED conserve parity, the full 80 × 80 ADM exhibits a 2 × 2 super-block structure, where each of the super-blocks has dimensions 40 × 40. The same refers to the matrices that describe the JMS→BMU basis transformation.
We choose the following ordering of the "non-flipped" operators in the subsequent sectors of the JMS basis (Appendix A.1) 7 In the chirality-flipped sectors, the ordering takes then the following explicit form: Obviously, the BMU basis (Appendix A.2) consists of the same number of operators. However, the VLL, VLR and the corresponding chirality-flipped sectors get subdivided into sets that do not mix under the QCD and QED renormalization. For the purpose of our transformation of the Wilson coefficients between the JMS and BMU bases in Section 3, we introduce the following reference order: SRR ∶ Q 25 , . . . , Q 40 .
In the chirality-flipped sectors (VRR, VRL and SLL), the ordering is analogous, up to shifting the operator subscripts (Q i → Q 40+i ).

Basis change JMS → BMU
The WET Wilson coefficient RG evolution is most conveniently performed in the BMU basis for which the corresponding NLO QCD ADMs are available from Ref. [33]. Many observables have already been calculated employing this very basis. Since the one-loop SMEFT to WET matching has been performed in the JMS basis, the obtained Wilson coefficients need to be transformed. At the tree level, it is sufficient to find linear relations between operators of the two bases in D = 4 spacetime dimensions. In our case, they are given by an 80 × 80 matrix R: Beyond the tree level, the basis choice is a part of the MS renormalization scheme definition, as the treatment of the Dirac algebra in D ≠ 4 dimensions does matter. In consequence, the Wilson coefficient transformation receives perturbative NLO corrections of the form 7 The Wilson coefficients themselves are determined starting with tree level, and subsequently adding higher orders in the running couplings α s (µ) and α em (µ), namely A complete evaluation of the tree-level contributions C (00) i in the matching of SMEFT onto WET is known in the JMS basis from Ref. [37]. It has been extended to the one-loop level in the SM gauge and Yukawa couplings in Ref. [44]. The latter calculation thus provides the NLO QCD correction C , too. In consequence, these corrections depend, in general, on the heavy boson and/or the top-quark masses. In the absence of tree-level contributions in certain cases, they constitute the first non-vanishing contributions, and higher orders in the perturbative expansion may become relevant for phenomenological purposes. We caution the reader that the one-loop matching results in Ref. [44] are not obtained for the mass-eigenstate basis -see details in section 6 of that paper. Hence, they require an additional rotation of the quark fields to the masseigenstate basis.
The matrices ∆r i in Eq. (15) can be systematically determined in perturbation theory. They depend on definitions of the physical and evanescent operators in both bases. In practice, their determination may require introducing supplemental evanescent operators to describe the basis transformation [45].
One of the complications beyond LO arises due to the Fierz relations required for the transformation between the JMS and BMU bases. Since they cannot be simply extended to D = 4 − 2 in the dimensional regularization beyond LO, Fierz-evanescent operators have to be introduced at the NLO in QCD. The adapted definitions for these operators in Refs. [33] and [44] are equivalent, after accounting for their different colour structures.
Apart from the Fierz-evanescent operators, other evanescent operators are encountered in Eq. (4.7) of Ref. [44]. They come with general O( ) terms that depend on arbitrary constants, of which only a ev , b ev and c ev are relevant here. These constants were kept arbitrary in the calculation of the SMEFT matching conditions for the WET Wilson coefficients in the JMS basis [44]. On the other hand, the corresponding evanescent operators in the BMU basis [33] are defined with the explicit choice a ev = b ev = c ev = 1, which coincides with the so-called Greek projection method [46,47]. In principle one could choose to calculate the NLO transformation matrices ∆r s and ∆r em keeping a ev , b ev and c ev arbitrary. Then the dependence on these constants would cancel between the Wilson coefficients in the JMS basis and the transformation matrices. However, we refrain from such a complication, and set a ev = b ev = c ev = 1 from the outset in our calculation of the transformation matrices. Therefore, also the matched Wilson coefficients in the JMS basis must be used with this particular choice.
Explicit examples of bases transformations involving Fierz-evanescent operators at the NLO in QCD are given in Ref. [33]. Very similar transformations are encountered here. Therefore, we shall not describe the calculation in detail but rather present below our final results only.
The 80 × 80 tree-level transformation matrixR in Eq. (15) readŝ with the following blocks: The factors N c arise from the identity δ αδ δ βγ = δ αβ δ γδ N c + 2 T A αβ T A γδ for the generators T A of the SU(3) c colour algebra. As one can see in the above equations,R SRR takes a simple block-diagonal form that is easy to determine from Fierz identities that relate operators of the same flavour content. On the other hand, the matricesR VLL andR VLR are somewhat more complicated, which is due to different guiding principles in organizing the JMS and BMU bases. The JMS one is organized according to the up-and down-type quark content of the quark currents, while the BMU one is organized according to the particular pattern of the QCD and QED mixing via penguin diagrams. In consequence, the QCD and QED penguin operators in the BMU basis contain (weighted) sums over the up-and down-type quark currents.
The NLO QCD correction ∆r s in Eq. (15) equals to where the 40 × 40 matrix ∆X has the following block-diagonal structure: Its upper-left 24 × 24 block corresponds to the VLL and VLR operators. Inside this block, the upper-left 8 × 18 sub-block ∆Â turns out to be a product of 8 × 1 and 1 × 18 matrices, namely The remaining lower-right 16 × 16 block of ∆X corresponds to the SRR sector. All its nonvanishing entries appear in the diagonal sub-blocks ∆B and ∆Ĉ whose explicit forms are found to be As we already have mentioned, non-vanishing entries in ∆r s arise whenever the tree-level transformation between the JMS and BMU bases involves Dirac algebra manipulations that are valid in D = 4 dimensions only. In such cases, evanescent operators need to be taken into account when working out the basis transformation in D ≠ 4 dimensions. The matrix ∆Â corresponds to the VLL and VLR sectors taken together. The transformations described by ∆B and ∆Ĉ correspond to the SRR-sector operators with two or three different flavours, respectively. To determine these matrices, we had to evaluate one-loop diagrams with insertions of various Fierz-evanescent operators. In the SRR sector case, no penguin diagrams can give non-vanishing contributions, and only the current-current diagrams matter. Conversely, no current-current diagrams with Fierz-evanescent operator insertions turn out to contribute to ∆Â once our conventions have been specified as described in the beginning of this section (a ev = b ev = c ev = 1). The simple structure of ∆Â in Eq. (22) stems from the fact that only penguin diagrams with Fierz-evanescent operator insertions contribute to this matrix.

Renormalization group equations
The Wilson coefficients satisfy the RG equation governed by the ADMγ(µ) of the WET operators. They are renormalized in the MS scheme, implying the proper treatment 8 of contributions from evanescent operators [47][48][49]. The ADM has the following perturbative expansion: where the ellipses stand for higher-order terms, i.e. O(α 3 s , α em α 2 s , α 2 em ). The ADM has a simple block structure in the BMU basis [33] because many sub-sectors of non-leptonic operators are closed under QCD and QED renormalization. Operators in such sub-sectors mix only among themselves, which we refer to as "self-mixing". The only exceptions are the sub-sectors VLL,u (or, alternatively, VLL,c), VLL,i + j, VLR,i + j and their chirality flipped counterparts where mixing into the QCD-and QED-penguin operators occurs in addition. In consequence, off-diagonal blocks arise in the corresponding parts of the ADM. Focusing on the chirality non-flipped case, one can writê where the diagonal blocks other thanγ P arise from self-mixing only, while the non-vanishing off-diagonal blocks are due to mixing into the QCD-and QED-penguin operators.
The one-loop ADMsγ (10) for QCD have been collected in Refs. [12,33,50] for the various sectors of operators. The two-loop QCD ADMγ (20) is known for the SM operators from Refs. [26,28] and for the BSM operators from Ref. [33]. 9 Higher order QCD corrections are numerically particularly important due to the large value of the strong coupling α s , especially when it evolves to lower scales µ ≲ m b . Moreover, at the NLO, non-trivial scheme dependences cancel for the first time.
On the other hand, QED corrections are negligible in general, with exceptions when the leading QCD contributions cancel in a specific observable. For what concerns the BSM subsectors (SRR,Q, SRR,i(j), SRL,Q, . . . ), the RG effects due to QED are generically small thanks to the block structure of the full ADM. Such a structure guarantees that they are not generated via operator mixing from other sectors under both the QCD and QED RG evolution. Therefore, in the self-mixing of each BSM sector, 10 the dominant effect comes from the LO and NLO QCD ADMsγ (10) andγ (20) . In this respect, the role of QED-penguin operators P QED , see Eq. (67), is special because they are the only ones that can be generated due to operator mixing at one-loop in QED from the sub-sectors VLL,u, VLL, i + j and VLR,i + j, as shown in Eq. (27). Still, whether the QED contributionsγ (01) andγ (11) to the ADM are really required depends on the observable and on the BSM scenario under consideration. We list therefore these additional matrices only for those sectors that mix into the QED-penguin operators as in Eq. (27).
In the SM, the QED-penguin operators receive direct matching contributions at µ ew that are one-loop and α em -suppressed. Therefore, the leading contributions to their Wilson coefficients at low scales µ had come actually from the RG mixing due to VLL,u operators that is given byγ (01) . Such contributions are enhanced by ln(µ ew µ had ) with respect to the direct matching contributions. On the other hand, in BSM scenarios that allow treelevel contributions to the penguin operators, the direct matching contributions may become relevant, even though they are usually suppressed by the corresponding NP scale Λ. In such cases, the numerically leading RG contribution may likely be due to the self-mixing of QED-penguin operators under the one-loop QCD ADMγ (10) . Whether the SM or the BSM contributions are numerically leading in the QED-penguin operators depends then strongly on the numerical values of the NP couplings w.r.t. α em , and has to be checked case by case. We emphasize that this is only of real concern for observables where the contributions of QED-penguin operators can compete numerically with the ones of charged-current or QCDpenguin operators, as, for example, in ε ′ ε. In view of that the SM matching contributions for QED-penguin operators are known immediately in the BMU basis, 11 the knowledge of the matrix ∆r em in Eq. (15) is only needed for BSM scenarios where the NP contributions arise as loop-suppressed QED corrections to guarantee full scheme independence. Since the RG equations are linear, the SM and BSM contributions can always be evolved separately from each other.
The solution to the RG equation (25) with several threshold crossings from N f = 5 to N f = 3 active quark flavours has the following general form: Here the scales are ordered as µ 4 < µ 5 < µ ew , and the threshold scales are comparable to the masses of the decoupled quarks: are vectors with dimensions that correspond to the number of operators present in the effective theories. The matricesM (N f ) account for threshold corrections that arise when going from the N f -to the (N f − 1)-flavour WET, namely They do not need to be quadratic matrices if the number of operators changes at the threshold. They are found by calculating partonic matrix elements of operators in the N f -flavour WET. In evaluating such matrix elements, we treat the N f -th quark as massive with mass m N f , and all other quarks as massless. Scheme dependences cancel in the product with the Wilson coefficients. Finally, the RG running of the Wilson coefficients is accounted for by the evolution operatorsÛ that multiply the initial conditions ⃗ where the ellipses indicate higher orders in α s and α em . The matrix V diagonalizes the transposed LO ADMγ (10) . To describe the matrices D and H in Eq. (31), we begin with defining where G (0) is diagonal, while G (1) usually contains off-diagonal terms. Then D is given by the diagonal matrix while the expression for H reads Here, β 0 and β 1 are the first and second coefficients of the beta function for α s , respectively. For processes/observables where the QED-penguin operators cannot be neglected, the above expressions need to be supplemented with non-leading QED corrections. Analytic solutions to the RG equations in such cases can be found in Refs. [30,51]. The LO and NLO ADMs for QCD are collected in Appendix B for arbitrary N f . As most of the sub-sectors of operators have only self-mixing, 12 the RG evolution in such cases can be performed for each sub-sector individually, except for the ones indicated in Eq. (27). In Appendix B we also provide the LO QED and the NLO QED×QCD ADMs, i.e.γ (01) and γ (11) , respectively.
In the remainder of this section, we discuss linear relations that arise among the BMU operators when going from N f = 5 to N f = 4, and further to N f = 3. In practice, they are only relevant for Kaon physics where the convention is usually d i = d and d j = s. 13 These relations reduce the number of linearly independent operators in the corresponding WETs. Some of the remaining operator Wilson coefficients are then affected by threshold corrections that arise already at the tree level. Such effects were taken into account in Refs. [52,53] where the LO BSM master formula for ε ′ ε was derived.
When the considered linear relations involve Fierz identities, a careful treatment beyond the LO is necessary, with inclusion of the proper evanescent operators. It turns out that all such cases involve the QCD-and QED-penguin operators. In this respect, it is a common practice in the literature to work for N f = 3, 4 with a redundant operator set in the framework of the SM [30], retaining the analogues of all the QCD-and QED-penguin operators that were present for N f = 5. In the next sections, we will comment on the consequences of such an approach for the ADMs and threshold corrections.
At the b-quark decoupling scale µ 5 , we remove operators that mediate d i → d j bb. There are six operators in the two sub-sectors SRR,b (4) and SRL,b (2) that, in principle, can give rise to threshold corrections to the QCD-and QED-penguin operators. Furthermore, in the QCD-and QED-penguin operators Q 3,...,10 , the sum runs now only over q = u, d, s, c. In consequence, some of our operators become linearly dependent: where F denotes a Fierz-evanescent operator that requires a careful treatment beyond the LO [30]. In the SM, as mentioned before, it is a common practice to retain Q 10 and work with a redundant set of operators. The relations in the VLL,i + j and VLR,i + j sub-sectors are valid also in D ≠ 4 dimensions, hence without any complication due to evanescent operators. Furthermore, the elimination of the VLL,i + j and VLR,i + j sub-sectors removes them as sources of mixing into the QCD-and QED-penguin operators, see Section 4. The operators VLL,i + j and VLR,i + j simply give tree-level threshold corrections to the involved QCDand QED-penguin operators. This gives a basis for N f = 4 with 40 − 10 = 30 operators in the "non-flipped" sector, and analogously in the chirality-flipped one. There are 31 operators in the redundant set though when Q 10 is retained. At the c-quark decoupling scale µ 4 , we remove operators that mediate d i → d j cc. There are six operators in the two sub-sectors SRR,c (4) and SRL,c (2). Furthermore, in the QCDand QED-penguin operators, the sum runs now only over q = u, d, s. In effect, two of the QCD-and QED-penguin operators become linearly dependent: In the SM, it is common to retain both operators and work with a redundant set. We note that F in Eqs. (35), (38) and (39) Both operators give rise to tree-level threshold corrections to the QCD-and QED-penguin operators. This gives additional 4 + 4 + 2 operators to be discarded, such that for N f = 3 the operator basis contains 30 − 10 = 20 operators in the "non-flipped" sector, and analogously in the chirality-flipped one. There are 23 operators though when Q 4,9,10 are retained.

Comments on scheme cancellations
It is well known that the two-loop anomalous dimensions depend on the renormalization scheme (RS) for the effective theory fields and parameters, including the Wilson coefficients. There is a correlation between RS dependences of the Wilson coefficients and matrix elements of the corresponding operators. In dimensional regularization, the MS scheme definition includes specifying the necessary evanescent operators, as well as choosing particular Dirac structures in the physical operator basis [3]. The reader is already familiar with these issues after studying our paper up to this point. The final physical amplitudes are RS-independent. In particular, they do not depend on definitions of evanescent operators. Such dependences are cancelled between ⟨ ⃗ Q BMU (µ had )⟩ and ⃗ C BMU (µ had ). The latter Wilson coefficients can be expressed in terms of the SMEFT ones at the new physics scale Λ as follows: Here, the matrixM JMS summarizes the JMS→BMU basis transformation in the WET as given in Eq. (15). The matrixK summarizes matching relations between the WET in the JMS basis and the SMEFT in the Warsaw basis [36]. It has been calculated including one-loop contributions in Ref. [44]. Explicitly Various contributions to Eq. (42) must be evaluated using the same RS to guarantee proper cancellation of the RS dependences. To describe it in more detail, let us factorize out the NLO contributions to the QCD evolution matrices 14 whereÛ (0) i are the RS-independent LO evolution matrices. On the other hand,Ĵ i stem from the RS-dependent two-loop ADMs, which makes them sensitive to the evanescent operator definitions. Explicit general expressions forÛ (0) i andĴ i can be found in Ref. [3]. Perturbative expansions ofM JMS (µ ew ),K(µ ew ) and ⃗ C SMEFT (Λ) take the form and Having all the above expressions at hand, we can easily trace out cancellations of the RS dependences: • The RS dependence of ⃗ C SMEFT is cancelled by the one inĴ SMEFT entering the last factor on the r.h.s. of Eq. (45).
• The one inĴ BMU entering the first factor on the r.h.s. of Eq. (44) is cancelled by the one in ⟨ ⃗ O BMU (µ had )⟩.
• Finally, the remaining RS dependences, including those related to the evanescent operator definitions cancel in the product In order to obtain RS-independent results, one has to ensure that evaluation of the matching matricesM JMS andK is consistent with the scheme dependences inĴ BMU andĴ SMEFT . This indeed has been done forM JMS in the present paper. As far asK in Ref. [44] is concerned, using precisely the same RS is recommended in the future calculation ofĴ SMEFT . Further details on RS dependences can be found in Refs. [29,33,45].

Numerical analysis
In this section, we study the impact of NLO QCD corrections on the four-quark operators in the following two cases: • non-leptonic B decays governed by b → d and b → s transitions, • the decays K → ππ governed by s → d transitions.
The central quantities are the RG evolution operators connecting the WET Wilson coefficients in the BMU basis at low-energy scales µ had to the ones of the JMS basis at the electroweak scale µ ew The basis transformation from the JMS to the BMU basisM JMS has been outlined in Section 3. The exact form of the evolution operatorÛ BMU (µ had , µ ew ) depends on the number of crossed quark thresholds. The most general form required for our purposes is given in Eq. (28). We study here the elements that depend explicitly on the two scales µ ew and µ had . In the s → d case, there is also an implicit dependence on the quark threshold scales µ 5 and µ 4 that cancels out up to residual higher-order effects. Apart from that, u ab (µ had , µ ew ) depend only on the gauge coupling constant initial values at µ = m Z = 91.1876(21) GeV [54] α (5) s (m Z ) = 0.1179(10), α as well as on the solutions to the RG equations for these couplings, mainly α s . Throughout our analysis, the RG evolution of α s is performed at the 3-loop level, which means that partial beyond-NLO effects are included, too. Furthermore, the MS quark mass values [54] m c (m c ) = 1.27(2) GeV, are used in the threshold corrections collected in Appendix C. We note that the renormalization scheme of the quark masses plays no role at the NLO level in QCD. The scales for the threshold crossings are set to The quantities u ab (µ had , µ ew ) are found by solving the general RG equation (25) for the Wilson coefficients. We consider two options: 1. A perturbative analytic solution at the NLO in QCD, as outlined in Eqs. (31)- (34). See also Refs. [30] and [51] for the coupled QCD×QED case.

2.
A direct numerical solution to the system of ordinary differential equations (25) together with the RG equation for the couplings α s and α em .
Both options are equivalent up to higher order corrections but numerically they differ by residual effects. Furthermore, numerical ambiguities of similar origin can arise from the treatment of products of quantities that have perturbative expansions in α s and α em . For example, products like (28) and (49) can either be expanded, thereby discarding higherorder contributions beyond the NLO, or kept unexpanded, thereby including partial NNLO contributions. Both results are equivalent at the NLO level. Numerical differences between them can be used to estimate uncertainties that arise due to lacking NNLO contributions. Our default choice is the numerical solution to the RG equation, unless stated otherwise. Products of the RG evolution operatorsÛ (N f ) and threshold correctionsM (N f ) are kept unexpanded.

Non-leptonic B decays
Hadronic B decays are governed by the four-quark operators specified in Eq. (1) and (2), with the flavour indices i = b and j = d or j = s. In the absence of quark flavour thresholds, the quantities of interest are The index b runs over all the 2 × 40 operators of the JMS basis, whereas a corresponds to all the 2 × 40 operators of the BMU basis in the N f = 5 theory, as described in Section 2 and Appendix A. The matrices u ab (µ had , µ ew ) consist of two 40 × 40 identical diagonal blocks corresponding the non-flipped and chirality-flipped operators. Many entries vanish due to the block structure of the matrixM JMS , as well as of the ADMs. As examples, we inspect closer the SRL, u, the SRR, i and the SRR,u sub-sectors given in Eqs. (72), (74) and (76), respectively. They consist of the following operators: in the JMS and BMU bases, respectively. We note that the results for the sectors SRL, Q and SRR, Q with Q = c, d k≠i,j are identical to the results of the SRL, u and SRR, u sectors, respectively. Furthermore, the results for the SRR, i sector hold for the SRR, j sector, too. Our scale choices in the numerical examples are µ ew = 160 GeV and µ had = 4.2 GeV. First, we investigate the impact of the NLO corrections on the transformation matrix M JMS between the JMS and BMU bases, as defined in Eq. (15). We find at µ ew numericallŷ at the LO and NLO in QCD, respectively. No evanescent operators are involved in the transformation of the SRL,u sector, which makes the corresponding NLO corrections vanish. 15 The NLO effects in the SRR, i sector amount to a relative change of +3% in the (1,2) element and −4% in the (2,2) element. In the SRR, u sector, the NLO corrections give rise to two non-vanishing entries (2,3) and (4,3). The largest modifications are in the elements (1,4) and (2,4), amounting to −16% and +4%, respectively. The corresponding u ab (µ had , µ ew ) are the ones with for the SRL, u and SRR, i sectors, respectively.. 16 Their µ-dependence is shown in Figure 1.  difference between the numerical and analytical solutions is larger at the LO, and becomes reduced at the NLO, which indicates reduction of renormalization-scheme dependences. The numerical effects we observe correspond to the particular renormalization scheme we have chosen for the Wilson coefficients and the running couplings. In particular, we always use the three-loop RG equations for the running of α s , irrespectively of whether the RG evolution of the Wilson coefficients is performed at the LO or NLO. Otherwise the shifts from the LO to the NLO would be much larger but unrelated to the corrections we have calculated. Let us note that differences between the numerical and analytical solutions of the NLO RG equations for the Wilson coefficients are formally of higher order but should not be misinterpreted as estimates of the overall higher-order effects.
The plots remind us impressively that the RG evolution is an important effect which leads to sizable changes of various Wilson coefficients when going from µ ew to µ had . The impact of NLO corrections for the numerical solution is summarized in Table 1 for µ had . We observe corrections in the ballpark of 10% to 25% w.r.t. the LO result for our choice of the renormalization scheme. The complete list of u ab (µ had , µ ew ) at the NLO for our choice of µ ew and µ had is given in Table 2. For comparison, we provide the LO results in Table 2 in parenthesis. One can see there that the NLO corrections in the SRR, Q sectors are even larger, being in the ballpark of 10% to 70%, with an extreme case of > 200% for one element. The size of NLO corrections in the VLL×VLL and VLL×VLR blocks range from a few percent up to > 100% for some entries. Large NLO corrections are observed for u ab entries whose absolute values at the LO are small. In the VLR×VLR block, the NLO corrections do not exceed 40%.
We note that inclusion of the QED and QED×QCD ADMs in the RG evolution modifies the existing entries in Table 2 at the level of 1% or less, with a few exceptions up to 3%. Furthermore, there are additional non-vanishing entries for the QED penguin operators Q 7,8,9,10 . They are all very small compared to the ones generated by the QCD mixing. As mentioned before, their numerical relevance strongly depends on both the observable and the BSM scenario under consideration.      Table 2: The elements u ab (µ had , µ ew ) of the NLO (LO) RG evolution operator in Eq. (53) for N f = 5. They connect the JMS Wilson coefficients at µ ew = 160 GeV to the BMU ones at µ had = µ 5 = 4.2 GeV.    Table 3: The elements u ab (µ had , µ ew ) of the NLO (LO) RG evolution operator in Eq. (59) for N f = 3 at µ had . They connect the JMS Wilson coefficients at µ ew = 160 GeV to the BMU ones at µ had = µ 4 = 1.3 GeV.

K → ππ and ε ′ ε
Non-leptonic Kaon decays are especially sensitive to CP-violating effects beyond the SM because of the strong suppression of the SM contribution. In particular, the observables ε and ε ′ measure indirect and direct CP violation in K 0 −K 0 mixing and K 0 decay into ππ, respectively. The world average from NA48 [55] and KTeV [56,57] collaborations for their ratio ε ′ ε reads The K → ππ hadronic matrix elements of the SM operators contributing to the ratio ε ′ ε from RBC-UKQCD Lattice QCD collaboration [10] together with refinements in the estimate of isospin-breaking effects [58,59] and QCD corrections [31,32,60] indicate that although the SM is compatible with the experimental world average, still significant room for NP is left in this ratio. The most recent analysis summarizing the present status of ε ′ ε in the SM implies [13,59] (ε ′ ε) SM = (13.9 ± 5.2) × 10 −4 .
Thus, modifications of ε ′ ε by NP as large as 10 −3 are presently not excluded. In fact, as argued on the basis of the large-N c QCD in Ref. [61], the values of ε ′ ε in the SM in the ballpark of 5 × 10 −4 are still possible. The ratio ε ′ ε puts serious constraints on CP-violating phases in BSM models. Previous analyses of NP effects in ε ′ ε were restricted to the LO RG evolution of the BSM contributions. For example, the derivation of a master formula for ε ′ ε [52] or the SMEFT anatomy of ε ′ ε [53] demonstrated possible interplay of ε ′ ε with other quark-flavour and collider observables. Our current work allows us to include also the NLO QCD corrections.
Below, we consider elements of the RG evolution operator to demonstrate the numerical impact of the NLO corrections. In Kaon physics, the b-quark and c-quark thresholds are included in the RG evolution. At each threshold, the number of operators is reduced, as described in Section 2, from 2 × 40 operators in the N f = 5 theory (index b) to 2×23 operators in the N f = 3 theory (index a). 17 A complete list of u ab (µ had , µ ew ) at the NLO for our choice of µ ew and µ had is given in Table 3. For comparison, we provide the results at the LO in Table 3 in parenthesis. The NLO corrections in the N f = 3 flavour theory are, in general, larger than in the N f = 5 flavour theory, as a consequence of fast RG evolution of α s at scales reaching µ had ≈ 1 GeV. Some of the u ab entries cross zero during their RG evolution. In certain cases, we observe suppressed magnitudes at the intermediate scale µ 5 but again sizable entries at µ had . In other cases, u ab cross zero in the vicinity of µ had -see, for example, u 19,V 8 and u 26,S8 in Table 1. Close to such points, the relative impact of NLO corrections becomes amplified. Thus, the actual impact of the NLO corrections strongly depends on the considered entries u ab . However, as in the case of N f = 5, their relative impact is usually the largest for u ab with small absolute values at the LO.
The impact of our results on the CP-violating ratio ε ′ ε is discussed in a separate publication [62] by most co-authors of the present article. The previous BSM master formula [52,53] is generalized there to NLO in QCD. It provides ε ′ ε in terms of the Wilson coefficients of the JMS basis at the electroweak scale µ ew . Depending on the considered JMS Wilson coefficient, the NLO QCD corrections in the RG evolution lead to an impact on the corresponding coefficient in ε ′ ε in the ballpark of (20 − 40)%, reaching 60% in certain cases. The magnitude of such corrections is thus similar or sometimes even larger than uncertainties of nonperturbative origin in K → ππ matrix elements of the nonleptonic operators.

Conclusions
In the present paper, we presented all ingredients that are necessary to write down the amplitude for any ∆F = 1 non-leptonic transition in terms of the electroweak-scale WET Wilson coefficients in the JMS basis. Our expressions include the NLO QCD RG evolution of the Wilson coefficients from the electroweak scale down to the low energy scale at which hadronic matrix elements are being calculated. As the latter evolution is most conveniently performed in the BMU basis, an involved NLO transformation from the JMS basis to the BMU one was necessary. While rather straightforward at the tree-level, it required a very careful treatment of evanescent operators to ensure that the corresponding scheme dependence cancels the one of the two-loop ADMs that govern the NLO QCD RG evolution.
The latter transformation, that we presented in detail in Section 3, allows us to obtain expressions for physical observables in terms of the JMS Wilson coefficients at the electroweak scale. The latter coefficients are related to the SMEFT ones at the NP scale Λ via whereK is the matching matrix calculated including one-loop contributions in Ref. [44], and U SMEFT (µ ew , Λ) is the SMEFT RG evolution matrix from Λ down to µ ew . Once the SMEFT two-loop ADMs are calculated in the future, our results will allow to obtain complete NLO expressions for amplitudes of all low-energy non-leptonic decays in any BSM scenario where the relevant new particles have masses of order Λ or larger.
One should again emphasize that only a combination of the NLO QCD RG evolution in WET with the one-loop QCD corrections in the matching of SMEFT onto WET, with the NLO QCD RG evolution within the SMEFT, and the corresponding one-loop matching of a BSM scenario onto SMEFT will enable a proper cancellation of RS dependences that appear at various places in a complete analysis.
In the present paper, we have restricted ourselves to illustrating the impact of NLO QCD effects on the WET Wilson coefficients that matter for non-leptonic Kaon and Bmeson decays. As a next step, in a separate paper, an NLO-upgraded master formula for the BSM contributions to ε ′ ε is going to be presented. The upgrade will improve the one presented in Ref. [52].

A.1 JMS basis
In the JMS basis [37], the full set of dimension-six four-quark operators is divided into four classes with different chirality structures: (LL)(LL), (RR)(RR), (LL)(RR) and (LR)(LR). The non-leptonic operators that contain at least one pair of down-type quarks are given as where p, r, s, t denote quark-flavour indices in the mass-eigenstate basis. Colour structures are expressed in terms of the generators T A of SU(3) c . The listed operators include both the "non-flipped" and the chirality-flipped ones. Some of the above operators do not equal their hermitian conjugates with permuted generation indices. Although we do not list them explicitly here, they have to be treated as independent operators. Note that also Q V 1,LR uddu and Q V 8,LR uddu require such hermitian conjugates. The same holds for the operators (LR)(LR). The JMS basis contains no operators with ψσ µν ψ ′ tensor currents.
A.2 BMU basis for N f = 5 The BMU basis [33] contains several sub-sectors of operators, depending on their quarkflavour content and Lorentz structures. They form separate blocks under RG evolution, thereby minimizing the operator mixing.
The first 10 elements of the BMU basis form the SM basis that consists of two d i → d j u kūk current-current-type and eight penguin-type operators. The current-current ones where the naming used in Ref. [33] is given as well. The penguin-type diagrams (see Figure 2) with insertions of the current-current and penguin operators can give rise to mixing into the penguin ones. The eight penguin-type operators with left-handed (d j γ µ P L d i ) vector currents that are present in the SM basis can be divided into the QCD-penguin sub-sector P QCD and the P QED one The sums over q = u, d, s, c, b include all the five active quark flavours whose electromagnetic charges are denoted by Q q : Q u = +2 3 and Q d = −1 3. Note that we use here the "traditional" convention [27,28,63] for the penguin operators that is suitable up to the NLO in QCD. Beyond the NLO, the issue of traces with γ 5 makes working with an alternative basis more convenient [45,64]. Ref. [29] provides NNLO QCD transformations between the two conventions for the VLL, c and P QCD Wilson coefficients. In summary, the SM basis contains ten operators in three sub-sectors, namely VLL, u (or VLL, c),  and Only the "i + j" operators mix into the penguin operators Q 3,...10 , whereas the "i − j" ones do not. Analogous operators for d i → d j u kūk are They are numbered as ).
The scalar operators of the SRR sectors for d i → d j d i d i and d i → d j d j d j have been chosen in Ref. [33] in such a way that they contain colour-singlet currents only, analogously to the ∆F = 2 case there. Due to the Fierz identities, there are only two operators in the SRR,i sub-sector and the SRR, j one The case of d i → d j QQ with Q = u, c, d k≠i,j comprises four operators per sub-sector ). (77) Here, we use the standard definition σ µν ≡ i 2 [γ µ , γ ν ], contrary to Refs. [33,34] where no overall factor of i was present. It implies that our BMU operators with tensor currents differ by a sign from those in the original BMU paper [33], which affects a few signs in the ADMs of Appendix B.
As far as the chirality-flipped operators are concerned, their numbering in the BMU basis is given by i.e. they are found by interchanging P L ↔ P R in the "non-flipped" operators.

B Anomalous dimensions
In this appendix, we collect results of the LO and NLO QCD ADMs. We set N c = 3 wherever the N c -dependence is provided in Ref. [33], but keep the N f -dependence explicit. Moreover, we provide the one-loop and two-loop ADMs involving QED effects for those sectors that mix into the QCD-and QED-penguin operators.
The VLL,u sub-sector (or, alternatively, VLL,c) contains two operators (Q VLL,u 1 , Q VLL,u 2 ). Their one-and two-loop QCD ADMs are given in Eqs. (3.2, 3.3) of Ref. [33]. For N c = 3 and arbitrary N f they read γ (10) Since the VLL,u operators are part of the SM operator basis, their QED ADMs are known, too [27,28,63,65]γ The sub-sectors VLR,i±j and VLR,u−c contain two operators (Q ). Their one-and two-loop QCD ADMs can be derived directly from the VLR sector given in Eqs. (3.4, 3.5) of Ref. [33], see also the explanation around Eq. (4.7) in Ref. [33]. For N c = 3 and arbitrary N f one findŝ for Q i ± Q j = {d i ± d j , u − c}. For completeness, we provide here also the LO and NLO ADMs involving QED effects for Q i + Q j . They read They have been derived from Ref. [27] using the fact that only the current-current insertions are needed, and considering the special case N d = 2 and N u = 0.
The sub-sectors VLL,i ± j mediate d i → d j d idi and d i → d j d jdj that both contain a single operator Q VLL,i+j 1 and Q VLL,i−j 1 , respectively. The corresponding LO and NLO ADMs in QCD are given by the ∆F = 2 result in Eq. (2.21) of Ref. [33] γ (10) VLL,i±j = (4) ,γ (20) VLL,i±j = −7 + For completeness, we provide here also the LO and NLO ADMs involving QED effects for Q VLL,i+j 1 . They read γ (01) The 8 × 8 ADMsγ (mn) P of QCD-and QED-penguin operators is known [26][27][28]66] for arbitrary numbers of active flavours N f = N u + N d at one-and two-loops, for all the contributions mn = 10, 20, 01, 11 in Eq. (26). 18 Similarly, the mixing of current-current operators VLL,u into the penguin ones is known from the SM calculations, too. We therefore refrain from repeating these results here. The mixing of the operators (Q 11 , Q 12 , Q 13 ) from the sub-sectors VLL,i + j and VLR,i + j into the penguin operators (Q 3 , . . . , Q 10 ) has been calculated in Ref. [33] up to the NLO in QCD, givinĝ  ). Their one-and two-loop QCD ADMs are given in Eqs. (3.8, 3.9) of Ref. [33]. For N c = 3 and arbitrary N f they read Again, overall sign differences w.r.t. to Ref. [33] show up in the off-diagonal 2 × 2 blocks due to the σ µν normalization issue.
The SRL,Q sub-sectors with Q = u, c, d k≠i,j contain two operators (Q SRL,Q 1 , Q SRL,Q 2 ). Their one-and two-loop QCD ADMs are given in Eqs. (3.6, 3.7) of Ref. [33]. For N c = 3 and arbitrary N f they read (93)

C Threshold corrections
In Section 4, threshold corrections that arise when going from the N f -to the (N f −1)-flavour WET have been discussed in general terms. They are described by the matricesM (N f ) (µ f ) defined in Eq. (29). These matrices are obtained by matching matrix elements of operators of the N f -and the (N f − 1)-flavour WET, which is marked by the superscript N f . They are not quadratic matrices because the number of linearly independent operators becomes smaller when one of the quarks gets decoupled, as discussed at the end of Section 4. Their perturbative expansion readŝ where the first term arises at the tree-level, while the remaining two are due to one-loop QCD and QED corrections, respectively. For the SM-generated operators Q 1 -Q 10 and their chirality-flipped counterparts, the 10 × 10 tree-level matching matrix blocks are equal to unit matrices, 19 i.e. [M (N f ) 00 ] Q i , Q j = δ ij . This is one of the advantages of using a redundant operator set in the SM-like sector of the BMU basis for N f = 3, 4. However, as far as the BSM operators are concerned, we shall remove all the linearly dependent ones. In particular, Eqs. (36) and (37) give rise to the following non-vanishing entries in [M Analogous relations hold in the chirality-flipped sector. They are obtained by simply replacing Q i → Q 40+i . Beyond the tree level, the band c-quark decoupling proceeds through insertions of operators that contain pairs of these quarks into penguin diagrams, as shown in Figure 2. In the VLL and VLR sectors, they arise mainly from the QCD-and QED-penguin operators. They are thus known from the ancient SM calculation [30]. The 10 × 10 matricesM where The corresponding results for the QED matricesM (5) 01 andM (4) 01 are given as [30] where P i = 0, 0, 0, 0, 0, 0, 1, 0, 1, 0 i .
Among the BSM operators in the BMU basis, only the VLR,u − c operator Q 17 gives one-loop threshold corrections when decoupling the charm quark. They are given by 10 ] i,Q 17 = in the QCD and QED cases, respectively. One should remember though that Q 17 above the charm threshold does not mix into the penguin operators. Consequently, its effect on Kaon physics originates solely from the one-loop matching corrections (104) and (105), as long as the physical amplitudes are calculated in the N f = 3 WET framework. The SRR,b operators with σ µν give rise to threshold corrections to the chromo-magnetic dipole operators only, see for example Refs. [67,68]. Other SRR,b and SRL,b operators will not generate threshold corrections. The same comment applies to the corresponding SRR,c and SRL,c operators.