On the importance of NNLO QCD and isospin-breaking corrections in ε′/ε

Following the 1999 analysis of Gambino, Haisch and one of us, we stress that all the recent NLO analyses of ε′/ε in the Standard Model (SM) suffer from the renormalization scheme dependence present in the electroweak penguin contributions as well as from scale uncertainties in them related to the matching scale μW and in particular to μt in mt (μt ). We also reemphasize the important role of isospin-breaking and QED effects in the evaluation of ε′/ε. Omitting all these effects, as done in the 2015 analysis by RBC-UKQCD collaboration, and choosing as an example the QCD penguin (Q6) and electroweak penguin (Q8) parameters B 6 and B (3/2) 8 to be B (1/2) 6 = 0.80 ± 0.08 and B 8 = 0.76 ± 0.04 at μ = mc = 1.3 GeV, we find (ε′/ε)SM = (9.4 ± 3.5) × 10−4, whereas including them results in (ε′/ε)SM = (5.6 ± 2.4) × 10−4. This is an example of an anomaly at the 3.3 σ level, which would be missed without these corrections. NNLO QCD contributions to QCD penguins are expected to further enhance this anomaly. We provide a table for ε′/ε for different values of B 6 and the isospin-breaking parameter ̂ eff, that should facilitate monitoring the values of ε′/ε in the SM when the RBC-UKQCD calculations of hadronic matrix elements including isospinbreaking corrections and QED effects will improve with time.

