Dissecting the $\Delta I= 1/2$ rule at large $N_c$

We study the scaling of kaon decay amplitudes with the number of colours, $N_c$, in a theory with four degenerate flavours, $N_f=4$. In this scenario, two current-current operators, $Q^\pm$, mediate $\Delta S=1$ transitions, such as the two isospin amplitudes of non-leptonic kaon decays for $K\to (\pi\pi)_{I=0,2}$, $A_0$ and $A_2$. In particular, we concentrate on the simpler $K\to\pi$ amplitudes, $A^\pm$, mediated by these two operators. A diagrammatic analysis of the large-$N_c$ scaling of these observables is presented, which demonstrates the anticorrelation of the leading ${\mathcal O}(1/N_c)$ and ${\mathcal O}(N_f/N_c^2)$ corrections in both amplitudes. Using our new $N_f=4$ and previous quenched data, we confirm this expectation and show that these corrections are $naturally$ large and may be at the origin of the $\Delta I=1/2$ rule. The evidence for the latter is indirect, based on the matching of the amplitudes to their prediction in Chiral Perturbation Theory, from which the LO low-energy couplings of the chiral weak Hamiltonian, $g^\pm$, can be determined. A NLO estimate of the $K \to (\pi\pi)_{I=0,2}$ isospin amplitudes can then be derived, which is in good agreement with the experimental value.


I. INTRODUCTION
Significant progress has been achieved recently in the lattice determination of K → (ππ) I=0,2 amplitudes and the CP violating observable / [1][2][3]. In particular, a large enhancement of the I = 0 amplitude over the I = 2 one has been reported, albeit with too large uncertainty to be considered a satisfactory first-principles determination of the ∆I = 1/2 rule.
In Ref. [4] an analysis of the different contributions was made and it was concluded that the main source of the enhancement lies in a strong cancellation of the isospintwo amplitude, as a result of a negative relative sign between the colour-connected and colour-disconnected contractions, with the two contributions adding up in the isospin-zero channel. In Refs. [5][6][7] we proposed to study the N c dependence of the amplitudes, because the two contributions scale differently in large N c and therefore can be rigorously disentangled in this limit. The enhancement, if explained in this fashion, seems to require unnaturally large-N c corrections with the appropriate sign.
Interestingly, the large-N c limit of QCD [8] has also inspired several phenomenological determinations of these observables [9][10][11][12][13][14] (for a recent discussion see [15][16][17]). It is well known, however, that the leadingorder large-N c prediction for the ratio of the amplitudes, lim Nc→∞ A 0 /A 2 = √ 2, i.e., no ∆I = 1/2 rule whatsoever. The subleading N c corrections should therefore be very large, which could be consistent with the previous hypothesis, but casts doubts on the phenomenological approaches that make use of large-N c inspired approximations: if we know that there must be significant large-N c corrections to explain the ∆I = 1/2, why should we trust approximations that neglect subleading N c terms?
The N c dependence can be studied from first-principles in lattice QCD by simply simulating at different number of colours [18][19][20][21][22]. In our previous work [5][6][7] we explored the related weak amplitudes K → π and K →K in the quenched approximation, and found no unnaturally large subleading N c corrections, although we confirmed the exact anticorrelation of these corrections in the two isospin channels. The quenched approximation introduces however an uncontrollable systematic error, which in practice is often found to be relatively small in most quantities. Since we are interested in subleading N c corrections, quenching effects are expected to enter at this order of the N c expansion and therefore need to be included. The main goal of this paper is to extend our previous study beyond the quenched approximation, which will allow us to determine from first-principles the subleading N c corrections to the ∆I = 1/2 rule, in a simplified setting with four degenerate flavours, m u = m d = m s = m c . This paper is organized as follows: in Section II we discuss our strategy for the lattice study of K → π transitions; in Section III we discuss the N c scaling of the amplitudes; Section IV deals with the necessary results in Chiral Perturbation Theory to connect to K → ππ; Section V describes the setup of our lattice computations; in Section VI we discuss our physics results; and we conclude in Section VII.

