Singularity-free Next-to-leading Order $\Delta S= 1$ Renormalization Group Evolution and $\epsilon_{K}^{\prime}/\epsilon_{K}$ in the Standard Model and Beyond

The standard analytic solution of the renormalization group (RG) evolution for the $\Delta S = 1$ Wilson coefficients involves several singularities, which complicate analytic solutions. In this paper we derive a singularity-free solution of the next-to-leading order (NLO) RG equations, which greatly facilitates the calculation of $\epsilon_K^{\prime}$, the measure of direct $CP$ violation in $K\to \pi\pi$ decays. Using our new RG evolution and the latest lattice results for the hadronic matrix elements, we calculate the ratio $\epsilon_{K}^{\prime}/\epsilon_{K}$ (with $\epsilon_{K}$ quantifying indirect $CP$ violation) in the Standard Model (SM) at NLO to $\epsilon_{K}^{\prime}/\epsilon_{K} = (1.06 \pm 5.07) \times 10^{-4} $, which is $2.8\,\sigma$ below the experimental value. We also present the evolution matrix in the high-energy regime for calculations of new physics contributions and derive easy-to-use approximate formulae. We find that the RG amplification of new-physics contributions to Wilson coefficients of the electroweak penguin operators is further enhanced by the NLO corrections: If the new contribution is generated at the scale of 1-10 TeV, the RG evolution between the new-physics scale and the electroweak scale enhances these coefficients by 50-100 %. Our solution contains a term of order $\alpha_{EM}^2/\alpha_s^2$, which is numerically unimportant for the SM case but should be included in studies of high-scale new-physics.


Introduction
The parameter K / K is the ratio of the measures of direct and indirect charge-parity (CP ) violation in the Kaon system. While indirect CP violation is a per-mille effect in the Standard Model (SM), K is smaller by another three orders of magnitude than K , with | K | ∼ O(10 −6 ). A strong suppression by the Glashow-Iliopoulos-Maiani (GIM) mechanism and an accidental cancellation of leading contributions in the Standard Model makes K / K highly sensitive to new physics. The first element of the SM prediction for K is the calculation of initial conditions for Wilson coefficients and their renormalization group evolution from the electroweak scale (of the order of W and top mass) down to the hadronic scale of order 1 GeV, at which hadronic matrix elements are calculated. These steps purely involve perturbative methods and have been carried out to leading order (LO) in the strong coupling constant α s in Refs. [1][2][3][4]. The next-to-leading order (NLO) involves the electromagnetic coupling α EM 1/128 [5][6][7][8], the next higher order in α s [9][10][11], and order α EM α s [11][12][13]. In terms of isospin amplitudes K is given by (see e.g. Ref. [14]) where A I ≡ (ππ) I |H |∆S|=1 eff |K 0 are isospin amplitudes and ω + = (4.53 ± 0.02) × 10 −2 (see Refs. [14,15] for the precise definition), | K | = (2.228 ± 0.011) · 10 −3 , and ReA 0 = (3.3201 ± 0.0018) × 10 −7 GeV are taken from experiment.Ω eff = (14.8 ± 8.0) × 10 −2 parameterizes isospin-violating contributions [15,16].
The |∆S| = 1 nonleptonic effective Hamiltonian for weak decays in the Standard Model is given by [13] where λ u = V * us V ud and τ = −V * ts V td / (V * us V ud ). The operator basis Q i comprises ten operators which are defined in Ref. [13]; the current-current operators Q 1 and Q 2 the QCD-penguin operators Q 3 to Q 6 and the QED-penguin operators Q 7 to Q 10 where V ∓ A represents γ µ (1 ∓ γ 5 ), α and β denote color indices, and e q is the electric charge of the quark q. The corresponding Wilson coefficients z i and v i (or y i ) serve as effective couplings to these effective operators. By virtue of the framework of effective theories, the parameter µ splits short distance from long distance scales, effectively separating the perturbative high energy regime from the non-perturbative realm of low energy QCD. Taking up the perturbative part of the calculation, the Wilson coefficients have been determined through matching calculations up to next-to-leading order at the scale M W [13]. The calculation of the hadronic matrix elements, being non-perturbative quantities, is a major challenge and has recently been performed on the lattice with unprecedented accuracy [17][18][19][20].
The combination of these calculations into a prediction for K / K requires a treatment within renormalization group (RG) improved perturbation theory to sum up large logarithms. However, it is known that the analytic determination of the required evolution matrix at the next-to-leading order suffers from singularities appearing in intermediate steps of the calculation, which make a computational evaluation highly laborious and complicated. The standard way to solve the NLO RG equations requires the diagonalization of the LO anomalous dimension matrixγ (0) s and the NLO correction involves fractions whose denominators contain the differences of eigenvalues ofγ (0) s . Some of these denominators vanish and are usually regulated in the numerical evaluation [11,21]. In Ref. [22] an analytic solution for the RG equations which is free of singularities is presented. This solution involves the diagonalization ofγ (0) s and gives explicit prescriptions to handle the different cases in which the formulae of Refs. [11,21] develop singularities.
In this paper, we present a new singularity-free solution which permits an easy and convenient numerical implementation. Instead of singularities our analytic formula has undetermined parameters. However, we will show that these spurious parameters cancel and leave the evolution matrix unambiguous. Unlike the solution of Ref. [22] our new formula requires neither the diagonalization ofγ (0) s nor a distinct treatment of the part of the RG evolution which involves the spurious singularities. Using our new RG evolution and the latest lattice results [17][18][19][20], we calculate the K / K in the Standard Model at next-to-leading order to find a value which is below the experimentally measured quantity by 2.8 σ.
The second objective of this paper is the derivation of a useful formula for the calculation of new physics contributions to K / K , in which we evaluate the evolution matrices for scales far above the electroweak scale. To this end we identify a contribution of order α 2 EM /α 2 s in the evolution matrix which can become relevant for studies of TeV-scale new physics, because α s decreases with increasing scale. We observe an approximately logarithmic behavior of the evolution matrix as a function of the energy scale above the electroweak scale. This paper is organized as follows. In Sec. 2, we briefly review the RG evolution of the |∆S| = 1 effective Hamiltonian at the next-to-leading order. We give a detailed analysis of the evolution matrix and its singularities and provide a new analytic solution without singularities. Then we evaluate K / K in the Standard Model at the next-to-leading order in Sec. 3. In Sec. 4, we work out the evolution matrices in the high-energy regime explicitly for calculations of new physics contributions. The last section is devoted to conclusions and discussion.