= 0.76 ± 0.04 at μ = m c = 1.3 GeV, we find (ε /ε) SM = (9.4 ± 3.5) × 10 −4 , whereas including them results in (ε /ε) SM = (5.6 ± 2.4) × 10 −4 . This is an example of an anomaly at the 3.3 σ level, which would be missed without these corrections. NNLO QCD contributions to QCD penguins are expected to further enhance this anomaly. We provide a table for ε /ε for different values of B (1/2) 6 and the isospin-breaking parameter eff , that should facilitate monitoring the values of ε /ε in the SM when the RBC-UKQCD calculations of hadronic matrix elements including isospinbreaking corrections and QED effects will improve with time.

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 [1]. In the SM ε /ε is governed by QCD penguins (QCDP) but receives also an important contribution a e-mail: jason.aebischer@tum.de b e-mail: christoph.bobeth@tum.de c e-mail: aburas@ph.tum.de from the electroweak penguins (EWP), pointed out already in 1989 [2,3], that entering ε /ε with the opposite sign to QCDP suppress this ratio significantly. The partial cancellation of these two contributions in addition to the significant uncertainties in the evaluation of the hadronic matrix elements of QCDP and EWP operators is the reason why until today a precise prediction for ε /ε in the SM is still missing.
The situation of ε /ε in the SM by the end of 2017 could be briefly summarized as follows: • The analysis of ε /ε by the RBC-UKQCD lattice QCD (LQCD) collaboration based on their 2015 results for K → ππ matrix elements [4,5], as well as the analyses performed in [6,7] that are based on the same matrix elements but also include isospin breaking effects [8,9], found ε /ε in the ballpark of (1−2)×10 −4 . This is by one order of magnitude below the experimental world average from NA48 [10] and KTeV [11,12] collaborations, However, with an error in the ballpark of 5 × 10 −4 obtained in these analyses, one could talk about an ε /ε anomaly of at most 3 σ . Simultaneously, we note that the 2015 RBC-UKQCD result for the ππ-strong-interaction phase δ 0 of the isospin I = 0 amplitude is in almost 3 σ conflict with the result from extrapolations in the chiral limit [13]. This suggests that there were methodical problems with the 2015 RBC-UKQCD calculation, which were meanwhile successfully addressed, as will be reported later. As a conclusion, one has to be aware that for I = 0 the K → ππ matrix elements, represented mainly by the parameter B (1/2) 6 , and hence also the 2015 RBC-UKQCD result for ε /ε suffer from an unaccounted systematic uncertainty.
• An independent analysis based on hadronic matrix elements from the Dual QCD (DQCD) approach [14,15] gave a strong support to the 2015 RBC-UKQCD result and moreover provided an upper bound on ε /ε in the ballpark of 6×10 −4 . However, in this approach the treatment of ππ strong interaction phases is presently problematic. • Chiral perturbation theory (ChPT) together with large-N considerations 1 are used in [17], leading to a SM prediction of ε /ε = (15 ± 7) × 10 −4 . The uncertainties are larger than in [6,7], reflecting in part the difficulties in matching long-distance and short-distance contributions in this framework, but are also of parametric origin due to low-energy constants. Consequently, the predicted central value is one order of magnitude larger compared to DQCD and lattice results of 2015, but with a small tension of 1.6 σ in view of the large uncertainties.
Recently progress towards an improved estimate of ε /ε in the SM has been made: • The RBC-UKQCD collaboration is expected to present this year new values of the K → ππ hadronic matrix elements, most importantly the parameter B (1/2) 6 . In particular the discrepancy in the prediction of the ππ-stronginteraction phase δ 0 has been identified [18][19][20][21] in the form of excited-state contamination requiring the introduction of additional ππ operators in the simultaneous fits to allow for a better isolation of the ground state. It can be expected that the statistical errors will decrease, though less dramatically as assumed before due to the additional operators. Unfortunately, the inclusion of isospin-breaking and QED effects will still take more time.
• An improved estimate of isospin-breaking corrections to ε /ε has been presented in [22] increasing moderately the role of these corrections in suppressing ε /ε. The updated ChPT analysis [22] resulted in the value in full agreement with the experimental world average (1). • The preliminary result on NNLO QCD corrections to QCDP contributions [23,24] demonstrates significant reduction of various scale uncertainties, foremost of μ c , and indicates an additional modest suppression of ε /ε.
In contrast to the expected RBC-UKQCD result, the ChPT analysis includes isospin-breaking and QED corrections but the known difficulties in matching long-distance and shortdistance contributions in this approach imply a large uncer-tainty. In particular, the absence of the so-called meson evolution in ChPT that suppresses ε /ε within the DQCD approach [15,16] is responsible for the poor matching and according to the latter authors responsible in part for the large value of ε /ε in (2). The DQCD analysis in [25] demonstrates through the example of BSM matrix elements in K 0 − K 0 mixing that the effects of meson evolution are included in the present LQCD calculations. As shown in [25], neglecting this evolution in the case of K 0 − K 0 mixing would miss the values of the relevant hadronic matrix elements by factors of 2 − 4, totally misrepresenting their values obtained by three LQCD collaborations [26][27][28][29][30]. The fact that in K 0 − K 0 mixing the FSI are absent allows to study the impact of meson evolution better than it is possible in K → ππ decays. Yet, as demonstrated in [14,15] these important QCD dynamics must also be present in K → ππ and is also required by the proper matching of long-distance and short-distance contributions. Therefore it is expected to suppress the result for ε /ε in (2). Now all the 2015 analyses of ε /ε and the one in [22] used the known Wilson coefficients at the NLO level [31][32][33][34][35][36] in the naive dimensional regularization (NDR) scheme [37]. But already in [38] and recently in [16,39] it has been pointed out that without NNLO QCD corrections to the EWP contribution the results for ε /ε are renormalization-scheme dependent and exhibit significant non-physical dependences on the scale μ t at which the top-quark mass m t (μ t ) is evaluated as well as on the matching scale μ W .
Fortunately, all these uncertainties have been significantly reduced in the NNLO matching at the electroweak scale performed in [38] and it is of interest to look at them again in the context of new analyses with the goal to improve the present estimate of ε /ε. Additional importance in such an analysis is the finding in [38] that these corrections further suppress ε /ε relative to the NLO results performed in the NDR scheme.
In view of the fact that LQCD calculations contain both the meson evolution and FSI, the optimal strategy for the evaluation of ε /ε as of 2019 appears to be as follows: , respectively-with improved values of ππ-stronginteraction phases δ 0,2 -but determine hadronic matrix elements of (V − A)⊗(V − A) operators from the experimental data on the real parts of the K → ππ amplitudes as done in [6,34]. 2. Use the result for isospin-breaking and QED corrections from [22], which are compatible with the ones obtained already 30 years ago in [40]. 3. Use the NNLO QCD contributions to EWP in [38] in order to reduce the unphysical renormalization scheme and scale dependences. 1 4. Include NNLO QCD contributions to QCDP from [23,24] in order to reduce left-over renormalization scale uncertainties.
In view of the fact that meson evolution and the remaining three effects tend to suppress ε /ε, whereas the aforementioned systematic uncertainties in the 2015 RBC-UKQCD calculation of B (1/2) 6 related to FSI effects discussed before in connection with δ 0 could tend to increase it, as does FSI in the case of ChPT, it could well happen that future LQCD predictions of ε /ε in the SM increase only moderately to end up in or a bit above the ballpark ε /ε ≈ (5 ± 2) × 10 −4 of the expectation based on the DQCD approach in [16].
The main goal of our paper is to illustrate how a future result from RBC-UKQCD would be affected by the inclusion of known isospin-breaking and QED corrections from [22] in point 2. and the NNLO QCD contributions to EWP in [38] in point 3. We also comment on the expected size of NNLO QCD contributions to QCDP from [23,24] in point 4. leaving a detailed analysis of them to these authors.
Our paper is organized as follows. In Sect. 2 we recall a number of basic formulae that will be used in the rest of our paper. In Sect. 3 we address the issue of scale and renormalization scheme dependences at the NLO level resulting from the EWP sector. We illustrate the size of these unphysical effects present at the NLO level that would increase the errors quoted in the existing NLO analyses but are significantly reduced when NNLO QCD corrections to EWP contributions from [38] are taken into account. To this end we use first as an example particular values for the hadronic parameters B (1/2) 6 and B (3/2) 8 quoted in the abstract of our paper and quantify also the role of isospin-breaking and QED corrections from [22]. As the values of B (1/2) 6 and the size of isospin-breaking corrections are expected to dominate the theoretical uncertainties in ε /ε for some time, we present in Sect. 4 a table of the SM values of ε /ε for different B (1/2) 6 and the isospin-breaking parameter eff that should facilitate monitoring the SM estimates of ε /ε when the LQCD calculations of hadronic matrix elements including isospinbreaking corrections and QED effects will improve with time. A brief summary and an outlook are given in Sect. 5.

