BSM Hadronic Matrix Elements for $\epsilon'/\epsilon$ and $K\to\pi\pi$ Decays in the Dual QCD Approach

We calculate for the first time all four-quark hadronic matrix elements of local operators possibly contributing to $K\to\pi\pi$ decays and in particular to the ratio $\epsilon'/\epsilon$ beyond the Standard Model (BSM). To this end we use the Dual QCD (DQCD) approach. In addition to 7 new mirror operators obtained from the SM ones by flipping the chirality, we count 13 BSM four-quark operators of a given chirality linearly independent of each other and of the aforesaid 14 operators for which hadronic matrix elements are already known. We present results in two bases for all these operators, one termed DQCD basis useful for the calculation of the hadronic matrix elements in the DQCD approach and the other called SD basis suited to the short distance renormalization group evolution above the 1~GeV scale. We demonstrate that the pattern of long distance evolution (meson evolution) matches the one of short distance evolution (quark-gluon evolution), a property which to our knowledge cannot be presently achieved in any other analytical framework. The highlights of our paper are chirally enhanced matrix elements of tensor-tensor and scalar-scalar BSM operators. They could thereby explain the emerging $\epsilon'/\epsilon$ anomaly which is strongly indicated within DQCD with some support from lattice QCD. On the other hand we do not expect the BSM operators to be relevant for the $\Delta I=1/2$ rule.


Introduction
The direct CP-violation in K → ππ decays, represented by the ratio ε /ε, plays a very important role in the tests of the Standard Model (SM) and more recently in the tests of its possible extensions. For recent reviews see [1,2]. In fact there are strong hints for sizable new physics (NP) contributions to ε /ε from Dual QCD approach (DQCD) [3,4] that are supported to some extent by RBC-UKQCD lattice collaboration [5,6]. Most recent SM analyses at the NLO level can be found in [7,8] and a NNLO analysis is expected to appear soon [9]. Most importantly, an improved result on ε /ε from RBC-UKQCD lattice collaboration is expected this summer.
This situation motivated several authors to look for various extensions of the SM which could bring the theory to agree with data . In most of the models the rescue comes from the modification of the Wilson coefficient of the dominant electroweak left-right (LR) penguin operator Q 8 , but also solutions through a modified contribution of the dominant QCD LR penguin operator Q 6 could be considered [15].
Here we want to emphasize that scalar-scalar and tensor-tensor four-fermion operators generated e.g. through tree-level exchanges of colour-singlet and colour-octet heavy mesons could also give significant contributions to ε /ε because they have, just like the Q 6 and Q 8 operators, chirally enhanced K → ππ matrix elements. However, whereas in the case of Q 6 and Q 8 operators significant progress in evaluating their matrix elements relevant for ε /ε by lattice QCD has been made [5,6], no lattice QCD calculations have been performed so far for these scalar-scalar and tensor-tensor operators although their two-loop anomalous dimensions have been known [31,32] for almost two decades. In fact, to our knowledge, there exist no analytic results for the matrix elements in question, even obtained by the vacuum insertion method which, in any case, as already demonstrated in several studies, totally misrepresents QCD.
We are aware of the fact that it will still take some time before lattice QCD will be able to provide K → ππ matrix elements for scalar-scalar and tensor-tensor operators. Yet, in view of the hints for NP in ε /ε, we think it is time to estimate their matrix elements in the framework of DQCD [33][34][35][36] which has been generalized in this decade [3,4,37,38] through the inclusion of vector meson contributions and improved through a better matching to short distance contributions. While not as precise as ultimate lattice QCD calculations, this approach offered over many years an insight in the lattice results and often, like was the case of the ∆I = 1/2 rule [35] and the parameterB K [36], provided results almost three decades before this was possible with lattice QCD. The agreement between results from DQCD and lattice QCD is remarkable, in particular considering the simplicity of the former approach compared to the very sophisticated and computationally demanding numerical lattice QCD one. The most recent example of this agreement was an explanation by DQCD of the pattern of values of B entering ε /ε obtained by lattice QCD [3,4] and of the pattern of lattice values for BSM parameters B i in K 0 −K 0 mixing [39]. This is also the case for hadronic matrix elements of the chromomagnetic operator presented recently in [40] that are in agreement with the result from the ETM collaboration [41].
Our paper is organized as follows. In Section 2 we recall the SM operators contributing to K → ππ decays and construct a basis of four-quark BSM operators. Consisting exclusively of 13 products of colour singlet scalar, vector and tensor bilinears, this complete basis is particularly useful for the calculations of hadronic matrix elements in the DQCD approach and will thus be called the DQCD basis. In Section 3 we recall very briefly the elements of DQCD relevant for our paper. In Section 4 we perform the evolution of BSM operators from a very low factorization scale up to scales µ = O(1 GeV), the so-called meson evolution, in the chiral limit. While this is a crude approximation, a recent analysis of BSM hadronic K 0 −K 0 matrix elements in this limit [39] was able to explain at a semi-quantitative level the pattern of the values of these matrix elements obtained by the ETM, SWME and RBC-UKQCD lattice QCD collaborations [42][43][44][45][46].
In Section 5 we present the formulae for the quark-gluon evolution in the DQCD basis while in Section 6 we demonstrate that the patterns of meson evolution of BSM operators presented in Section 4 and of the quark-gluon evolution of Section 5 are compatible with each other, assuring us that the matching of hadronic matrix elements evaluated in DQCD and of their Wilson coefficients will be satisfactory. In Section 7 we calculate the matrix elements of all BSM operators at leading order in the DQCD basis and, using the results for their meson evolution of Section 4, we obtain their values at the scale µ = O(1 GeV).
In Section 8 we introduce a different basis of 13 BSM operators which turns out to be particularly suited to the usual short distance (SD) QCD evolution. For this reason we call this basis the SD basis. We establish the relation between the DQCD and SD bases which, using the results of Section 7, allows us to obtain hadronic matrix elements of all BSM operators in the SD basis at µ = O(1 GeV). For completeness we give in Appendix D also their values in the large N limit. The 13 BSM operators of a given chirality already mentioned are all allowed by the SU(3) c × U(1) em invariance. In Section 9, following [47,48], we emphasize that in the SM effective field theory (SMEFT), based on the full SM gauge symmetry SU(3) c × SU(2) L × U(1) Y , only 7 four-quark BSM operators of a given chirality are allowed. We identify these operators in the DQCD and SD bases.
In Section 10 we calculate the K → ππ matrix elements of all BSM operators in two bases in question for values of µ to be explored one day by lattice QCD. The results in the DQCD basis demonstrate once again that the pattern of meson evolution agrees with the one of SD evolution. The ones in the SD basis can now be used in the BSM phenomenology of ε /ε and ∆I = 1/2 rule.
A brief summary of our results is given in Section 11. In a number of appendices we collect useful auxiliary material. Phenomenological implications of our results will be presented elsewhere.