Renormalization Group Evolution of the ∆S = 1 Hamiltonian
In this section, we review the singularities in the RG evolution of the |∆S| = 1 effective Hamiltonian at the next-to-leading order. Then we generalize the analytic ansatz of the RG evolution given in the literature and present a solution, which is finite at all stages of the calculation. Our solution contains free parameters, which we show to cancel from the evolution matrix, and compare our singularity-free solution with the standard results from the literature.

Singularities in the Evolution Matrix
The evolution of the Wilson coefficients v i and z i from the W boson mass and the charm mass respectively to the hadronic scale µ are given by whereÛ f (µ 1 , µ 2 ) is the RG evolution matrix from µ 2 down to µ 1 and f is the number of the active flavors between these two energy scales. The matricesM c,b represent matching matrices between effective theories with different numbers of flavor and are given in Ref. [13]. Although the effect of the running of α EM is numerically negligible for K / K in the Standard Model [13], we consider this effect to cover new-physics scenarios with largely separate scales. The general form of the evolution matrix is given by [23,24], with the g s -ordering operator T gs and the anomalous dimension matrixγ and the QCD β function. The expansions of the latter two quantities and α EM up to NLO read:γ where β 0 = 11 − 2f /3, β 1 = 102 − 38f /3, β se 1 = −8/9(u + d/4), and β e 0 = −4/3(4u/3+ d/3 + ) are the leading and next-to-leading coefficients of the QCD and QED beta functions, and u, d, are the numbers of the active up-type-quark, down-type-quark, and charged-lepton flavors (f = u + d).γ (0) s is the LO QCD anomalous dimension matrix, and the NLO corrections consist of the three remaining matrices,γ (0) e ,γ (1) s , andγ (1) se , which are the leading QED, next-to-leading QCD, and combined QCD-QED anomalous dimension matrices, respectively.
The ansatz for the NLO evolution matrix (with µ 1 < µ 2 ) is given by [11,21] whereK and the LO evolution matrix where the QED contributions to the beta functions (β se 1 , β e 0 ) are discarded in this subsection 2.1.
The matricesK(µ 1 ) andK (µ 2 ) encode the NLO corrections and depend on the number of active flavors through the beta function and the anomalous dimension matrices. The matricesĴ e ,Ĵ s andĴ se govern the leading electromagnetic, next-toleading strong, and next-to-leading combined strong-electromagnetic contributions to the RG evolution.
Differentiating Eqs. (15) and (11) with respect to g s (µ 1 ) yields the following differential equation forK(g s (µ 1 )) [9,23], The traditional ansatz in the literature is to takeĴ e ,Ĵ s andĴ se as constant matrices for any fixed number of flavors. The differential equation (19) then implies the following equations for the matricesĴ e ,Ĵ s andĴ se [11], It is well known, however, that Eqs. (20) and (21) develop singularities in the case of three flavors. Furthermore, Eq. (22) is even singular for any number of flavors.
We now show how these singularities arise. For this purpose, it is instructional to transform Eqs. (20)- (22) into the diagonal basis ofγ (0)T s . This is a common procedure in the literature since it allows to isolate the singularities and remove them "by hand". We stress that this is only for the purpose of a better understanding of the origin of these singularities. A numerical evaluation of our solution does not require the diagonalisation ofγ We find singular solutions if the difference of two eigenvalues ofγ (0)T s is equal to 2β 0 , which is the case for three flavors:γ (0)T s,D has the elements 2 and −16 and 2β f =3 0 = 18, so that one denominator in Eq. (23) vanishes with a generally non-zero numerator. When we transform Eq. (22) into the same basis we find singular results for i = j and also for degenerate eigenvalues. Nonetheless, once all relevant terms have been joined together, all these singularities cancel and the evolution matrixÛ f (µ 1 , µ 2 ) becomes finite [11]. This procedure, however, requires taking care of each singularity by hand by adopting the aforementioned diagonal basis, then regularizing the singularities and keeping track of them until the end of the calculation. Indeed, Buras et al. have regulated some of the singularities by a logarithmic term [13]. Subsequently, Adams and Lee have proposed a systematical solution for all singularities [25], which, however, still requires the adoption of a certain diagonal basis. The freedom of choosing the order of the eigenvalues on the diagonal ofγ (0)T s,D involves an ambiguity. This can pose a problem in computational implementations, since it is absolutely necessary to use the same diagonal basis as Adams and Lee do, which is not the one which orders eigenvalues by their numerical value. The solution in Ref. [22] follows the same line, after diagonalizinĝ γ (0)T s,D several different cases must be considered: whenever two eigenvalues differ by an integer multiple of 2β 0 a special implementation is required. In the next subsection we propose a solution which does not rely on a specific basis and permits a much faster, easier and, in particular, more stable computational algorithm.
We systematically include O(α 2 EM /α 2 s ) corrections in the RG evolution. This contribution has not been considered in the literature. Although appearing as O(α 2 EM ), these terms can become sizable at high energies because of the awkward 1/α 2 s dependence, making them numerically comparable to O(α s ). We note that this contribution does not receive contributions from higher orders of the anomalous dimension matrix in Eq. (12), but only appears at the next-to-leading order.
With these generalizations we can now solve the differential equation in Eq. (19). Inserting our ansatz into Eq. (19) we obtain the following nine matrix equations for the nine constant matricesĴ: These equations yield finite solutions forĴ. As an effect of the constant matriceŝ J s(,e,se),1 , the analytic singularities of Eqs. (20)- (22) do not occur, because for the problematic matrix elements now both sides of the equations are zero. We stress that one can solve Eqs. (29) to (37) without diagonalizingγ (0)T s ; these equations are mere systems of linear equations for the 100 elements ofĴ s,e,ee,0,1 andĴ se,0,1,2 each, which are quickly solved by computer algebra programs [27]. However, there are multiple solutions in some of the inhomogeneous equations, because the corresponding homogeneous equations have a non-trivial null space. As a consequence, these solutions for J depend on arbitrary parameters, e.g. there are 16 undetermined components in the case of three active flavors. These parameters, however, do not produce any ambiguity in physical results. In the next subsection, we will show that they completely drop out after combining terms of the same order and the evolution matrix in Eq. (15) does not depend on these parameters. Therefore, one can set them to arbitrary values from the beginning. In our calculation of K / K we kept the parameters arbitrary as a crosscheck of the consistency of our calculation.
The procedure to determine the evolution matrix from µ 2 to µ 1 requires algebraically solving the matrix equations (29)-(37) for a given number of active flavors and inserting the solutions into the full evolution matrix in Eq. (15). We use 10 × 10 anomalous dimension matricesγ    (1) se [10][11][12]24]. The solutions for the matricesĴ in the case of three active flavors (with two active leptons) in naive dimensional regularization (NDR) scheme with MS subtraction, are given as follows: where t s , t e , and t se1,2,...,14 are the arbitrary parameters of the matrix equations. Our convention for the matrixV is (γ 10 . Although Eq. (42) makes explicit reference to the the diagonal basis, the term involvingV completely drops out from the evolution matrix (see next subsection), and thereby our solution for the latter does not require any matrix diagonalisation. Our Eqs. (26)-(37) hold in any operator basis. Moreover, if an ordinary four-dimensional basis transformation is applied to Eqs. (4)- (8), the corresponding RG matricesĴ ... can be simply found by transforming those in Eqs. (38)- (46) in the same way asγ (0)T s . If the basis transformation is D-dimensional, meaning that it involves evanescent operators, theĴ ... matrices undergo an additional scheme transformation [26,28]. We collect the solutions for more than three active flavors in Appendix A.

Cancellation of Spurious Parameters
We now present some details of the cancellation of the arbitrary parameters. First, we take a look at the O(α s ) part of the evolution matrix in Eq. (47), In the three-flavor regime, the matrixĴ s,0 in Eq. (38) contains an undetermined component t s . Since the first and second term ofÛ QCD in Eq. (52) depend on different scales, one naively could argue that the cancellation of any dependence has to take place for each term independently of the other. However, we will show that this is not the case.
We locate the undetermined parameter in [Ĵ s,0 ] 8,7 = t s . The matrix product J s,0Û0 (α 1 , α 2 ) naturally contains a dependence on t s in the 8th row. Actually, this dependence does cancel for all elements except for [Ĵ s,0 U 0 (α 1 , α 2 )] 8,7 ⊃ (α 2 /α 1 ) 1/9 t s . The matrix productÛ 0 (α 1 , α 2 )Ĵ s,0 in the second term ofÛ QCD naturally obtains the parameter t s in the 7th column, and again the product consistently cancels this dependence for all entries except for [Û 0 (α 1 , α 2 )Ĵ s,0 ] 8,7 ⊃ (α 2 /α 1 ) −8/9 t s . The full cancellations is thus only achieved by taking both terms of the first line of Eq. (52) into account and takes the form The reason that causes the singularity to arise -eigenvalues ofγ (0)T s differing by 2β 0 in Eq. (23) -is also responsible for the cancellation of the undetermined parameter between the high and low scales. The difference of two eigenvalues ofγ (0)T s by 2β 0 causes a difference of 1 in the exponents of (α 2 /α 1 ) and indeed the spectrum of γ (0)T s /2β 0 contains both 1/9 and −8/9 as eigenvalues. Thus, this difference allows the prefactors α 1 and α 2 of the first two terms in Eq. (52) to exactly cancel these terms between the different scales and entirely independent on the actual size of the scales.
Next, we focus on the arbitrary parameter t e which appears in the matrixĴ e,0 in Eq. (40) in the three flavor regime and must cancel in theÛ QED part of the evolution matrix. Let us denote the t e -dependent piece ofĴ e,0 witht e , where [t e ] 7, and the other components are zero. Using the matrixV it can be written ast e =Vt eV −1 , where [t e ] 10,1 = −t e and the other components are zero. Then, in the evolution matrix, the t e dependence takes the following form: All components except for (10, 1) of the parenthesis in Eq. (54) are zero trivially. The cancellation of the (10, 1) component then proceeds in the same way as in the QCD case: Therefore, the t e dependence ofÛ QED vanishes. The cancellation of the parameters t se1,2,...,14 in the second matrix product of Eq. (42) is more trivial. Let us define the second matrix product asVt seV −1 . In the evolution matrix, the matrixt se appears only in theÛ QCD-QED part and the cancellation can be understood in the following way: where we use the fact that (γ On the contrary, the cancellation of t s arising inÛ QCD-QED and t e inÛ QED-QED is highly non-trivial. The t s dependence, for example, resides inĴ s,0 ,Ĵ se,0 andĴ se,1 which appear in the matrixÛ QCD-QED . Logarithmic α s terms are accompanied bŷ J se,1 and by the matrix productsĴ sÛQED andÛ QEDĴs . Although we do not give an analytic explanation for these cancellations in this paper, we have checked that taking the sum of all terms in Eqs. (50) and (51) eliminates any t s and t e dependence ofÛ QCD-QED andÛ QED-QED . Now we have shown that the evolution matrix in Eq. (47) is independent of the undetermined parameters, so that we can set them to arbitrary values from the beginning. These parameters are directly related to the singular components in Eqs. (23), (24) of the standard solution in the literature. Therefore, our method automatically regularizes all singularities and these parameters correspond to the choices of the finite pieces of the regulated expressions, which can therefore be viewed as scheme parameters.
We have also found that the cancellation of the parameters occurs between the high and low scales. This insight is especially important when considering new physics at a high scale. The Wilson coefficients for a given model are typically calculated at leading order only. In the evolution to the scale of 1 GeV appropriate for Kaon physics one then usually neglects the corrections toK in Eq. (15) justified by the smallness of α s (µ 2 ) compared to α s (µ 1 ). In the typical applications in flavor physics, which do not involve corrections of order α EM , this procedure is scheme-independent. We here show that such a treatment is inconsistent in view of the cancellation of the singularity regulating scheme parameters.
This inconsistency does not appear in the QCD and QED parts which are nonsingular at f = 4, 5, 6. However the combined QCD-QED part, in which singularities persist for all numbers of flavors, will yield results depending on unphysical arbitrary scheme parameters if parts of the evolution matrix are discarded in the described way. Instead, the pieces ofK which depend on the scheme parameters t se must be consistently retained.

Validation of the Logarithmic Contribution
Finally, let us comment on the logarithmic contributionsĴ s,1 andĴ e,1 . At the O(α s ) part, we have the following logarithmic contributions to the evolution matrix, In theĴ s,1 matrix, the only nonzero component is [Ĵ s,1 ] 8,7 = −10/27. Using a calculation parallel to the one in the previous subsection, we find that the only nonzero component in the matrix productĴ s, Then, the (8, 7) component in the parenthesis in Eq. (60) becomes −(10/27)α 1 (α 2 /α 1 ) 1/9 ln(α 1 /α 2 ), and we arrive at Eq. (61). We find that this result is consistent with Eq. (40) of Ref. [25], where, in order to regulate the singularity, a small regulator is introduced in the eigenvalues ofγ (0)T s . With a similar calculation for the O(α EM /α s ) part we obtain the following term, This logarithmic contribution is also consistent with Eq. (2.28) of Ref. [13].
2.5 Higher orders in α EM and comparison with Ref. [22] The RG evolution in the pioneering papers [5,6,13] discards all terms which are quadratic or higher-order in α EM . Our solution in Eq. (47) is correct to order α 2 EM /α 2 s , but neglects terms of order α 2 EM /α s and higher. The extra term is numerically unimportant for the SM analysis, but matters in studies of new-physics contributions generated at very high scales, where α s is small. We come back to this point in Sec. 4.2. The RG evolution derived in Ref. [22] considers terms quadratic in α EM , including terms of order α 2 EM /α s which we neglect. In particular, the µ dependence of α EM affects the RG evolution at order α 2 EM /α 2 s and is therefore also included in Ref. [22]. While Ref. [22] addresses B decays, the derived formulae equally apply to K and were used in Ref. [14]. We argue that the inclusion of α 2 EM /α s terms in the RGE does not improve the prediction of K / K , because other terms of the same order are not included in the standard NLO solution: For instance, at this order the two-loop pure QED anomalous dimension matrixγ (1) e must be added toγ (g s (µ)) in Eq. (12).
Another issue are the ∆I = 1/2 operators which are generated by electroweak box diagrams, so that their Wilson coefficients are of order α EM . In agreement with Ref. [6] we find a small impact of these operators, contributing (−0.07 × 10 −4 ) to K / K . Furthermore, this contribution dominantly comes from A 2 which is entered by Q 11,12 through RG mixing triggered byγ (0) e and is thus O(α 2 EM /α s ) and to be discarded. While the contribution of Q 11,12 to A 0 is formally part of the NLO solution for K / K , it is numerically completely negligible (contributing −0.01 × 10 −4 ).
We close this section by comparing our solution of the RG equations in Eqs. (25)-(37) to the one in Ref. [22]. Actually, the latter also regulates all the singularities by logarithmic terms, and uses the diagonalisation ofγ In this form one can most easily compare our result with Eq. (47) of Ref. [22]. Thê [22], respectively. We have checked that our formulae of the RG evolution matrices are numerically equivalent to those in Ref. [22]. We find that our solution is easier to implement and leads to a faster numerical evaluation.

3
K / K in the Standard Model at Next-to-Leading Order In this section, we evaluate K / K in the Standard Model at next-to-leading order, using the evolution matrix derived in the previous section.
We calculate the Wilson coefficients v i and z i in Eqs. (9) and (10) with the methodology of Ref. [13]. Throughout this paper, the MS-NDR regularization scheme is used. For the next-to-leading order RG evolution of the Wilson coefficients, we use the singularity-free evolution matrix in Eq. (47) and systematically discard higher-order contributions. Table 1 shows our result of the Wilson coefficients at 18 GeV, and µ c = 1.4 GeV, which is the threshold scale between three and four flavor effective theories in Eqs. (9) and (10). Note that we include ln(m 2 c /µ 2 c ) contributions in the charm quark threshold correction z i (µ c ) in Eq. (10), where we use m c = 1.275 GeV [29]. To calculate α s (µ) we use RunDec:v1.0 with two-loop accuracy [30].  Next we take the hadronic matrix elements from a recent lattice QCD calculation [17][18][19][20], using the real parts (CP -conserving parts) of the isospin amplitudes A I=0,2 = (ππ) I=0,2 H |∆S|=1 eff K 0 as additional constraints [13]. These amplitudes have been measured very precisely [19], Since the real parts are dominated by Standard-Model tree-level coefficients z 2 (see Table 1), they can be used to fix one of the hadronic matrix elements (ππ) I |Q i (µ)| K 0 ≡ Q i (µ) I . Q 2 0 dominates the real part of A 0 , but contributes to the imaginary part only through the operator Fierz relations #1 Q 1 0 is the second largest contribution and the remaining matrix elements are almost negligible. The situation is more handy in the case of A 2 , where the real part is parameterized entirely by Q 2 2 due to the fact that Q 1 2 = Q 2 2 in pure QCD [6,13]. In our analysis we derive values of Q 2 0 and Q 2 2 at the scale µ from the experimental measurements of ReA 0 and ReA 2 , respectively #2 .
#2 On the other hand, once one introduces the ratio The decay amplitude of K → (ππ) I=0 has been computed using a 2 + 1 flavor lattice QCD simulation at the renormalization scale µ = 1.531 GeV [20]. In order to combine these matrix elements with the Wilson coefficients evaluated in the threeflavor regime -that is, at a scale below the charm quark mass-we need to evolve the hadronic matrix elements down to a scale below µ c . The isospin amplitude is given as where µ 1 < µ 2 and C i (µ) ≡ z i (µ) + τ y i (µ). In the final line, we use the fact that the physical amplitude A I is independent of the renormalization scale, so that In practice, we first evaluate the hadronic matrix elements for the I = 0 states at µ = 1.3 GeV from the lattice results [20] using a three flavor evolution matrix, cf. Eq. (69). Here we use α  [20]. Then we determine Q 2 (µ) 0 (and Q 4,10 (µ) 0 through Eq. (66)) from the experimental value of ReA 0 using the Wilson coefficients shown in Table 1. We have taken the CKM parameters from CKMfitter [31]. The results are shown in Table 2(a).
The decay amplitude of K → (ππ) I=2 has also been computed using a 2 + 1 flavor lattice QCD simulations, albeit at the scale µ = 3.0 GeV [17][18][19]. According to Ref. [18], one can extract the lattice results in an operator basis renormalized by the MS-NDR regularization scheme. From Ref. [19], which is the latest lattice QCD calculation for I = 2, we obtain where the results of the (q /, q /) intermediate scheme are taken as central value, while the results of the (γ µ , γ µ ) scheme are taken as uncertainty. Using the three flavor evolution one can calculate ImA I /ReA I without using the fit of Q 2 I to the data. Ref. [14] uses this strategy with the parameter range 0 ≤ q ≤ 0.1. Basically, the difference with our method (corresponding to the q-dependent terms in Ref. [14]) only affects numerically subleading contributions (the i = 3, 4, 9, 10 components of ImA 0 /ReA 0 ). In either method the hadronic uncertainties are reduced compared to the choice to take all matrix elements from lattice. Table 2. The hadronic matrix elements (a), (b) and B parameters (c) extracted from the lattice calculations for I = 0 [20] and I = 2 [19]. The experimental values of the real parts of the amplitudes have been used [19]. The large errors result from the quoted lattice errors on the hadronic matrix elements. The experimental errors are small in comparison. We take µ = 1.3 GeV. s (3 GeV) = 0.24544 [18]. Then, from the experimental value of ReA 2 we determine Q 2 (µ) 2 (and Q 1,9,10 (µ) 2 through Eq. (66), Q 1 (µ) 2 = Q 2 (µ) 2 and Q 9 = 1 2 (3Q 1 − Q 3 ) which is a Fierz relation). The results are shown in Table 2(b). Note that through the evolution matricesÛ QED andÛ QCD-QED this procedure generates small nonzero values of Q 3-6 (µ) 2 , which are regarded as non-electroweak penguin contributions to ImA 2 . Since the lattice simulations have not calculated them at 3.0 GeV, one should not use them at the lower hadronic scale µ. On the other hand, they have been calculated with chiral perturbation theory [15,16] and are included in the isospin-violating correctionsΩ eff of Eq. (1) #3 . Therefore, we have decided to omit these contributions at the hadronic scale µ.
To compare with the literature, we also extract B parameters from the hadronic matrix elements in Table 2(c). These B parameters are defined as in Ref. [14]: All other B parameters are defined in Ref. [13]. For running quark masses, we use the lattice results m s (2 GeV) = 93.8(2.4) MeV and m d (2 GeV) = 4.68 (16) MeV with the three-flavor RG evolution [32]. Since the uncertainty from the strange quark mass is already included in the lattice results of Q i I as one of the systematic errors, we do not include it in the estimation of uncertainties of the B parameters. The B parameters are consistent with Ref. [14], and we also confirmed the almost µ-independent behavior of B (1/2) 6 (µ) and B (3/2) 8 (µ) [13]. Note that in the following analysis we will directly use the hadronic matrix elements Q i I rather than the B parameters.
Finally we combine the short-distance and long-distance contributions. The master equation of K / K is given in Eq. (1). Since the isospin-violating correction by the electroweak penguins to ImA 0 are already subtracted fromΩ eff as Q 7-10 0 , one should evaluate the last term in Eq. (1) as with the two terms representing the contributions from Q 3-6 0 and Q 7-10 0 , respectively [14]. In addition, the experimental values of ReA 0 in Eq. (64) and | K | = 2.228 × 10 −3 [29] are used. Our result for K / K in the Standard Model at the nextto-leading order is The first error originates from the lattice-QCD simulations [19,20] and is dominated by the uncertainty stemming from Q 6 0 (which is ±4.52×10 −4 ) (see Figure 2(c)). The uncertainties from Q 3 0 through Eq. (66) and from Q 8 2 are subleading (±0.77×10 −4 and ±0.56 × 10 −4 , respectively). The second uncertainty comes from perturbative higher-order corrections, which we estimate in two ways. Firstly, we estimate uncertainties from higher-order corrections to the Wilson coefficients by calculating the RG evolution of the Wilson coefficients with a different method. Instead of using the analytic evolution matrices formulated in Sec. 2, we solve the corresponding set of differential equations numerically.
Since this RG evolution contains higher-order (namely O(α 2 s , α s α EM )) corrections, the result is interpreted as a conservative estimate of the uncertainty in the short-distance contributions. As a result, we find that the Wilson coefficients are shifted by about 10 percent compared with Table 1, and we obtain K / K = −0.32 × 10 −4 . Hence, we estimate that the uncertainty from higher-order corrections is ±1.38×10 −4 . Secondly, we have investigated the µ c and µ dependences of K / K . In Fig. 1(a), we show the µ c dependence of K / K in the range 1.3 < µ c < 3.0 GeV with fixed µ = 1.3 GeV. In Fig. 1(b), we vary µ with µ c fixed at 1.4 GeV. We find that the µ dependence is small, ±0.77 × 10 −4 , while the µ c dependence is slightly larger, ±1.09 × 10 −4 . The scale µ enters the prediction in three ways: First, the decomposition of the isospinviolating corrections in Eq. (75) is imposed at this scale. Second, the omitted nonelectroweak penguin contributions to ImA 2 depend on µ, and third, the experimental values of ReA 0 and ReA 2 to fix Q 2 (µ) 2 and Q 2 (µ) 0 are imposed at the hadronic scale µ. In this process, we double-count the uncertainty from the isospin-violating contributions, however, we find that these uncertainties are very small compared with the uncertainties stemming from lattice and thus we have not investigated them any further. We show the µ dependences of ImA 0 (and not the µ dependence of (1 −Ω eff )ImA 0 ) and ImA 2 in Figs. 1(c) and (d), respectively. We add the three uncertainties in quadrature. Strictly speaking, this double-counts some pieces of the unknown higher-order corrections.
The third uncertainty in Eq. (76) stems from isospin-violating corrections [15,16], such as strong isospin violation (m u = m d ), non-electroweak penguin transitions in the I = 2 state and ∆I = 5/2 corrections [33,34]. The uncertainty is dominated by the non-electroweak penguin contributions to ImA 2 , however, the uncertainty in The last uncertainty in Eq. (76) comes from the running mass of the top quark m t (m t ) = 163.3 ± 2.7 GeV [35]. Since the other uncertainties we have not elaborated here are negligibly small according to Ref. [14], we have omitted them in our error estimate.
which is consistent with Refs. [14] and [20]. On the other hand, it is well-known that the experimental value is much larger [36][37][38][39][40][41]. The current world average is [29], We observe that our prediction of K / K in the Standard Model is 2.8 σ below the experimental value. This small Standard Model prediction and thus the large tension is supported by the large-N c "dual QCD" approach [42][43][44][45][46][47], which is an entirely different approach to low energy QCD than lattice gauge theory. There has been a dispute concerning the role of final-state interactions (FSI) for the size of Q 6 0 , with y 6 Q 6 y 4 Q 4 y 9 Q 9 y Q y 5 Q 5 3 Q 3 y 10 Q 10 y 7 Q 7 y the chiral perturbation community favouring an enhancement of Q 6 0 by FSI [48] and an opposing view of the large-N c community [47]. Modern lattice calculations do include FSI [49] and will speak the final word on FSI. Since the main uncertainty of the SM prediction for K / K comes from statistical and systematical errors in the lattice calculation of the hadronic matrix elements for A 0 , the expected progress in this field will sharpen the Standard Model prediction in the near future [20].
We note that in absence of a lattice result for the hadronic matrix element and the smallness of the corresponding Wilson coefficient, we omit the contribution from the chromomagnetic penguin operators Q 8g = m s g s /(16π 2 )sT a σ µν (1 − γ 5 )dG µν a (and the opposite-chirality analogueQ 8g ). According to Ref. [14], chromomagnetic penguins contribute |0.2-0.7| × 10 −4#4 to K / K , which rather small compared with the QCD-penguin and QED-penguin contributions (see Figure 2(c)). Even if we add this contribution as +0.7 × 10 −4 to the central value (to the higher-order uncertainty) of K / K , the discrepancy still persists at 2.7 (2.8) σ.
In Fig. 2 we show the composition of ImA 0 , ImA 2 and K / K with respect to the operator basis. We observe that the positive dominant contribution to K / K comes from Q 6 while Q 9 is subdominant. The dominant negative contribution comes from Q 8 while Q 4 is subdominant. Remarkably, their sum almost cancels at next-to-leading order. This leads to an extremely small central value of the Standard Model prediction for K / K .
Although the results of the Wilson coefficients by themselves are slightly different #4 The sign depends on the sign of the hadronic matrix element. The preliminary lattice calculation of π|Q 8g |K [50] and calculations in the chiral quark model [51][52][53] imply that a contribution to K / K is positive at the leading order. However, next-to-leading order contributions to (ππ) I=0 |Q 8g |K 0 are expected to mess up the leading order estimate because of a parametric enhancement ∝ 1/N c · m 2 K /m 2 π [54,55].
when compared to the result of Ref. [14], the products with the hadronic matrix elements are well consistent #5 . The main difference between this reference and our analysis is in the subleading contributions. In Ref. [14], the hadronic matrix elements Q 3 (µ) 0 , Q 5 (µ) 0 and Q 7 (µ) 0 are set to be 0 as central values, while we have evaluated them from the lattice data. The numerical difference in K / K is ∼ −1 × 10 −4 . We also find that the contribution of O(α 2 EM /α 2 s ) terms, which has not been considered in the literature so far, only contributes to K / K as little as −0.10 × 10 −4 . This term, however, can be relevant in new-physics models with TeV-scale isospin violation.

Preliminaries
Upon integrating out heavy degrees of freedom in models of new physics, new contributions to Wilson coefficients of the Standard Model operators Q i (and their oppositechirality analoguesQ i ) arise.
As we have shown in the previous section, the Standard Model prediction of K / K is significantly below the experimental data. Although the discrepancy is only 2.8 σ at present, its confirmation with higher significance by future lattice results may establish a footprint of new physics. Indeed, several new physics models can alleviate the K / K tension, like generic flavor-violating Z and Z models [56][57][58], 331 models [59][60][61], the Littlest Higgs model with T -parity [62], flavor-violating additional pseudo-scalar models [63], and the Minimal Supersymmetric Standard Model [64,65].
Since K / K is linear in the Wilson coefficients, the SM and new-physics contributions are simply additive: (80) Using the following effective Hamiltonian for the new physics contributions, where the opposite-chirality operatorsQ i are found from Q i by interchanging V −A ↔ V + A, the new physics contribution is given by Indeed, the values of y 6 Q 6 0 and y 8 Q 8 2 are in good agreement with Ref. [14].
where the isospin-violating correction in Eq. (75) is and we employed Q i (µ) The evolution matrix in Eq. (82) is given bŷ Since the matching matrices depend only on the difference of the number of active upand down-type quark flavors, we takeM t (m) =M c (m). Note that the RG evolution of the opposite-chirality operators is the same as for the Standard Model operators and that these two sets of operators do not mix with each other. We also note that the chromomagnetic operators are omitted in our analysis.
In this section, we give a useful formula for the new physics contributions to K / K considering the analytic solutions of the next-to-leading order evolutions matrices and the hadronic matrix elements we derived. We note that we omit the weak boson exchanges in the RG evolutions from µ NP to M W , where µ NP represents the matching scale between the new physics and the effective Hamiltonian in Eq. (81). Like the photon exchanges one should treat weak boson exchanges as next-to-leading contributions. Note that large isospin violation in new-physics models enters K / K through the initial conditions of the Wilson coefficients and not through the RG evolution.
We also should comment on the running of α EM . Above M W scale, we use e(µ NP ) = g(µ NP )g (µ NP )/ g 2 (µ NP ) + g 2 (µ NP ), and β e 0 = β g 0 / cos θ 2 W (M Z ), where β g 0 = −53/9 (µ < m t ) or −41/6 (µ > m t ). Strictly speaking, we have to consider the running of θ W for consistency. However, we have checked that the numerical effect for an O(10 TeV) scale of new physics is small. Therefore we use a fixed value: sin 2 θ W = 0.231.

Counting of Orders
In a full next-to-leading order estimation, we have to consider the leading order term O(1) arising from the one-loop QCD RG evolution as well as the terms defined as next-to-leading order, which are: the one-loop QED correction O(α EM /α s ), the QCD two-loop correction O(α s ), and the two-loop term including a photon and a gluon at O(α EM ). The next-to-leading order RG evolution matrix has an additional O(α 2 EM /α 2 s ) correction, which appears only at this order. Hereafter, we will always refer to these orders when labelling perturbative quantities of the Wilson coefficients and the evolution matrices as s 0 , s e , s s , s se , s ee andÛ 0 ,Û e ,Û s ,Û se ,Û ee , respectively.
When we multiply two quantities which are given by a perturbation series, we have to carefully keep track of and consistently discard higher orders of the perturbative series. This is a subtle and cumbersome feature which complicates mathematical expressions. In this context, equations of the RG evolution should be more of a symbolic character which are exact in the limit of expanding the corresponding quantities to all orders. Since we necessarily truncate the perturbation expansion of the Wilson coefficients as well as the evolution matrices at some point, a product of them at next-to-leading order is represented as follows: +Û e s e +Û ee s 0 +Û 0 s ee Here we have suppressed the opposite-chirality coefficients s and the arguments of U (µ, µ NP ) and s(µ NP ) for better readability. This procedure defines s NLO (µ) as a next-to-leading order quantity, where higher orders have been discarded consistently. In view of undetermined Wilson coefficients, it is beneficial to arrange the terms above according to the Wilson coefficients evaluated at the new physics scale as where we have again suppressed s and the arguments ofÛ (µ, µ NP ) and s(µ NP ). For given numerical values for the hadronic matrix elements at a low scale and with our evolution matrices connecting µ NP with the low scale µ, we can determine the weights which multiply the Wilson coefficients Im[ s(µ NP ) − s(µ NP )] in Eq. (87) for any chosen scale of new physics.

Evolution Matrices at the TeV scale
Above the electroweak scale we observe an approximately logarithmic behavior of the evolution matrixÛ (µ, µ NP ) in Eq. (85) with increasing energy scale. This observation allows us to derive an approximation for the evolution matrix in the high energy region, which has an error of only a few percent. We give approximate functions for all components of the evolution matrix linking the new physics scale to the hadronic scale. Cast in the form U 0,e,s,se,ee (µ, we combine them in terms of Eq. (87). Using the analytic evolution matrices evaluated in Sec. 2 and the next-to-leading order matching matricesM c,b,t , we obtain U 0 (µ, µ NP ) +Û e (µ, µ NP ) +Û s (µ, µ NP ) +Û se (µ, µ NP ) +Û ee (µ, µ NP ) for the O(1) Wilson coefficients at the µ NP scale, and for the O(α EM /α s ), O(α s ), O(α EM ) (or O(α 2 EM /α 2 s )) Wilson coefficients at the µ NP scale, respectively. Here µ = 1.3 GeV and µ c = 1.4 GeV are taken, and the fitting matricesÛ f it are given in Appendix B. We find that these approximate evolution matrices are highly accurate in the range of 500 GeV-10 TeV.
In order to estimate which Wilson coefficients are expected to gain large enhancements through the RG evolution, we calculate weights for the Wilson coefficients at the µ NP scale. We regard the coefficients of Q K (µ) T iÛ i (µ, µ NP ) ( s(µ NP ) − s(µ NP )) in Eq. (87) as weights of the Wilson coefficients.
In Table 3, we list the coefficient Q K (µ) T (Û 0 +Û e +Û s +Û se +Û ee ) for the O(1) Wilson coefficients at the scale µ NP = 1, 3, 5 and 10 TeV in units of (GeV) 3 , where the hadronic matrix elements of Table 2 Tables 4, 5, and 6, respectively. Note that these values are not obtained by fitting but using the exact analytic evolution matrices. We observe that these values are of course dominated byÛ 0 , with the sub-dominant contribution stemming fromÛ e because of the 1/ω + enhancement andÛ s . We also find, that the largest weights come in the 7 and 8 components, and they are further enhanced through the RG evolution in the high energy regime. Compared with the coefficients at the weak scale,   Tab. 6), which shows the impact of the NLO corrections on these elements. Although the enhancement factor from the RG evolution has been pointed out before in Ref. [58,60] within a leadingorder analysis, it has not been considered in most of the literature. We emphasize that this factor should be included when one studies TeV-scale new-physics contributions to the QED-penguin operators in order to alleviate the K / K discrepancy.

Conclusions and Discussion
Based on the first complete lattice calculation of the hadronic matrix elements for the K → ππ decay, we have evaluated the Standard-Model prediction of K / K at the next-to-leading order. It is well known that the analytic RG evolution matrices for the ∆S = 1 nonleptonic effective Hamiltonian at the next-to-leading order contains singularities in intermediate steps of the calculation. These singularities make practical calculation laborious even though appropriate regulators disappear from the final (physical) result. In this paper, we have generalized the analytic ansatz of the Table 3. The coefficient Q K (µ) T (Û 0 +Û e +Û s +Û se +Û ee ) for the O(1) Wilson coefficients at the scale µ NP in units of (GeV) 3 , where µ = 1.3 GeV.  Roma group [11,21] to solve the RG equations and derive a singularity-free solution by adding logarithmic terms to the ansatz. As a novel feature of our solution compared to Refs. [22,25] we do neither require the diagonalization of the LO anomalous dimension matrix nor case-by-case implementations for different eigenvalues of this matrix. Instead, the different cases are encoded in theĴ matrices given in Eqs. (38)- (46) and Appendix A. The singular nature of the RG equations leads to the presence of spurious parameters which cancel between the high-scale and low-scale NLO terms in the RG evolution matrix and thereby do not produce any ambiguity and play the role of scheme parameters with respect to the regularization of the singularities. Thus we have explicitly proven that all singularities are automatically treated in the proper way without the need for a manual regularization of the evolution matrix. This feature also leads to a subtlety whenever the NLO evolution matrix is combined with LO initial conditions for the Wilson coefficients, as one usually does in studies of new-physics contributions to K . Using the improved RG evolution matrices and applying the recent lattice results, we have calculated K / K in the Standard Model at the next-to-leading order. Our final results is K / K = (1.06 ± 5.07) × 10 −4 , which is 2.8 σ below the measured value. Our result is consistent with the recent literature and highlights a tension between the Standard-Model prediction and experiment. The uncertainty is dominated by the lattice result of (ππ) 0 |Q 6 |K 0 . Therefore, upcoming improvements of lattice calculations will reveal whether this tension really calls for new physics or not.
We have also evaluated the evolution matrices in the high energy region for calculations of new physics contributions to K / K . To this end we have further obtained an easy-to-use approximate formula for the RG evolution matrices in the TeV region at the next-to-leading order and have also calculated the weights for each of the Wilson coefficients at the scale of new physics. We observe that the largest weights come in the 7 and 8 components of the Wilson coefficients and that they are further enhanced through the RG evolution between electroweak and TeV scales. Here we confirm the feature noticed at LO in Ref. [58,60] and find a further enhancement by the NLO corrections to the evolution matrices. Especially the Wilson coefficients of the QED-penguin operators at the scale of 1-10 TeV increase by 50-100 % compared with the Wilson coefficients at the weak scale.
acknowledges support from the DFG-funded doctoral school KSETA.

A Solutions for the matricesĴ
In this appendix, we summarize the solutions for the matricesĴ of Eqs. (29)- (37).
Here we set all arbitrary parameters to be zero, which does not affect the evolution matrix in Eq. (47). We find that the matricesĴ s,1 ,Ĵ e,1 ,Ĵ se,2 andĴ ee,1 are zero matrices in the case where the active number of flavours is four, five or six.
In the case of four active quark and three active lepton flavors, the matricesĴ are given as follows: In the case of five active flavours, the matricesĴ are given as follows: Above the scale M W in the f = 5 case onlyĴ ee,0 is replaced bŷ In the case of six active flavours, the matricesĴ are given as follows:

B Approximation of Evolution Matrices
In this appendix we list the approximate evolution matricesÛ f it of Eqs.