An analytic formula
As in [6], our starting expression is the formula where [9,22] Here a and eff summarize isospin-breaking corrections and include strong isospin violation (m u = m d ), the correction to the isospin limit coming from I = 5/2 transitions and electromagnetic corrections as first summarized in [8,9] and recently updated in [22]. Our eff , defined by differs from eff in [8,9,22] as in contrast to these papers it does not include EWP contributions to Im A 0 , summarized in these papers by 0 . This is indicated here by 0 | α=0 , which contains the remaining contributions only. We find it more natural to calculate Im A 0 including both QCD and EWP contributions as this allows to keep track of NP contributions to Im A 0 . The dominant EWP contribution to ε /ε is of course present in Im A 2 . In fact the RBC-UKQCD collaboration includes EWP contributions to Im A 0 as well. We note also that the latest central value for IB = 0.25 ± 0.08 from [22] agrees perfectly with the one obtained already 30 years ago in [40]. The real parts of the isospin amplitudes A 0,2 in (3) are then extracted from the branching ratios on K → ππ decays in the isospin limit. In the limit a = 1 and eff = 0 the formula in (3) reduces to the one used in RBC-UKQCD [4,5], where all isospin-breaking corrections except for EWP contributions at the NLO level have been set to zero.
Using the technology in [6] we arrive at the formula with the coefficients a (1/2) 0,6 and a (3/2) 0,8 given in Table 1 at NLO and NNLO from EWPs as discussed below. Explicit formulae for a (1/2) 0,6 and a (3/2) 0,8 in terms of Wilson coefficients and Re A 0,2 are given in [6]. The values of the Wilson coefficients used by us are collected in Appendix A, whereas λ t = V td V * ts is the relevant CKM combination.
As an example we will first use the values to be compared with the 2015 values B While we do not expect significant modification of the RBC-UKQCD result for B from ChPT into account [22], we allow for an enhancement of B , which however is still consistent with the arguments in [15] that the suppression of B (1/2) 6 by meson evolution below unity is stronger than its enhancement by FSI. We emphasize that the choice of B (1/2) 6 > 1.0, in the spirit of [22], are considered. We anticipate a significant reduction of the error on B (1/2) 6 in the new results of RBC-UKQCD collaboration relative to its 2015 analysis so that the expectations from [15] and [22] will be tested.