Preliminaries
The isospin amplitudes A I in K → ππ decays are introduced through where the parameter h distinguishes between various normalizations of A 0,2 found in the literature. We use h = 1 but the RBC-UKQCD collaboration uses h = 3/2 implying that their amplitudes A 0,2 are by a factor 3/2 larger than ours. This difference cancels of course in all physical observables. All matrix elements listed below should be multiplied by h in case h = 1 is not used.

SM Operators
We begin by recalling the SM operators: Current-Current: QCD Penguins: Electroweak Penguins: Here, α, β denote colour indices and e q the electric quark charges reflecting the electroweak origin of Q 7 , . . . , Q 10 . Finally, (sd) V ±A ≡s α γ µ (1 ± γ 5 )d α . As we are only interested in hadronic matrix elements in this paper, the summations are only over u, d, s quarks.

Chiral Fierz Identities
The following 16 Dirac bilinears form the appropriate chiral basis for a systematic classification of all ∆S = 1 weak operators beyond the Standard model (BSM): with Within these conventions, the corresponding dual basis obeys the orthogonality property If we substitute the matrix indices by parentheses () and brackets [ ] such that each parenthesis/bracket represents a different spinor index, the completeness relation leads then to the so-called chiral Fierz identities [49] ( Anticipating the fermion field anticommutation, we can turn all the four-quark operators into products of colour singlet bilinears according to four classes of Fierz identities: Class B: Class C: Class D: Mirror operators are obtained through an obvious chirality-flip (L ↔ R).

Illustration with ∆S = 2 Operators
Tree-level neutral meson exchanges can lead to various BSM ∆S = 2 transitions. As a consequence, we have to consider the four classes of operators for a single set of Fierz-conjugate flavour indices {ab; cd} = {sd; sd}: Using this technology, we succeeded [32,39] in turning the SUSY basis [50,51] represented by 1 SM operator O 1 and 4 BSM operators O 2−5 into the SD basis of [32].

Application to the SM ∆S = 1 Operators
In the case of tree-level charged or neutral meson exchanges leading to various ∆S = 1 transitions we now have to consider four different sets of flavour indices {ab; cd}, namely for each class of operators. The last set of operators with a (ss) bilinear does not contribute to the K → ππ matrix elements at the factorization scale, the starting point of our calculation. Yet, they may contribute above this scale through the evolution into other pure ∆I = 1/2 operators. According to our generic Fierz classification and up to parity-transformations (L ↔ R), this complete basis gives rise to 5 × 4 = 20 linearly independent operators written as the product of two colour-singlet bilinears. But now, we want to build the optimal basis for a transparent evolution of the BSM operators orthogonal to the SM Q i displayed in (4)- (6). For that purpose let us first isolate the latter ones by selecting the appropriate combinations of operators in (23), with sums over q = u, d, s understood.
Class A: In fact this class covers all SM (V − A) × (V − A) operators since the charge matrix e q appearing in the electroweak weak penguins in (6) can be decomposed in the following way: and implies Class B: Here we use the fact that {ss; sd} is the mirror partner of {sd; ss}.
Class C: So we end up with a total of 7 four-quark SM operators scattered in the first three classes, all of them being invariant under the discrete symmetry CPS which is the product of ordinary CP with d ↔ s switch [52].

Application to the BSM ∆S = 1 Operators
As a consequence of our counting in the previous subsection, by orthogonality we are left with 20 − 7 = 13 four-quark BSM operators linearly independent from the SM ones and violating CPS symmetry, namely -one in class A: -two in class B: -two in class C: -eight in class D: Finally we want to emphasize the following virtue of the chosen basis of BSM operators in which, in contrast to the SM basis in (5) and (6), no summation over quarks is performed. Indeed when considering various extensions of the SM it often turns out that a new heavy mediator, vector or scalar, couples at tree level only to up-quarks or down-quarks but not to both. Moreover right-handed up-and downquark couplings to new mediators, not related by SU(2) L symmetry, could differ by much from each other. In this manner this basis, to be called DQCD basis, is not only appropriate to the evaluation of hadronic matrix elements but also useful for model building.

Dual QCD basics
The explicit calculation of the contributions of pseudoscalars to hadronic matrix elements of local operators is based on a truncated chiral Lagrangian describing the low energy interactions of the lightest mesons [33,34,53] where is the unitary chiral matrix describing the octet of light pseudoscalars and transforming as U → g L U g † R under the chiral symmetry SU (3) L × SU (3) R . The parameter F is related to the weak decay constants F π ≈ 130 MeV and F K ≈ 156 MeV through so that Λ χ ≈ 1.1 GeV. The diagonal mass matrix m involving m u , m d and m s is such that The flavour-singlet η 0 meson decouples due to the large mass m 0 generated by the non-perturbative U (1) A anomaly. Consequently the matrix Π in (40) reads In order to calculate the matrix elements of the local operators in question we need meson representations of colour-singlet quark bilinears. Only currents and densities are directly extracted from the effective Lagrangian in (39). They are given respectively as follows with U turned into U † under parity. As a matter of fact, the chiral correction to densities is meaningful for the Q 6 operator only, as seen in eq. (189). For Q 8 and the other density-density operators, O(p 4 ) mass terms with additional low-energy constants should be introduced in (39). For the sake of consistency, we will thus work in the chiral limit for all of them. The lowest-order chiral realization of tensor bilinears requires two derivatives to get the correct Lorentz structure. It thus involves yet another dimensionful lowenergy constant [56]: In the large N limit, all four-quark operators factorize into two colour singlet bilinears. As a consequence, at O(p 2 ) the only relevant colour-singlet bilinears are and their parity partners, with

Meson evolution in the DQCD basis
The formulae (49)-(52) apply to the strict large N limit and, consequently, to very low scales at which the factorization of matrix elements is valid. In order to obtain the matrix elements at scales O(1 GeV) we have to evolve them with the help of meson evolution. Like in our recent paper on BSM hadronic matrix elements for K 0 −K 0 mixing [39], we will work in the O(p 2 ) chiral limit. While this is a rough approximation, it has been demonstrated there that the pattern of matrix elements evaluated in DQCD at a scale Λ = (0.65 ± 0.05) GeV agrees well with the pattern at µ = 1 GeV obtained from lattice QCD results at µ = 3 GeV through the usual perturbative quark QCD evolution. This makes us confident that the results presented here could also be in the ballpark of the future lattice QCD results for the matrix elements in question. Let us first emphasize that all K → ππ matrix elements of chirally-flipped (relative to the SM operators) mirror operators denoted by a prime can be obtained from the results of RBC-UKQCD collaboration by just reversing their signs. What remains for us is the calculation of the matrix elements which cannot be expressed in terms of SM ones.
The flavour SU (n) generators already introduced in (40) for n = 3 are normalized such that If we again substitute the matrix indices by parenthesis () and brackets [ ] such that each parenthesis/bracket represents now a different flavour index, then the completeness relation among matrices of the fundamental representation of SU (n) simply reads With the background field technology used in [57], the relevant operator evolutions can also be classified according to these identities: with M 2 , an SU (3)-breaking mass term equal to m 2 K for ∆S = 2 transitions and of order (m 2 K − m 2 π ) for ∆S = 1 ones. Had we worked with a nonet of light pseudoscalars, the 1/n term resulting from the purely non-perturbative axial anomaly would have been absent (m 0 = 0).
As already mentioned below (23), the set of operators {sd; ss} may evolve into other ∆I = 1/2 ones. Such is indeed the case in the meson evolution where only the one in Class A induces the SM Q 4 operator. We also observe a reordering of the flavour indices with mixing between operators of classes B and C. Eventually, the mirror operators are again obtained through a parity transformation (U ↔ U † ).
Applied to the 13 BSM operators classified in section 2.3.4, these non-factorizable meson evolutions imply: They agree respectively with the evolution derived in [39] for the O 1,4,5 {sd; sd} operators defined in (18)- (20) and for the O 2 {sd; sd} operator extracted from (21)- (22). However, to infer the meson evolution of the four tensor-tensor operators D * i above the factorization scale, we now have to rely on the SD running pattern.

Quark-gluon evolution in DQCD basis
In order to study the short distance RG evolution from µ = 1 GeV to higher scales we collect here the fundamental equations in the leading order approximation. To see the pattern of this evolution we keep first only the first leading logarithms. One has then for µ 2 > µ 1 where X i denote generically the operators in the DQCD basis. Strictly speaking, only the matrix elements of the X i are scale dependent and the notation X i (µ) is shorthand for with the last expression usually encountered in the literature. We group them in classes I-IV with no mixing under renormalization between various classes. This will also allow us a transparent comparison with the so-called SD basis that we discuss in Section 8. The anomalous dimension matrix (ADM) for all SD operators is then given in the DQCD basis as follows (in units of α s /4π): Class IIγ with the same matrix for the operators B 2 , C 2 , with the same matrix for the operators D 4 , D * 4 . The numerical values of the elements of these ADMs correspond to N = 3 but their explicit N dependence will turn out to be very useful soon.
Of particular interest here are the large entries in the elements (3,2), and (4, 1) in (71) and (2,1) in (72) that imply large mixing of the scalar-scalar operators into tensor-tensor operators. We will see soon that this feature enhances the matrix elements of tensor-tensor operators in the process of O(1/N ) meson evolution. This feature has some analogy to the observation made in [58,59] where the QED short distance RG evolution of NP contributions to charged-current induced leptonic and semileptonic meson decays has been presented, focusing on chirality-flipped operators at the quark level. It has been pointed out that the large mixing of the tensor-tensor operators into the scalar-scalar ones has an important impact on the phenomenology. Recently this aspect has also been discussed in the context of R(D ( * ) ) anomalies in [60,61].
In fact the one-loop QED diagrams responsible for this mixing are the same as the one-loop QCD diagrams with gluon replaced by photon, QCD coupling replaced by QED one and colour matrices replaced by charge ones. Even if the RG evolution of the QED coupling constant is different from the QCD one, the pattern of mixing analysed in [58,59] is very similar to the one in (71) 1 .
The reason why in [58,59] tensor operators have the impact on the scalar ones, as opposed to the case discussed by us, is simply related to the known fact that while the evolution of the matrix elements of operators is governed by the ADM of operators, the evolution of their Wilson coefficients, analysed in [58,59], is governed by the corresponding transposed matrix.
While in [58][59][60][61] the large mixing in question had impact on the phenomenology of B-meson decays, in our case it will have significant impact on ε /ε.

BSM matrix elements in DQCD basis
Having established the meson evolution from the factorization scale (corresponding in the chiral limit to Λ = 0) to Λ = O(1 GeV), what remains to be done is the calculation of the matrix elements of all 13 four-quark BSM operators in the large N limit that here will be generically denoted by with I = 0, 2 being strong isospin and X i = A, B 1,2 , .... When calculating the non-zero matrix elements of the BSM operators in Class B and D, we have to take into account the fact that the partial ππ|U ds |K 0 contribution to on-shell K → ππ decay amplitudes is precisely canceled by a non-local pole diagram involving the strong K 0 → ππK 0 vertex followed by theK 0 annihilation into the vacuum through the non-vanishing K 0 |U ds |0 weak matrix element. Considering charge conservation and Lorentz invariance, the recipe is to simply neglect any contribution from the identity when expanding the U uu,dd,ss components of density-density operators in Classes B and D. In giving the values of the matrix elements, we drop the overall +i factor that is immaterial for physical applications. In the large N limit, the non-vanishing BSM matrix elements in the DQCD basis are then given as follows (h = 1): with the chiral enhancement factor r 2 (m 2 K − m 2 π ) as seen from (42) and F ∼ F π as seen from (41). Those of the mirror operators differ by sign only.
The matrix elements for the operator D 4 , containing three s-quarks, as well as for the tensor-tensor operators D * 1,2,3,4 involving at least four derivatives vanish: We are thus left with 8 four-quark BSM matrix elements for a given chirality and isospin, each of them expressed in terms of either (m 2 K − m 2 π ) or r 2 (µ) in the chiral limit considered.

Summary of hadronic matrix calculations
We have now completed the calculation of 13 BSM hadronic matrix elements evaluated at the cut-off scale Λ, which is governed by the general formula with the coefficients a ij to be extracted from (60)-(66), (83) and (84) and with the matrix elements X j (0) I collected above. In evaluating these matrix elements one should set n = 3.
In the numerical evaluation of matrix elements we will deal with two scales, Λ explicitly seen in (95) and µ in r(µ) hidden in X j (0) I . In this context it is useful to make the following comments: • The Λ dependence is present only in the non-factorizable part of the matrix elements as given above. The scale µ present in r(µ) is at this stage not related to Λ.
• As far as meson evolution in the chiral limit is concerned, there is no distinction between the matrix elements X i 0 and X i 2 so that this distinction is fully described by the values of these matrix elements in the large N limit, that is X j (0) I .
Concerning the value of Λ we will set it at 0.7 GeV. Evaluating then r(µ) at µ = 1 GeV, we will interpret the resulting values of matrix elements as valid at µ = 1 GeV. From there on we will use the standard renormalization group evolution as summarized in Section 5, thereby summing this time leading logarithms to all orders of perturbation theory.

Number of BSM matrix elements to be calculated
In principle one has to evaluate 13 matrix elements for a given chirality and isospin at the factorization scale. They are given in (86)-(94). However by definitionss bilinears do not contribute to the K → ππ matrix elements at the factorization scale. Consequently the matrix elements of D 4 and D * 4 vanish, while the matrix elements of A and B 2 and C 2 can then be expressed in terms of the SM ones as follows wherever sΓ A s = 0. Therefore, in the general case the number of BSM matrix elements for a given chirality and isospin one has to evaluate at the factorization scale is reduced to 8. These are the matrix elements of However, in the chiral limit this number reduces to 3. Indeed, at the factorization scale the matrix elements of D * 1−3 vanish and we have the relations: such that it is sufficient to calculate the matrix elements of But above the factorization scale the situation changes as can be seen from the meson evolution. In particular, the relation between matrix elements of B 1 and D 3 is violated while the one between D 1 and D 2 is preserved. Most importantly the matrix elements of tensor-tensor operators D * 1−3 become non-zero.

SD basis
While the short distance evolution for scales above µ = 1 GeV can be performed in the DQCD basis as demonstrated above, for short distance renormalization group evolution a different basis of 13 BSM operators based on [32] is more useful 2 . That basis, to be termed SD basis in what follows, is closer to bases used in various computer codes present in the literature (see Appendices A.3 and A.4). So, for completeness we would like to present our results in that basis as well. In the SD basis, like in the DQCD basis, no summations over quark flavours are performed, but in contrast to the DQCD basis, colour non-singlet operators are present. While this is a disadvantage with respect to the DQCD basis as far as calculations of hadronic matrix elements are concerned, it turns out to be more suitable for quark-gluon evolution. The 13 BSM independent operators in the SD basis are also linearly independent from the SM ones, in particular none of them mixes into QCD-and QED-penguin operators Q 3,...10 . Using the notation of [32] they are given as follows: Class I: Class II Class III Class IV The connection between the operators in the DQCD and SD bases is rather simple and given in Appendix A. To obtain these relations, Fierz identities in (13)- (17) have been used. Having these relations and the expressions for meson evolution in the DQCD basis, it is straightforward to obtain analogous evolution in the SD basis. It is given in Appendix B. The large N hadronic matrix elements of the 13 BSM operators in the SD basis are collected in Appendix D and one-loop anomalous dimension matrices in Appendix C.
Similar to the discussion in Section 7.3, the number of matrix elements one has to calculate at the factorization scale is reduced for the following reasons. Asss bilinears do not contribute to the matrix elements at the factorization scale, the matrix elements of Q SLL,s wherever sΓ A s = 0. Therefore, in the general case the number of BSM matrix elements for a given chirality and isospin one has to evaluate at the factorization scale is reduced as in the DQCD basis to 8: However in the chiral limit, analogous arguments to those presented in Section 7.3 imply that only the matrix elements of the following operators have to be evaluated This can be easily verified by using the result in (101) together with the connection between the operators in the DQCD and SD bases given in Appendix A.

SMEFT view of BSM operators
Following [47,48], it should be emphasized that while generally 13 BSM operators are consistent with the SU(3) c × U (1) Q symmetry, only 7 operators are consistent with the full SM gauge symmetry SU(3) c × SU(2) L × U (1) Y [64, 65] 3 . As in [48] a different operator basis has been used, we identify here these 7 operators in the SD and DQCD bases. Beginning with the SD basis, the operators in (103) and (104) and in class IV violate the full SM gauge symmetry SU(3) c × SU(2) L × U (1) Y [64,65]. In particular one can easily check that they do not conserve hypercharge. Consequently, if no new particles close to electroweak scale exist and the SMEFT is the correct description between the electroweak scale and the NP scale, then the Wilson coefficients of the operator (104) and those in class IV must vanish. Though explicitly carrying a nonzero weak hypercharge (Y = 2), the operator (103) can in principle be generated through the dimension-six gauge-invariant operator i(φ † D µ φ)(ūγ µ P R d) frozen at the vacuum expectation value of the Higgs doublet φ (Y = 1). Such would for example be the case in an SU (2) L × SU (2) R × U (1) B−L UV completion of the SM via the tree-level W L − W R mixing. If we neglect this modified W L coupling [64], the full set of contributing linearly independent operators contains not 40 but 28 operators: • 7 SM operators and the corresponding 7 mirror operators with L and R interchanged, • 7 BSM operators, those in classes I and III and the operators in (105) and (106), and the corresponding 7 mirror operators obtained again by interchanging L and R.
Similarly, in the DQCD basis the imposition of invariance under the full gauge group of the SM reduces the number of contributing BSM operators of a given chirality to 7. In Table 1 Table 1: SD and DQCD BSM operators generated in SMEFT.
As discussed already at the end of Sections 7.3 and 8, the number of matrix elements one has to evaluate at the factorization scale for given chirality and isospin is reduced to 8 when one takes into account thatss bilinears do not contribute to the matrix elements at this scale. They are listed in (99) However, in the chiral limit even without the imposition of SMEFT, the number of matrix elements to be calculated at the factorization scale in each basis is reduced to 3. They are given in (101) and (119). Consequently, once SMEFT is imposed only one matrix element at the factorization scale in each basis for a given chirality and isospin has to be evaluated. For instance Once meson evolution is turned on the picture is much reacher as some relations between matrix elements valid at the factorization scale are broken. Most important, the matrix elements of colour singlet tensor-tensor operators do not vanish any longer.

Numerical Results
Having all these results for hadronic matrix elements in DQCD and SD bases at hand, we will next present their numerical values. To this end we will assume that values of the matrix elements at Λ = 0.7 GeV give an adequate representation of their values at 1 GeV. This treatment was rather successful in the case of our analysis of K 0 −K 0 matrix elements in [39] and we expect that it is a reasonable approximation for the time being. This is furthermore supported by the fact, as discussed at the end of Section 4, that the meson evolutions for operators discussed in the present paper have similar structure to the ones in [39].
Next two features should be noticed: • The BSM hadronic matrix elements of classes B and D, involving the factor r 2 (µ), are chirally enhanced with and the corresponding operators have the highest potential to have an impact on ε /ε without requiring large values of their Wilson coefficients.
• Among the chirally enhanced matrix elements of classes B and D, those contributing to the isospin amplitude A 2 are most important as the contribution of the amplitude A 0 to ε /ε is automatically suppressed by a factor 1/22 relative to the one of A 2 since Consequently as pointed out in [15], in order to obtain a significant contribution to this ratio from NP contributing to A 0 the imaginary part of the corresponding Wilson coefficient should be larger than in the case of A 2 in order to compensate this suppression. This in turn could lead, in certain NP models, to the violation of the present bounds on rare decays.
• The values of I = 0 and I = 2 matrix elements of all BSM operators are at µ = O(1 GeV) similar to each other so that we do not expect these new operators to be relevant for the ∆I = 1/2 rule. This is an important result as it implies that either this rule is fully governed by the SM dynamics or NP contributions to the A 0 amplitude, at the level of (10 − 20)%, come from modifications of the Wilson coefficients of the SM operators. See [10] for a detailed analysis.
Now among the hadronic matrix elements of SM penguin operators calculated by lattice QCD and DQCD, the most important ones are the matrix elements of Q 6 and Q 8 which are given as follows [33,66,67] with the values of B from RBC-UKQCD collaboration and similar results from DQCD [3]. It will be then of interest to compare the values of the matrix elements of BSM operators with these two matrix elements for A 0 and A 2 , respectively. In doing this we have to take into account that Q 6 and Q 8 have the (V −A)⊗(V +A) Dirac structure, while the BSM operators considered by us involve the P L and P R chiral projectors. In order to compensate for this, we will give the results for For the numerical analysis of the matrix elements SD evolution we use WCxf [68] and wilson [69] as well as the results from the previous section. Below the hadronic scale, the running of the matrix elements is given by the meson evolution formulae given in previous sections.
We use the following input  Tables 2 and 3 we show the results in the DQCD basis for operators allowed by SMEFT and forbidden by it, respectively. The corresponding results for the SD basis are given in Tables 4 and 5. Let us then extract the most important lessons from these tables.
• Concentrating first on I = 2 matrix elements, that according to (123) could have larger impact on ε /ε than I = 0 matrix elements, we pin down several BSM operators for which the values of matrix elements are in the ballpark of the value of Q 8 (µ) 2 for µ = O(1 GeV). As the Wilson coefficient of Q 8 is O(α e ), it is conceivable that some of these operators could have significant impact on ε /ε. In the DQCD basis this is the case of the pure tensor-tensor operators with the first two allowed by SMEFT.
In the SD basis this is the case of the operators with all three allowed by SMEFT.
• Looking then at I = 0 matrix elements, we bring out a number of matrix elements which are comparable to or larger than the one of Q 6 (µ) 0 . In the DQCD basis this is the case not only for the operators in (128) but also for B 1,2 and D 3 . In the case of the SD basis the additional large I = 0 matrix elements are those of Q SLR,u 2 and Q SLL,d 2 .
Explicitly we have at µ = 1 GeV. In the DQCD basis the meson evolution from the factorization scale to scales O(1 GeV) is primarily responsible for this result. It should however be noticed that while the matrix elements of Q 6 , Q 8 and generally scalar-scalar operators increase with increasing µ, the ones of tensor-tensor operators are only weakly dependent on µ for µ > 1 GeV. Therefore the numerical factors in (130) and (131) generally decrease significantly with increasing µ. An exception is the operator Q SLL,u 3 which, among tensor-tensor operators in both operator bases, is the only one which is colour-non-singlet. We will return to this point at the end of the present section.
It is of interest to understand why the matrix elements of the remaining tensortensor operators exhibit such a weak dependence on µ. This is most transparently seen by studying the RG formula (67) for X i = D * 1 together with the ADM in (71), although it can already be suspected at the one-loop level, from (78) and (83) which imply D 2,1 − 1 4 D * 2,1 = 0 at Λ = 2πF ≈ 0.8 GeV. Now the matrix element D * 1 (µ 1 ) with µ 1 = 1 GeV is generated, as seen in (83), in the process of meson evolution, through mixing with the scalar-scalar operator D 2 . Because this mixing is very large and further enhanced through short but fast meson evolution D * 1 (µ 1 ) is, as seen in Table 2, much larger than D 2 (µ 1 ) . The values of D * 1 (µ 2 ) , for µ 2 > µ 1 are now governed first of all by the self-mixing of D * 1 , the (3,3) in (71) and the mixing of D * 1 with D 2 given by the entry (3,2) in (71). But while the (3, 3) entry is much smaller than (3,2), the matrix element of D * 1 is much larger than that of D 2 as mentioned above. Using this information in   Table 2: Matrix elements X i 0,2 of BSM operators in the DQCD basis allowed by SMEFT contributing to the isospin amplitudes A 0,2 in units of GeV 3 for four values of µ. In the last column we give the values of the matrix elements Q 6 (µ) 0 and Q 8 (µ) 2 for comparison.   Table 3: Matrix elements X i 0,2 of BSM operators in the DQCD basis allowed by SU(3) c × U (1) Q but not by SMEFT contributing to the isospin amplitudes A 0,2 in units of GeV 3 for four values of µ. In the last column we give the values of the matrix elements Q 6 (µ) 0 and Q 8 (µ) 2 for comparison. (67) we find that these effects cancel each other to a large extend leaving a very weak µ dependence of D * 1 (µ). This cancellation even improves with increasing µ because, while the diagonal evolution of D * 1 slowly decreases its matrix element, the one of D 2 governed by the large entry (2, 2) in (71) and having opposite sign to (3,3) increases this matrix element. The evolution of D * 1 is governed in addition by the entry (3, 1), the mixing of D * 1 and D 1 , and by the entry (3, 4), the mixing of D * 1 and D * 2 . Taking into account that the matrix elements of D * 2 and D * 1 are equal to each other and the same applies to matrix elements of D 1 and D 2 , one can then easily check using (67) that these two additional effects cancel each other to a large extent and have only a very small impact on the evolution of D * 1 . This result remains true after performing the summation of leading logarithms to all orders of perturbation theory, which all the numerical values in the tables are based on.
In order to further demonstrate that the pattern of meson evolution agrees with  the SD evolution, it is useful to normalize our results for BSM matrix elements of scalar-scalar operators to Q 8 (µ) 2 and consider the ratios with X i denoting any scalar-scalar operator either in DQCD or SD basis. This has the virtue that in the case of chirally enhanced matrix elements of scalar-scalar operators the dominant µ-dependence present in r(µ) cancels out, exhibiting the non-factorizable µ dependence of R I for such operators. Inspecting the Tables 6 and 7 for the DQCD basis and Tables 8 and 9 for the SD basis, we can make the following observations: • The ratios R I for scalar-scalar matrix elements governed by r 2 (µ), that is of B 1,2 and D 1,2,3 in the DQCD basis, all decrease with increasing µ so that the SD evolution matches well the meson one. In particular, the SD evolution of B 1,2 is in the range 1 − 3 GeV by roughly a factor of two slower than that of D 1,2,3 in agreement with the meson evolution equations of Section 4. The matrix elements of A and C 1,2 are not chirally enhanced and the ratios R I are not useful in this case for exhibiting the proper matching of SD and meson evolutions. However, inspecting the SD and meson evolutions, also in this case the matching is good. in the SD basis, all governed by r 2 (µ), exhibit also the SD pattern as expected from the matching of SD evolution to the meson evolution summarized in Appendix B. All decrease with increasing µ with the speed expected from meson evolution. The comments made on matrix elements not enhanced by r 2 (µ) in the DQCD basis applies also to the SD basis.
In the case of tensor-tensor operators the diagonal evolution of these operators differs from the one of scalar-scalar operators simply because their anomalous dimensions, as seen from (80), are very different and read for N = 3 where Q T denotes any colour singlet tensor-tensor operator. Consequently, the diagonal SD evolutions of Q 8 (µ) and Q T (µ) are rather different with β 0 = 11−2f /3. Moreover, as we discussed above, the evolution of tensor-tensor operators is not governed by their diagonal evolution but rather by a complicated ADM. Therefore, the ratios in (132) are not useful in this case for the demonstration of the matching of non-factorizable evolutions and we do not show them in Tables  6-8. An exception is the colour non-singlet operator Q SLL,u 3 , as already mentioned above. In this case, in order to exhibit better non-factorizable SD evolution of matrix elements of this operator, it appears to be more appropriate to consider the ratio  with γ 3 extracted from the (3, 3) entry in (186). We show this ratio in Table 8.
One can easily check that the decrease of this ratio with increasing µ matches well the meson evolution of this operator. It is also interesting to note that at the one-loop level (79) together with (84) imply D 3 − D * 3 /4 = 0 at a higher value of Λ = √ 6πF = 1 GeV compared to (78) and (83) for which a similar cancellation occurs at Λ = 2πF ≈ 0.8 GeV.

Summary and outlook
Motivated by the hints for NP contributing to the ratio ε /ε we have calculated hadronic matrix elements of 13 BSM four-quark operators in DQCD, including the meson evolution in the chiral limit. This is the first calculation of these matrix elements to date, thereby opening the road to the investigations of ε /ε and K → ππ decays beyond all BSM analyses found in recent literature  in which NP affected only Wilson coefficients of SM operators.
Our main messages to take home are the following ones: • The pattern of long distance evolution (meson evolution) matches once more the one of short distance evolution (quark-gluon evolution), a property which to our knowledge cannot be presently achieved in any other analytical framework. It should be emphasized that this important result has been obtained for 13 operators without any free parameter except possibly the physical cutoff Λ which in any case has to be chosen in our framework in the ballpark of 0.7 GeV. See Section 10 for details.
• Several matrix elements and in particular those of tensor-tensor operators have values at µ = O(1 GeV) in the ballpark of the ones of the dominant electroweak penguin matrix element Q 8 (µ) 2 . Therefore they could have large impact on ε /ε.
• The mixing of the scalar-scalar operators into tensor-tensor operators in the process of meson evolution is responsible for this result so that the role of     scalar-scalar operators should not be underestimated. They can be most easily generated through tree-level exchanges of heavy colourless and coloured scalars.
Clearly, in order to identify the most important operators, also their Wilson coefficients in a given NP scenario must be known. But having to disposal large K → ππ matrix elements identified in our paper will, in some NP scenarios, facilitate the removal of ε /ε anomaly without violating other constraints. A detailed analysis of ε /ε with and without the imposition of SMEFT and concentrating rather on the structure of Wilson coefficients of BSM operators resulting from the renormalization group effects is presented in an accompanying paper [48] while a master formula for ε /ε beyond the SM has been recently presented in [47].
Our calculation of meson evolution has been performed in the chiral limit but, as we already mentioned in the Introduction, DQCD could reproduce some lattice results in this approximation. Still it is desirable to extend our work beyond chiral limit, a goal which could be reached before lattice QCD calculations for the hadronic matrix elements in question are expected to be available.
In the case of confirmation of ε /ε anomaly predicted in DQCD by more precise lattice QCD calculations, our results will play an important role in selecting those NP models that can provide sufficient upward shift in ε /ε in order to explain the data.
On the other hand, if future lattice QCD calculations within the SM will confirm the data on ε /ε, our results will put strong constraints on the parameters of a multitude of NP models.
was done and financially supported by the DFG cluster of excellence "Origin and Structure of the Universe".

A Basis transformations A.1 From SD to DQCD basis
Here we list the basis transformation between the SD and DQCD bases. (136) Class C Class D A.2 From DQCD to SD basis Class II Class III Class IV with ' on C 1 , B 1 meaning mirror partner (L ↔ R) and * on D 1,2,3,4 meaning tensortensor partner of P L − P L .

A.3 From Flavio to SD basis
To use wilson for the SD matrix element evolution we need to translate the used SM and BSM operators into a predefined basis. We choose the flavio basis [70], where the operators are defined in [71]. The transformation from the flavio basis into SM basis reads:

A.4 From SD to Flavio basis
Here we report the inverse transformation of the relations in the previous subsection.

B Meson evolution in the SD basis
Using the notationΛ = Λ/(4πF ) for short, we have the following evolution: Class I: Class IV with analogous equations for Q SLL,s 1,2 .

C Quark-gluon evolution in the SD basis
The anomalous dimension matrices for all BSM operators are then given in the SD basis as follows (in units of α s /4π) [32,58].

Class Iγ
Class IIγ where N denotes the number of colours with N = 3.

Class III
In with the same matrix for the operators Q SLL,s 1,2 .

D Hadronic matrix elements in the SD basis
In the large-N limit the matrix elements of the two most important SM operators are given as follows [33,66,67] Q 6 (µ) 0 = − r 2 (µ)(F K − F π ) , Q 6 (µ) 2 = 0 , The matrix elements of BSM operators are listed below. We omit for brevity the argument Λ = 0.