II. STRATEGY
The Operator Product Expansion allows to represent CP-conserving ∆S = 1 transitions by an effective Hamiltonian of four-fermion operators. At the electroweak scale, µ M W , we can neglect all quark masses, and the weak Hamiltonian takes the simple form: Only two four-quark operators of dimension six can appear with the correct symmetry properties under the flavour symmetry group SU(4) L × SU(4) R , namelȳ where J µ is the left-handed current J ij µ = (ψ i γ µ P − ψ j ); i, j are quark flavour indices; P ± = 1 2 (1 ± γ 5 ); and parentheses around quark bilinears indicate that they are traced over spin and colour. Z ± Q (µ) is the renormalization constant of the bare operator Q ± (x) computed in some regularization scheme as, for example, the lattice. There are other operators that could mix with those above: however, they vanish in the limit of equal up and charm masses, that we refer to as the GIM limit [23]. From the lattice point of view the GIM limit is very advantageous, not only for the simpler operator mixing, but also because no closed quark propagator contributes to the amplitudes. Even though the presence of a heavy charm was argued long ago to be at the origin of the ∆I = 1/2 rule via the mixing with penguin operators [24], the relevance of penguin contributions has been found to be small in non-perturbative studies [1,25]. If we want to test the primary mechanism of the ∆I = 1/2 enhancement proposed in [4], the GIM limit may be good enough.
The operatorsQ σ (µ) are renormalized at a scale µ in some renormalization scheme, being their µ dependence exactly cancelled by that of the Wilson coefficients k σ (µ). It is also possible to define renormalization group invariant (RGI) operators, which are defined by cancelling their µ dependence, as derived from the Callan-Symanzik equations,Q whereḡ(µ) is the running coupling and β(g) = −g 3 n b n g 2n , γ σ (g) = −g 2 n γ σ n g 2n are the β-function and the four-fermion operator anomalous dimension, respectively. The one-and two-loop coefficients of the βfunction, and the one-loop coefficient of the anomalous dimensions, are renormalization scheme-independent. Their values for the theory with N f flavours are [26][27][28][29] and for the operators Q ± [30,31] The normalization ofĉ σ (µ) coincides with the most popular one for N c = 3, whilst using the 't Hooft coupling λ = N cḡ 2 (µ) in the first factor instead of the usual coupling, so that the large-N c limit is well-defined.
Defining similarly an RGI Wilson coefficient we can rewrite the Hamiltonian in terms of RGI quantities, which no longer depend on the scale, so that where µ is a convenient renormalization scale for the non-perturbative computation of matrix elements of Q ± , which will be later set to the inverse lattice scale a −1 .
The factor U σ (µ, M W ) =ĉ σ (µ)/ĉ σ (M W ), therefore, measures the running of the renormalized operator between the scales µ and M W . Ideally one would like to evaluate this factor non-perturbatively, as has been done for N c = 3 [32,33], but such a challenging endeavour is beyond the scope of this paper. We will instead use the perturbative results at two loops in the RI scheme [34,35] to evaluate the Wilson coefficients k σ (M W ), the running factors U σ (µ, M W ), andĉ(µ). This implies relying on perturbation theory at scales above µ ≥ a −1 ∼ 2.6 GeV. Similarly we will also use lattice perturbation theory to estimate the renormalization factors Z ± Q , that are known to one loop [36,37].
We are interested in considering K → π amplitudes in the two isospin channels, that we can extract from ratios of three-point correlators where and the two-point correlators with A ij 0 (x) ≡ψ i (x)γ 0 γ 5 ψ j (x). From these correlators we define the bare lattice ratios: x,y C ± 3 (y, z, x) x,y C du 2 (y, z)C us 2 (x, z) .
The renormalization factors for these ratios, Z ± , are obtained from the ratio of the renormalization factors of    [36,37], whereas U σ and k σ are computed using the two-loop MS coupling. The star labels the simulation points with finer lattice spacing, a ∼ 0.065 fm.
the four fermion operators, and the current normalization factors that appear in the denominator. From the renormalized ratios we can obtain the RGI normalized ratioŝ and the K → π amplitudes, written either in terms of the RGI or the renormalized ratios, as All the required factors to reconstruct the physical amplitudes are summarized in Table I for N f = 4 (this work), and in Table II for the quenched case [5,6].
A simple diagrammatic analysis of the three and two point correlators of Eqs. (10,12) shows a clear pattern of the large-N c scaling, and demonstrates the expected anticorrelation of the leading large-N c corrections of the A ± amplitudes.
After integration over fermion fields, the correlators are obtained from the gauge averages of the colourdisconnected and colour-connected contractions of Fig. 1, corresponding to the operator insertion O suud and O sduu ,    [5,6]. Z σ (a −1 ) have been computed at one loop in tadpole-improved perturbation theory using the results in [36,37], whereas U σ and k σ are computed using the two-loop MS coupling. Note that the values of Z σ (a −1 ) differ from those in Refs. [5,6], where bare lattice perturbation theory was used. Furthermore, the values of k σ and U σ also supersede the ones in Refs. [5,6]. In the evaluation ofĉ σ (a −1 ) we have used Λ MS as described in Ref. [5]. respectively. In Figs. 2 and 3 we show the scaling with N c of the lowest-order diagrams contributing to these correlators. The leading N c dependence of both the renormalized and bare correlators are therefore of the form: where all the coefficients a − f in these expressions (each of them related to one or more diagrams in Figs. 2 and 3) are independent of N c and N f . These relations imply (e) They are therefore fully anticorrelated in the ± correlators. Importantly, the anticorrelated terms include the leading fermion loop corrections, O(N f ). These relations also imply the following scaling of the renormalization factors: and a similar one for the Wilson coefficients, k σ . This dependence can be explicitly checked in the perturbative coefficients known up to two loops in the MS scheme [34,35]. These results imply the following scaling of the amplitudes: where the coefficientsã −d are combinations of the coefficients a − f in Eq. (17), and are also independent of N c and N f , and a natural expectation is that they are O(1). Not only the leading corrections N −1 c are, therefore, fully anticorrelated in the ratios, but also the leading effects of dynamical quarks, O(N f ). Note that this analysis does not predict the sign of the different terms, i.e., the sign of theã −d coefficients, only the (anti)-correlation between the two isospin channels. This way, a negative sign ofã andb results into an enhancement of the ratio 't Hooft vs. Veneziano scaling As we will see the number of active flavours, N f , plays a relevant role in the 1/N c expansion of the K → π amplitudes. The scaling in N f is in fact the difference between the 't Hooft and Veneziano limits of QCD. While the former keeps N f constant when taking N c → ∞, the latter keeps the ratio N f /N c constant. From Eq. (20), it is then clear thatã andb have the same scaling in the Veneziano limit (the same holds forc andd). In our simulations, we will be studying the 't Hooft limit, since we keep N f fixed, but the quantity N f /N c is large (ranging from 4/3 to 2/3, depending on N c ), so its contribution may be very significant even for naturally largeã −d coefficients.
The chiral dependence of the ratios in Eq.(13) can be studied within the framework of Chiral Perturbation Theory (ChPT) with N f = 4 active flavours. An extensive discussion of this framework can be found in Refs. [23,39]. Here we just summarize the required formulae, and refer to those references for details.
The weak Hamiltonian in Eq. (1) can be translated to an effective weak Hamiltonian in terms of meson fields preserving the flavour symmetries. Since the operators Q + andQ − transform under representations of SU (4) L of dimension 84 and 20, their ChPT counterparts must be constructed accordingly. At leading order, there are only two terms, with couplings g ± , that need to be de-termined non-perturbatively: with where U is the chiral meson field, i, j, k, l are flavour indices, and c σ ijkl are Clebsch-Gordan coefficients (see Appendix A in Ref. [23]).
By means of the chiral weak Hamiltonian in Eq. (21) and the standard NLO ChPT Lagrangian, the chiral predictions for the amplitudes in Eq. (16) are found to be: where L r ± are the NLO counterterms. The NLO corrections in Eq. (23) are fully anticorrelated. Extrapolating the ratios in Eq. (13) to zero pion mass, one can determine the leading low-energy couplings (LECs) of the chiral weak Hamiltonian: The extracted values of g ± can then be used to make predictions of other observables, such as the K → ππ decay amplitudes. We now turn to the analysis of the combined chiral and N c dependence. First, we note that Eq. (20) should hold at any pion mass, and therefore we expect: Furthermore, by comparing the chiral dependence in Eq. (23) with the N c scaling in Eq. (20) we can see that both L r + and L r − must be O(N 0 c ), and identical at this order. The next term in the 1/N c expansion for L r ± could in principle differ: Hence, the combination of Eq. (23) with Eqs. (25,26) can be used to do global fits including different meson masses and values of N c . It will be convenient to also study the chiral and N c dependence of the product of A + A − . The reason is that the leading chiral and N c corrections cancel out, which leads to a more robust chiral extrapolation. The chiral corrections for this quantity are with where α and β depend on the coefficients a χ − d χ .
B. Relation to K → ππ amplitudes Once the effective couplings g ± have been extracted from the chiral extrapolations of the ratios A ± , they can be used to compute the K → ππ weak decay amplitudes. The two pions in the final state can be in a state with total isospin I = 0 or 2: The ratio of the two amplitudes can be calculated at leading order in ChPT using the Hamiltonian in Eq. (21) [23,40]: The measured hierarchy of ∼ 22 between A 0 and A 2 must then be translated into a large ratio of the couplings g ± . Note that for g + = g − = 1, the expected large-N c result is recovered, A 0 /A 2 = √ 2. Large 1/N c corrections in the g − /g + ratio could therefore be the origin of the ∆I = 1/2 rule.
We have also derived the ChPT NLO result for the non-degenerate case in which we send the pion mass to zero, while keeping the kaon mass at its physical value 1 . As we are forced to work in the exact GIM limit, we must also send the charm quark mass to zero with the up quark mass. The calculation for m s > m u = m d = m c = 0 yields: where Λ eff is an unknown scale that contains information of the NLO LECs of the effective Chiral Lagrangian and the effective weak Hamiltonian. We note that the NLO effect tends to enhance (reduce) the ratio for Λ eff > M K (Λ eff < M K ).