Scale uncertainties at NLO
It should be emphasized that although the NLO QCD analyses of ε /ε in [31][32][33][34][35][36] reduced renormalization scheme dependence in the QCDP sector, the dependence of ε /ε on the choice of μ t in m t (μ t ) remained. This dependence can only be removed through the NNLO QCD calculations, but in the QCDP sector it is already weak at the NLO level because of the weak dependence of the QCDP contributions on m t . On the other hand, as pointed out already in [38], the EWP contributions at the NLO level suffer from a number of unphysical dependences.
• First of all there is the renormalization-scheme dependence with ε /ε in the HV scheme, as used in [35,36], generally smaller than in the NDR scheme used in [31][32][33][34]. In what follows we will consider only the NDR scheme as this is the scheme used by the RBC-UKQCD collaboration and other analyses listed above. • The dependence on μ t , which is much larger than in the QCDP sector because the EWP contributions exhibit much stronger dependence on m t . Increasing μ t makes the value of m t smaller, decreasing the EWP contribution and thereby making ε /ε larger. At NLO there is no QCD correction that could cancel this effect. • The dependence on the choice of the matching scale μ W .
It turns out that with increasing μ W in the EWP contribution, the value of ε /ε decreases.
One should note that the scales μ W and μ t can be chosen to be equal or different and they could be varied independently in the ranges illustrated in Fig. 1 implying significant uncertainties in the NLO prediction for ε /ε as demonstrated in [38]. In obtaining the values in Table 1 we provide the two settings from [38]: i) μ W = μ t = m W as well as ii) μ W = m W and μ t = m t . For example ii) has been used in [6]. Other choices of these scales would significantly change the NLO values of ε /ε with significantly reduced change when NNLO corrections to EWPs are included.
We next evaluate ε /ε for the values of B 1 that is a value by a factor of 7 larger than the 2015 result from the RBC-UKQCD collaboration. The quoted error is a guess estimate based on the uncertainties in (12) and scale uncertainties as well as omission of isospin-breaking effects ignoring the known signs of these effects. But as we will see soon its precise size is irrelevant for the point we want to make. The result in (13) is compatible with experiment (1) with a tension of 1.7 σ . At first sight it would appear that this result confirms the claims in [17] and [22] as (13) is quite consistent with (2). But such a conclusion would be false as we will illustrate now.
Indeed as stated above at the NLO level significant dependences on μ W and μ t are present and the impact of a nonvanishing eff is very significant. In order to exhibit these dependences we vary in Fig. 1  Fortunately all these uncertainties have been significantly reduced in the NNLO matching at the electroweak scale performed in [38]. In the NDR scheme, used in all recent analyses, these corrections enhance for (i) μ W = μ t = m W the EWP contribution by roughly 7 % and for (ii) μ W = m W and μ t = m t by 16 %. Thereby they imply a negative shift in ε /ε that depends on B which compared with the experimental value in (1) signals an anomaly at the level of 3.3 σ . In Table 2 below, we have set Im λ t = 1.4 × 10 −4 . For the result in (13) and (14) we have used Im λ t = (1.43 ± 0.04) × 10 −4 , based on recent analyses of the unitarity triangle by the bayesians ("UTfitter") and frequentists ("CKMfitter") that can be found in [43] and [44], respectively. The error budget, discussed in Appendix B and summarized in Table 6, would imply the parametric theoretical error of 2.3 × 10 −4 . We increased it in order to take into account left-over scale uncertainties both in the EWP sector discussed here and in the QCDP sector analyzed at NNLO in [23,24]. But one should keep in mind that the central value in (14) will be shifted down by NNLO QCD corrections to QCDP by about 0.5 × 10 −4 as indicated in the preliminary plots in [23,24] without modifying the error in (14). We are looking forward to the final results of these authors. Our NNLO central value in (14), represented in Fig. 1 by the black points at μ t = {m W , m t , 300 GeV}, is much less dependent on μ t . This exercise shows also the importance of isospin-breaking corrections. They are significantly larger than the NNLO QCD corrections to EWP contributions.
It should be emphasized that in [38] complete O(α W α s ) and O(α W α s sin 2 θ W m 2 t ) corrections, with α W = α/ sin 2 θ W , to the Wilson coefficients C 7−10 (μ) of EWP operators at μ = m c have been calculated. In particular as demonstrated in Sect. 3 of that paper no three-loop anomalous dimensions of involved operators are necessary to find these corrections. See formula (3.14) of that paper. In order to complete the NNLO analysis of EW contributions one should calculate m tindependent O(α W α s sin 2 θ W ) corrections, which as argued in [38] are much smaller than the ones included here.
Much more difficult is the NNLO analysis of QCD penguin contributions which in addition to two-loop calculations requires three-loop anomalous dimension matrices [23,24], obtained fortunately already in [45].
Inspection of the formulae in [6] together with the numbers in Table 1 shows that the NNLO matching corrections lead mainly to an enhancement of the coefficient a 07 (1.16) due to y 8 (μ), whereas the NNLO impact on ε /ε through y 7,9,10 is negligible due to the smaller matrix elements multiplying them. The size of the enhancement depends on the choice of the matching scale μ W and the μ t scale in m t (μ t ) in the NLO expressions. The implications of these uncertainties for ε /ε are clearly seen in Fig. 1.
What is still missing are NNLO QCD corrections to QCDPs which on the basis of [23,24] are expected to further suppress ε /ε, albeit the effect appears to be smaller than the one of NNLO QCD contributions to EWPs. One could in principle question the inclusion of the latter contributions while leaving out NNLO corrections to QCDPs. Yet these two different NNLO contributions do not have anything to do with each other. In particular while NLO QCD corrections to QCDPs remove already some scale and renormalization scheme dependences present at the LO, in the EWP sector these unphysical scheme dependences are first removed at the NNLO level [38].

