BSM hadronic matrix elements for ε′/ε and K → ππ 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 → ππ decays and in particular to the ratio ε′/ε 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 ε′/ε 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 ΔI = 1/2 rule.

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) Q 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

2)
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:

5)
Electroweak penguins (2.6) 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.

BSM operators 2.3.1 Chiral Fierz identities
The following 16

JHEP02(2019)021
Within these conventions, the corresponding dual basis {Γ A } = P L , P R , γ µ P R , γ µ P L , 1 2 σ µν , (2.9) 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 (2.11) 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: 14) Class C:

15)
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 Fierzconjugate flavour indices {ab; cd} = {sd; sd}:

JHEP02(2019)021
Using this technology, we succeeded [32,39] in turning the SUSY basis [50,51]  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 (2.4)-(2.6). For that purpose let us first isolate the latter ones by selecting the appropriate combinations of operators in (2.23), with sums over q = u, d, s understood.
Class A: (2. 25) 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 (2.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}.

JHEP02(2019)021
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 (2.5) and (2.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 downquarks but not to both. Moreover right-handed up-and down-quark 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.

JHEP02(2019)021 3 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] is the unitary chiral matrix describing the octet of light pseudoscalars and transforming as 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 (3.2) 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 (3.1). 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. (D.1). For Q 8 and the other JHEP02(2019)021 density-density operators, O(p 4 ) mass terms with additional low-energy constants should be introduced in (3.1). 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 low-energy 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 (3.11)-(3.14) apply to the strict large N limit at which the factorization of matrix elements is valid. In this limit it is not possible to determine the scale associated to these matrix elements. To this end one has to calculate non-factorizable contributions represented by loops in the meson theory. The factorization scale is then the scale at which these non-factorizable contributions vanish. Quite generally the factorization scale is found to be at very low scales O(m π ) and in order to obtain the matrix elements at scales O(1 GeV) one has to evolve them with the help of the meson evolution. As in our paper this evolution will be performed in the chiral limit, the factorization scale is simply at zero momentum. 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 with This makes us confident that the results presented here on the basis of the same 1/N counting and the same Λ-to-µ identification have rather similar uncertainties. 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 (3.2) for n = 3 are normalized such that Tr (λ α λ β ) = 2δ αβ . 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) With the background field technology used in [57], the relevant one-loop operator evolutions from the factorization scale taken at zero momentum (denoted 0) to the cut-off-momentum (denoted Λ) can also be classified according to these identities:

JHEP02(2019)021
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 pseudo-scalars, the 1/n term resulting from the purely non-perturbative axial anomaly would have been absent (m 0 = 0).
As already mentioned below (2.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 (2.18)-(2.20) and for the O 2 {sd; sd} operator extracted from (2.21)-(2.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.

JHEP02(2019)021 5 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. 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π): with the same matrix for the operators B 2 , C 2 ,

Class IIIγ
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 (5.4) and (2, 1) in (5.5) 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 tensortensor operators in the process of O(1/N ) meson evolution. This feature has some analogy

JHEP02(2019)021
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 (5.4). 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 ε /ε.

JHEP02(2019)021
• dropping the subleading 1/N terms in the ADMs (5.3)-(5.5): From (4.10)-(4.16), these non-factorizable evolutions are compatible with the longdistance ones in the nonet approximation (m 0 → 0) with since the tensor-tensor operators vanish at the factorization scale at O(p 2 ). Comforted by such a consistent matching of mixing pattern, we extend the octet meson evolution (m 0 → ∞) to the non-factorizable tensor-tensor operators as follows:

Large N limit
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 onshell 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 JHEP02(2019)021 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 (3.4) and F ∼ F π as seen from (3.3). 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 (4.10)-(4.16), (6.11) and (6.12) 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 (7.11) 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 Λ.

JHEP02(2019)021
• 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 (7.2)-(7.10). 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 A I = 1 6 ( Q 3 I − Q 9 I ), (7.12) 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 JHEP02(2019)021 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 independent BSM 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:

9)
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 (2.13)-(2.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.

JHEP02(2019)021
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 (7.17) together with the connection between the operators in the DQCD and SD bases given in appendix A.  coupling [64], the full set of contributing linearly independent operators contains not 40 but 28 operators:

SMEFT view of BSM operators
• 7 SM operators and the corresponding 7 mirror operators with L and R interchanged, Once meson evolution is turned on the picture is much richer 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 JHEP02(2019)021 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]  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

JHEP02(2019)021
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 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 (10.7) 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    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.
significantly with increasing µ. An exception is the operator Q SLL,u 3 which, among tensortensor 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 tensor-tensor operators exhibit such a weak dependence on µ. This is most transparently seen by studying the RG formula (5.1) for X i = D * 1 together with the ADM in (5.4), although it can already be suspected at the one-loop level, from (6.6) and (6.11) 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 (6.11), 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 (5.4) and the mixing of D * 1 with D 2 given by the entry (3,2) in (5.4). 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 (5.1) 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 (5.4) 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 (5.1) 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.

JHEP02(2019)021
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 (6.8), 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 (10.11) 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

JHEP02(2019)021
with γ 3 extracted from the (3, 3) entry in (C.5). 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 (6.7) together with (6.12) imply D 3 − D * 3 /4 = 0 at a higher value of Λ = √ 6πF = 1 GeV compared to (6.6) and (6.11) 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 cut-off Λ 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 ε /ε.

JHEP02(2019)021
• 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 JHEP02(2019)021  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]. Furthermore the hadronic matrix elements computed in this paper are implemented in the open-source codes flavio [70] and smelli [71], which allow for a numerical analysis including ε /ε.
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.

A.1 From SD to DQCD basis
Here we list the basis transformation between the SD and DQCD bases.
Class C Class III Class IV with on C 1 , B 1 meaning mirror partner (L ↔ R) and * on D 1,2,3,4 meaning tensor-tensor 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 [72]. The transformation from the flavio basis into SM basis reads: Class III  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]. 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 , (D.1) The matrix elements of BSM operators are listed below. We omit for brevity the argument Λ = 0.

JHEP02(2019)021
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.