A. Simulation and matching of sea and valence sectors
Our lattice setup is the same as the one presented in Ref. [22], and we refer to it for details on the simulations and scale setting. We use ensembles with N f = 4 dynamical fermions for an SU (N c ) gauge theory, with N c = 3 − 6. They have been generated using the HiRep code [42,43]. We have chosen the Iwasaki gauge action (following previous experience with 2+1+1 simulations [44]) and clover Wilson fermions for the sea quarks, with the plaquette-boosted one-loop value of c sw . The simulation parameters are shown in Table III. The lattice spacing is found to be a ∼ 0.075 fm for all values of N c (see also Ref. [22]). In addition, we have produced two ensembles with a finer lattice spacing, a ∼ 0.065 fm, to estimate discretization effects.
In order to achieve automatic O(a) improvement 2 [47] and avoid the mixing of different-chirality operators for weak decays, we employ maximally twisted valence quarks [48], i.e., the mixed-action setup [49] previously used in Refs. [45,46]. Working in twisted quark field variables, maximal twist is ensured by tuning the untwisted bare valence mass m v to the critical value for which the 2 As discussed in [45,46], there are residual O(a) cutoff effects from virtual sea quarks, which are proportional to am s and carry coefficients that are O(α 2 s ) in perturbation theory. These effects are expected to be numerically very small and thus irrelevant for the discussion below. It is also worth stressing that using the one-loop value of csw will also lead to residual effects of O(a α 2 s ).
valence PCAC mass is zero: The bare twisted mass parameter µ 0 is tuned such that the pion mass in the sea and valence sectors coincide, M v π = M s π . Since twisted mass already provides O(a) improvement, the clover improvement parameter c sw can be chosen to be an arbitrary value in the valence sector. We choose c sw = 0 in the valence sector 3 for this work, our main motivation being that this minimizes the isospin breaking effects coming from the twisted-mass action. In addition, this will allow for a partial crosscheck of the systematics due to the use of perturbative renormalization constants, by comparing the latter to the nonperturbative determination in Ref. [50] for N c = 3 (see below). Finally, we also observe that c sw = 0 leads to smaller statistical errors.
In Table IV we present our measurements for the ensembles used in this work. We have achieved good tuning to maximal twist, with the PCAC mass being zero within 1 or 2σ. In addition, the valence and sea pion masses are matched also within 1 or 2σ. The bare results for the ratios are also presented in the same table, together with the chiral parameter ξ = M 2 π /(4πF π ) 2 , that will be used for the chiral extrapolations We conclude the discussion of the simulation setup by mentioning that we will compare the new results with dynamical fermions to the ones in Refs. [5,6]. Those results used quenched simulations, with plaquette gauge action and twisted mass fermions. The lattice spacing was a ∼ 0.093 fm and the the pion mass was fixed at around M π = 550 − 590 MeV for N c = 3 − 8 and 17.
In this work, we perform a reanalysis of these quenched data.