Numerical analysis
Our analysis shows that the largest remaining uncertainties in the evaluation of ε /ε are present in the values of Q 6 (m c ) 0 (or B (1/2) 6 ) and eff . In Table 2 we give ε /ε as a function of these two parameters for B ≥ 0.7 also the dependence on eff is significant.
Finally, it is of interest to ask how large a central value of B (1/2) 6 should be in order to reproduce the central experimental value in (1). It turns out to be B (1/2) 6 = 1.40, in total disagreement with (12). The central value in (2) is obtained for B (1/2) 6 = 1.24.

Summary and outlook
Our analysis and in particular the comparison of the results in (13) and (14) as well as the Table 2 demonstrates the importance of NNLO QCD corrections and of isospin-breaking effects. Anticipating that the new RBC-UKQCD analysis will find B (1/2) 6 (m c ) < 1.0 as hinted by DQCD, the values of ε /ε in the SM will be significantly below the data. Our example with B (1/2) 6 (m c ) in the ballpark of 0.80 ± 0.08 illustrates a significant anomaly in ε /ε of about 3.3 σ . If confirmed by new RBC-UKQCD analysis this would turn out to be one of the largest anomalies in flavour physics present in a single observable and comparable to the anomaly in the flavour conserving (g − 2) μ . Moreover, this would be presently the only significant anomaly as far as CP-violation is concerned with the possible exception of the one present in B → π K decays [46], as recently reviewed in [47][48][49][50].
However, even if our expectations for ε /ε in the SM would be confirmed by new RBC-UKQCD results, in order to obtain a better assessment which NP is responsible for this anomaly it is very important to perform a number of the following steps: • Obtain satisfactory precision on Q 6 (m c ) 0 or B (1/2) 6 .
• Reduce the error on eff . In particular isospin-breaking and QED effects should be taken into account in LQCD calculations. • Even if the insight from DQCD allowed us to identify the dynamics (meson evolution) responsible for this anomaly, at least a second lattice QCD collaboration should calculate K → ππ matrix elements and ε /ε. • Further reduce the short-distance uncertainties, in particular in the QCD penguin sector. But the subleading NNLO QCD contributions to the electroweak penguin sector should be also evaluated. • Calculation of BSM K → ππ hadronic matrix elements of four-quark operators by lattice QCD that presently are known only in the DQCD [51].
There have been numerous BSM analyses of ε /ε which we collect in Table 3. Here we just mention that the leptoquark models, with possible exception of the vector U 1 model, are not capable in explaining this anomaly because of the constraints from rare Kaon decays [52]. This shows how crucial correlations of ε /ε with other observables in a given NP scenario are. As indicated in Table 3, they have been analyzed in other NP scenarios. In particular, very recently, a correlation of hinted anomalies in ε /ε and B → π K decays has been pointed out in the context of models with U (2) 3 flavour symmetry in [53].
Also the lessons gained from the SMEFT analysis in [39] should be very helpful in identifying NP behind hinted ε /ε anomaly. Such a general analysis allows to take the constraints from other processes, in particular from electroweak precision tests and collider processes into account. To this end the master formula for ε /ε [54] valid in any extension of the SM should facilitate the search for the dynamics behind the anomaly in question.
Acknowledgements This research was done and financially supported by the DFG cluster of excellence "Origin and Structure of the Universe".

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This manuscript has no associated data, all results are given within the main text.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

A Wilson coefficients
Here we summarize the S = 1 Wilson coefficients at the scale μ = m c = 1.3 GeV in the NDR-MS scheme using the NLO RG evolution from [34]. The numerical input is fixed to values in Table 4   We will use the results in [38] to demonstrate the numerical impact of the dominant NNLO matching corrections that resolve the NLO renormalization scheme ambiguities for the two choices μ t = m W and μ t = m t . As given in Table  2 (Table 3) Table 5, which we adapt in the numerics. The prime in this Table indicates that still small m t -independent O(α W α s sin 2 θ W ) corrections are not included and NNLO corrections to y 3,4,5,6 are neglected as well.

B Error budget
We summarize the error budget leading to the result in (14) in Table 6. The scale uncertainties after the inclusion of NNLO corrections to both QCDP and EWP are not shown as they have negligible impact on the final error. Table 6 Table of the absolute error of ε /ε for benchmark point (12) with input parameters from (4) and (12). The absolute error of ε /ε from these parametric uncertainties becomes 10 4 ×δ(ε /ε) = 2.3 when added in quadrature