B. Comments on systematics
We conclude this section by discussing the systematic errors that can affect our results.
We start with finite-volume effects. Our ensembles have M π L > 3.8 in all cases so we expect finite-volume effects to be small, and suppressed as 1/N c . Still, we find that for the observable ξ they can be of O(1%) and thus we correct for them, as explained in Ref. [22], following Refs. [51,52]).
Since B K andR + differ by a volume-independent proportionality factor, we can use the results in Ref. [53], where the finite-volume effects of B K have been calculated. In addition, it is known that the finite-volume and  The value of the lattice spacing is a 0.075 fm for the "A" ensembles (see Ref. [22]), whereas it is a 0.065 fm for "B" ensembles. We provide the pion mass in the valence sector, aM v π , and the PCAC mass, am v pcac . We also include the results for the ratios in Eq. (13), and in the last column, the chiral parameter ξ ≡ M 2 π /(4πFπ) 2 . Moreover, ξL labels ξ corrected by finite-volume effects as explained in the main text. chiral corrections ofR + andR − are fully anticorrelated [39]. Thus, we find: The correction for these quantities is numerically negligible for our ensembles. While additional finite-volume effects could be present (see Ref. [52]) we observe that a factor of two increase or decrease of these finite-volume corrections alters our results well within the statistical precision.
Concerning discretization effects, we have included the results from two ensembles with a finer lattice spacing at N c = 3. Assuming O(a) improvement, we expect that the finer lattice spacing should reduce by ∼ 30% the O(a 2 ) discretization effects. We observe no significant difference for these data points in Fig. 6, so we see no sign of sizeable discretization errors within our statistical uncertainty. We stress however that a more extensive study is needed for a robust estimate of the discretization error.
The largest systematic error that we have found is related to the renormalization constants, which we have estimated by one-loop perturbation theory. We have first compared the non-perturbative renormalization constants of Ref. [50] to the one-loop perturbation theory results in their setup (they used c sw = 0). The difference is roughly ∼ 5% for N c = 3. On the other hand, we have computed the ratios using c sw = 1.69 in the valence sector for the 3A10 ensemble. Using the perturbative renormalization constants for this new value of c sw we get a result that differs from our c sw = 0 result by roughly 20% in the ratio. Since it is unlikely that this effect can be accounted for by discretization effects, given the tests in a finer lattice mentioned above, we conclude that there must be significant non-perturbative effects on renormalization constants for the larger c sw (the perturbative one-loop corrections are also significantly larger for the larger value of c sw ). This is a large error, and probably a conservative estimate, but it is comparable to the statistical error we achieve, as it will be seen later.

A.
Nc scaling of K → π amplitudes The physical amplitudes A ± can be obtained, as explained in Eq. (16), from the bare ratios in Table IV,  and the renormalization coefficients in Tables I and II. As explained above, a rigorous way to isolate the (anti-)correlated contributions to the ratios consists on taking the half-sum and half-difference of the ratios. By doing so, the two contributions can be fitted independently since: In the following, we compare the results of the fits to Eq. (35) in three different scenarios: 1. Quenched results (N f = 0) at a heavy pion mass ∼ 570 MeV.
The results for the coefficientsã −d for the three scenarios are presented in Table V and Fig. 4. The coefficients are all of O(1) and therefore of natural size. Importantly the sign of theã andb coefficients is the same and negative. This implies both terms contribute to reduce the A + amplitude and enlarge, in a correlated way, the amplitude A − . The fact thatb,d ∼ O(1) implies a very large unquenching effect in the large-N c scaling, and the ratio A − /A + , which is however compatible with the expansion in Eq. (35). Specifically, it is due tob andd being absent for N f = 0. The other two coefficients,ã andc, are comparable in size in the quenched and dynamical theories. We note however that uncertainties only include statistical errors, and relative discretization errors and the systematics of the perturbative renormalization constants may be significant. Finally, we observe that the mass dependence for the N f = 4 results seems to affect mostly the coefficientã, which is consistent with the chiral dependence in Eq. (23), and goes also in the direction of enhancing the ratio A − /A + towards the chiral limit.

B. Kaon B-parameter (BK )
The kaon B-parameter, B K , is defined from the matrix element of the ∆S = 2 operator that mediates neutral kaon oscillations at physical kinematics: It is customary to quote the renormalization group independent (RGI) version, labelled asB K . Its value at the physical point has been computed accurately in N f = 2, 2 + 1, and 2 + 1 + 1 simulations [50,[54][55][56][57][58] (see Ref. [59] for a review). In our setup,B K coincides with the renormalized ratioR + up to a normalization. Specifically, we havê whereĉ + can be read off Table I. There are two essential differences in our setup: all meson masses are degenerate, in particular M K = M π , and we have an active light charm quark. Both can significantly affect the value of B K . We show our results in Fig. 5. We observe a very significant N c dependence ofB K for N f = 4, and a much milder one for N f = 0. For N c = 3, the quenched result agrees with the standard value ofB K , while the N f = 4 result is about 25% smaller. We have included as bands the Buras-Bardeen-Gerard (BBG) phenomenological model prediction from Ref. [15], using inputs on meson masses from our own simulations in both casesquenched and dynamical. We find that our results are reasonably compatible with the BBG prediction, in particular regarding the suppression ofB K in the presence of a light charm.
To conclude this subsection, we can use the scaling in N c to infer a value ofB K with three active flavours and quasi-physical kinematics. For this, we use the coefficientsã −d in Table V for the case of N f = 4 and M π = 560 MeV, and so predict the value of A + with N c = 3 and N f = 3 at the same value of the pion mass, degenerate with the kaon. We can the get the RGI valuê B K as in Eq. (37), extractingR + and using theĉ + (a −1 ) for three-flavour QCD 4 . We find including statistical error, and a ∼ 10% error due to the systematics of the renormalization constants. We also quote a "fit" error that we estimate by using the N c scaling derived from a direct fit of the half-sum and difference ofR ± instead of A ± . We have not found results in the literature for the degenerate case that we can compare to. On the other hand, ChPT relates the value ofB K in the degenerate case, to the quasi-physical (QP) situation with M π = 0 and M K at its physical value: where Λ B K eff labels an unknown scale that parametrizes the effect of the unknown LECs. For Λ B K eff > M K ,B QP K is larger thanB K and could be compatible with the existing results at the physical point from N f = 2 + 1, N c = 3 simulations [50,[54][55][56][57][58].  [5,6]), and N f = 4 (this work). Error bars are only statistical errors. We also include the predictions from Ref. [15], where the band indicates the values obtained when varying the (arbitrary) matching scale M from 600 to 1000 MeV.

C.
Extraction of the effective couplings g ± The main goal of this work is to compute the ratio g − /g + by extrapolating A ± to the chiral limit. For the required chiral extrapolation, we follow the same strategy as in Ref. [40]. We extract g + from a chiral fit to A + , and the product g + g − from that of the product A + A − . The ratio can then be evaluated as This approach results in a milder chiral extrapolation, that will hopefully introduce a smaller systematic error. We have performed two kinds of fits. In Fit 1, we use all data points with N c = 3 − 6 in a simultaneous chiral and N c fit using Eqs. (23) and (27), incorporating the 1/N c expansion of the couplings as in Eqs. (25,26,29). In Fit 2, we fit using only the data with N c = 3, and extract the effective couplings for this theory. This way, for N c = 3 we find: Fit 1: g + = 0.187 (21), g + g − = 0.91(4), Fit 2: g + = 0.190 (27), g + g − = 0.80 (6).
The complete results of these fits are shown in Tables VI, and VII, and also in Fig. 6. From these results, we obtain for the ratio of couplings at N c = 3: where errors are only statistical, but correlations are taken into account.
-3(4) 7(7) 2.4(8) -11(4) 12.0/11    Table IV. Empty squares for Nc = 3 indicate a finer lattice spacing. Solid lines indicate a simultaneous chiral and Nc fit as in Eq. (23). Dashed lines represent the chiral extrapolation of the data points for Nc = 3 following Eqs. (23) and (27). Errors are only statistical. K → ππ decay for N c = 3. In Fig. 7, we show this prediction as a function of an unknown effective scale Λ eff . This prediction, valid for M π = M D = 0 and physical M K , shows small NLO effects in a wide range of values of the effective scale.
We are now in the position to quote a final result for the ratio of isospin amplitudes: In the previous equation, the various error sources originate as follows : (i) statistical error, (ii) systematic error from the difference between fit 1 and 2 in Eq. (42), (iii) a 20% error from the renormalization constants -see Section V B -, and (iv) a 10% error from the NLO effects -see Fig. 7. Combining all error sources in quadrature results in a ∼ 30% uncertainty on the total result, which is dominated by systematics. We also stress that this is a result in the theory with a light charm quark. Interestingly, this indirect computation yields a value compatible with the experimental result for the ∆I = 1/2 enhancement. We use the input of Fit 2 in Eq. (42). This prediction is valid for Mπ = MD = 0, and MK at its physical value. The shaded area represents the statistical error associated to the ratio of couplings -see Eq. (42). As a guideline, we also show the experimental value for the ratio of amplitudes (in blue).

VII. CONCLUSIONS
We have presented the first non-perturbative study of the scaling of ∆S = 1 weak amplitudes with the number of colours, N c = 3 − 6, in a theory with four degenerate light flavours N f = 4. These results have been obtained from dynamical simulations with clover Wilson fermions, at a 0.075 fm and a 0.065 fm and pion masses in the range 360 − 570 MeV. We have analysed the K → π amplitudes A ± , mediated by the two current-current operators Q ± of the ∆S = 1 weak Hamiltonian in Eq. (1).
The diagrammatic analysis of the large-N c scaling of these observables presented in Sect. III allows to classify the subleading N c corrections, and demonstrates the anticorrelation of the leading O(1/N c ) and O(N f /N 2 c ) con-tributions in the A ± amplitudes. Our numerical results confirm this expectation and show that these corrections are naturally large in the Veneziano scaling limit, i.e., the coefficients of both corrections are O(1). They can nevertheless explain the large enhancement of the ratio A − /A + for N c = 3 with respect to the N c → ∞ limit. This involves an unprecedentedly large unquenching effect in this ratio, that is nevertheless compatible with natural size O(N f /N 2 c ) corrections. The amplitudes A ± in the chiral limit can be matched to their ChPT counterparts, which depend on the leading low-energy couplings, g ± , of the chiral effective weak Hamiltonian. From a chiral extrapolation of the combinations A + and A + A − , we have then extracted the couplings g ± , which are finally used to predict in ChPT the K → (ππ) I=0,2 amplitudes. In particular, we have obtained an indirect prediction of the ratio of isospin amplitudes, A 0 /A 2 , by this procedure which seems to largely account for the elusive "∆I = 1/2 rule". Our estimate for this ratio in the theory with a light charm is which suggests that the enhancement may indeed be largely dominated by intrinsic QCD effects.