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

Abstract

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\rightarrow (\pi \pi )_{I=0,2}$, $A_0$ and $A_2$. In particular, we concentrate on the simpler $K\rightarrow \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 \rightarrow (\pi \pi )_{I=0,2}$ isospin amplitudes can then be derived, which is in good agreement with the experimental value.

1 Introduction

Significant progress has been achieved recently in the lattice determination of $K \rightarrow (\pi \pi )_{I=0,2}$ amplitudes and the CP violating observable $\epsilon '/ \epsilon$ [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 $\Delta I=1/2$ rule.¹

In Ref. [5] an analysis of the different contributions was made and it was suggested that the main source of the enhancement lies in a strong cancellation of the isospin-two 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. [6,7,8] 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 [9, 10] has also inspired several phenomenological determinations of these and related observables [11,12,13,14,15,16,17,18,19] (for a recent discussion see [20,21,22]). It is well known, however, that the leading-order large-$N_c$ prediction for the ratio of the amplitudes, $\lim _{N_c \rightarrow \infty } A_0/A_2 = \sqrt{2}$, i.e., no $\Delta 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 $\Delta 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 [23,24,25,26,27]. In our previous work [6,7,8] we explored the related weak amplitudes $K \rightarrow \pi$ and $K \rightarrow \bar{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 $\Delta 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 Sect. 2 we discuss our strategy for the lattice study of $K\rightarrow \pi$ transitions; in Sect. 3 we discuss the $N_c$ scaling of the amplitudes; Sect. 4 deals with the necessary results in Chiral Perturbation Theory to connect to $K\rightarrow \pi \pi$; Sect. 5 describes the setup of our lattice computations; in Sect. 6 we discuss our physics results; and we conclude in Sect. 7.

2 Strategy

The Operator Product Expansion allows to represent CP-conserving $\Delta S=1$ transitions by an effective Hamiltonian of four-fermion operators. At the electroweak scale, $\mu \simeq M_W$, we can neglect all quark masses, and the weak Hamiltonian takes the simple form:

$$\begin{aligned} H_{\mathrm{w}}^{\Delta S=1} = \int d^4x~\frac{g_{\mathrm{w}}^2}{4M_W^2}V_{us}^*V_{ud}\sum _{\sigma =\pm } k^\sigma (\mu ) \, \bar{Q}^\sigma (x,\mu ), \end{aligned}$$

(1)

where $g_{\mathrm{w}}^2=4\sqrt{2}G_{\mathrm{F}} M_W^2$. Only two four-quark operators of dimension six can appear with the correct symmetry properties under the flavour symmetry group $\mathrm{SU}(4)_{\mathrm{L}} \times \mathrm{SU}(4)_{\mathrm{R}}$, namely

$$\begin{aligned} \begin{aligned} {{\bar{Q}}}^\pm (x,\mu )&= Z_Q^\pm (\mu ) \, \big ( J_\mu ^{su}(x)J_\mu ^{ud}(x) \pm J_\mu ^{sd}(x)J_\mu ^{uu}(x) \\&\quad -[u\leftrightarrow c]\big ), \end{aligned} \end{aligned}$$

(2)

where $J_\mu$ is the left-handed current $J_\mu ^{ij} = ({{\bar{\psi }}}^i\gamma _\mu P_-\psi ^j)$; i, j are quark flavour indices; $P_\pm ={1\over 2} (\mathbf {1}\pm \gamma _5)$; and parentheses around quark bilinears indicate that they are traced² over spin and colour. $Z_Q^\pm (\mu )$ is the renormalization constant of the bare operator $Q^\pm (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 [28]. 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 $\Delta I=1/2$ rule via the mixing with penguin operators [29], the relevance of penguin contributions has been found to be small in non-perturbative studies [1, 30].³ If we want to test the primary mechanism of the $\Delta I=1/2$ enhancement proposed in [5], the GIM limit may be good enough.

The operators ${{\bar{Q}}}^\sigma (\mu )$ are renormalized at a scale $\mu$ in some renormalization scheme, being their $\mu$ dependence exactly cancelled by that of the Wilson coefficients $k^\sigma (\mu )$. It is also possible to define renormalization group invariant (RGI) operators, which are defined by cancelling their $\mu$ dependence, as derived from the Callan-Symanzik equations,

$$\begin{aligned} \hat{Q}^\sigma \equiv \hat{c}^\sigma (\mu ) {{\bar{Q}}}^\sigma (\mu ), \end{aligned}$$

(3)

with

$$\begin{aligned} {{\hat{c}}}^\sigma (\mu )\equiv&\left( {N_c \over 3} {\bar{g}^2(\mu )\over 4 \pi } \right) ^{{\gamma ^\sigma _0\over 2 b_0}} \nonumber \\&\times \exp \left\{ -\int _0^{\bar{g}(\mu )}\mathrm{d}g \left[ {\gamma ^\sigma (g)\over \beta (g)} - {\gamma ^\sigma _0\over b_0 \, g}\right] \right\} , \end{aligned}$$

(4)

where $\bar{g}(\mu )$ is the running coupling and $\beta (g)=-g^3 \sum _n b_n g^{2 n}$, $\gamma ^\sigma (g) = -g^2 \sum _n \gamma ^\sigma _n g^{2 n}$ are the $\beta$-function and the four-fermion operator anomalous dimension, respectively. The one- and two-loop coefficients of the $\beta$-function, and the one-loop coefficient of the anomalous dimensions, are renormalization scheme-independent. Their values for the theory with $N_f$ flavours are [31,32,33,34,35,36]

$$\begin{aligned} b_0= & {} \frac{1}{(4\pi )^2}\left[ \frac{11}{3}N_c-\frac{2}{3}N_f\right] \,, \end{aligned}$$

(5)

$$\begin{aligned} b_1= & {} \frac{1}{(4\pi )^4}\left[ \frac{34}{3}N_c^2-\left( \frac{13}{3}N_c-\frac{1}{N_c}\right) N_f\right] \,, \end{aligned}$$

(6)

and for the operators $Q^\pm$ [37, 38]

$$\begin{aligned} \gamma _0^\pm = \frac{1}{(4\pi )^2}\left[ \pm 6 - \frac{6}{N_c}\right] \,. \end{aligned}$$

(7)

The normalization of $\hat{c}^\sigma (\mu )$ coincides with the most popular one for $N_c=3$, whilst using the ’t Hooft coupling $\lambda = N_c \bar{g}^2 (\mu )$ 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

$$\begin{aligned} \hat{k}^\sigma \equiv {k^\sigma (\mu )\over \hat{c}^\sigma (\mu )}, \end{aligned}$$

(8)

we can rewrite the Hamiltonian in terms of RGI quantities, which no longer depend on the scale, so that

$$\begin{aligned} \begin{aligned} \hat{k}^\sigma \, \hat{Q}^\sigma&= \left[ { k^\sigma (M_W) \over \hat{c}^\sigma (M_W) } \right] \, \left[ {{\hat{c}}}^\sigma (\mu )\, \bar{Q}^\sigma (\mu ) \right] \\&= k^\sigma (M_W) \, U^\sigma (\mu ,M_W) \, \bar{Q}^\sigma (\mu ), \end{aligned} \end{aligned}$$

(9)

where $\mu$ is a convenient renormalization scale for the non-perturbative computation of matrix elements of $Q^\pm$, which will be later set to the inverse lattice scale $a^{-1}$. The factor $U^\sigma (\mu , M_W) = {{\hat{c}}}^\sigma (\mu )/ {{\hat{c}}}^\sigma (M_W)$, therefore, measures the running of the renormalized operator between the scales $\mu$ and $M_W$. Ideally one would like to evaluate this factor non-perturbatively, as has been done for $N_c=3$ [39, 40], 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 [41, 42] to evaluate the Wilson coefficients $k^\sigma (M_W)$, the running factors $U^\sigma (\mu ,M_W)$, and ${\hat{c}}(\mu )$. This implies relying on perturbation theory at scales above $\mu \ge a^{-1} \sim 2.6~\mathrm{GeV}$. Similarly we will also use lattice perturbation theory to estimate the renormalization factors $Z_Q^\pm$, that are known to one loop⁴ [43, 44].

We are interested in considering $K\rightarrow \pi$ amplitudes in the two isospin channels, that we can extract from ratios of three-point correlators

$$\begin{aligned} \begin{aligned}&{{\mathcal {C}}}_3^\pm (y,\,z,x) \\&\quad \equiv \langle P^{du}(y) [O^{suud}(z)\pm O^{sduu}(z)] P^{us}(x) \rangle , \end{aligned} \end{aligned}$$

(10)

where

$$\begin{aligned} P^{ij}(x) \equiv {\bar{\psi }}^i(x)\gamma _5 \psi ^j(x),\;\;O^{ijkl} \equiv {\bar{\psi }}^i \gamma _\mu \psi ^j {\bar{\psi }}^k \gamma _\mu \psi ^l, \end{aligned}$$

(11)

and the two-point correlators

$$\begin{aligned} {{\mathcal {C}}}^{ij}_{2}(y,z)\equiv & {} \langle P^{ij}(y) A_0^{ji}(z) \rangle , \end{aligned}$$

(12)

with $A_0^{ij}(x) \equiv {\bar{\psi }}^i(x)\gamma _0 \gamma _5 \psi ^j(x)$.

From these correlators we define the bare lattice ratios:

$$\begin{aligned} R^\pm = \lim _{ \begin{array}{c} z_0-x_0\rightarrow \infty \\ y_0-z_0\rightarrow \infty \end{array}} \frac{\sum _{{{\mathbf {x}}},{{\mathbf {y}}}}{\mathcal C}_3^\pm (y,z,x)}{\sum _{{{\mathbf {x}},{\mathbf {y}}}} {\mathcal C}^{du}_2(y,z) {{\mathcal {C}}}^{us}_2(x,z)}, \end{aligned}$$

(13)

which are proportional to the $K \rightarrow \pi$ matrix elements with a convenient normalization. The renormalization factors for these ratios, $Z^\pm$, are obtained from the ratio of the renormalization factors of the four fermion operators, and the current normalization factors that appear in the denominator.

From the renormalized ratios

$$\begin{aligned} \bar{R}^\sigma = Z^\sigma R^\sigma , \end{aligned}$$

(14)

we can obtain the RGI normalized ratios

$$\begin{aligned} \hat{R}^\sigma = \hat{c}\left( a^{-1}\right) Z^\sigma R^\sigma , \end{aligned}$$

(15)

and the normalized⁵ $K\rightarrow \pi$ amplitudes, written either in terms of the RGI or the renormalized ratios, as

$$\begin{aligned} A^\sigma = \hat{k}^\sigma \hat{R}^\sigma = k^\sigma (M_W) U^\sigma (a^{-1},M_W) \bar{R}^\sigma . \end{aligned}$$

(16)

All the required factors to reconstruct the physical amplitudes are summarized in Table 1 for $N_f=4$ (this work), and in Table 2 for the quenched case [6, 7].

Table 1 Perturbative renormalization constants and RG running factors for the ensembles with $N_f=4$. $Z^\sigma \left( a^{-1}\right)$ have been computed at one loop in tadpole-improved perturbation theory using the results in [43, 44], whereas $U^\sigma$ and $k^\sigma$ are computed using the two-loop $\overline{\mathrm{MS}}$ coupling. The star labels the simulation points with finer lattice spacing, $a \sim 0.065$ fm. In the evaluation of $\hat{c}^\sigma \left( a^{-1}\right)$ we have used $\Lambda _{\overline{\text {MS}}}(N_f=4)= 298$ MeV from Ref. [45]

Table 2 Perturbative renormalization constants and RG running factors for the runs with $N_f=0$ of Refs. [6, 7]. $Z^\sigma (a^{-1})$ have been computed at one loop in tadpole-improved perturbation theory using the results in [43, 44], whereas $U^\sigma$ and $k^\sigma$ are computed using the two-loop $\overline{\mathrm{MS}}$ coupling. Note that the values of $Z^\sigma (a^{-1})$ differ from those in Refs. [6, 7], where bare lattice perturbation theory was used. Furthermore, the values of $k^\sigma$ and $U^\sigma$ also supersede the ones in Refs. [6, 7]. In the evaluation of $\hat{c}^\sigma (a^{-1})$ we have used $\Lambda _{\overline{\text {MS}}}$ as described in Ref. [6]

3 large-$N_c$ scaling of $K\rightarrow \pi$ amplitudes

3.1 Diagrammatic expansion of $A^\pm$

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^\pm$ amplitudes.

After integration over fermion fields, the correlators are obtained from the gauge averages of the colour-disconnected and colour-connected contractions of Fig. 1, corresponding to the operator insertion $O^{suud}$ and $O^{sduu}$, 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:

$$\begin{aligned} \langle P^{ij} J_\mu ^{ji} \rangle= & {} N_c \left( a + b {N_f\over N_c}\right) +\ldots , \nonumber \\ \langle P^{du} O^{suud} P^{us} \rangle= & {} \langle P^{du} J_\mu ^{ud} \rangle \langle P^{su} J_\mu ^{us} \rangle + c+ d {N_f\over N_c} + \ldots ,\nonumber \\ \langle P^{du} O^{sduu} P^{us} \rangle= & {} N_c\left( e + f {N_f\over N_c}\right) + \ldots , \end{aligned}$$

(17)

Fig. 1

Left diagram: $O^{suud}(x)$ insertion or colour-disconnected contribution to $\mathcal{C}^\pm _3$ in Eq. (10). Right diagram: $O^{sduu}(x)$ insertion or colour-connected contribution to $\mathcal{C}^\pm _3$ in Eq. (10)

Fig. 2

$N_c, N_f$ scaling of various contributions to the colour-disconnected contraction, corresponding to the $O^{suud}(x)$ insertion

Fig. 3

$N_c, N_f$ scaling of various contributions to the colour-connected contraction, corresponding to the $O^{sduu}(x)$ insertion

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 that the leading $N_c$ corrections in the ± correlation functions of Eq. (10) are of ${{\mathcal {O}}}(N_c^2, N_f N_c)$, but factorizable. On the other hand, the leading non-factorizable corrections are of ${{\mathcal {O}}}(N_c)$ and ${{\mathcal {O}}}(N_f)$, and cancel in the sum of the ± correlators:

$$\begin{aligned} {{\mathcal {C}}}_3^+ + {{\mathcal {C}}}_3^-= & {} \mathrm{disconnected}+ {\mathcal O}(N_c^0) + {{\mathcal {O}}}\left( {N_f\over N_c}\right) + \cdots ,\nonumber \\ {{\mathcal {C}}}_3^+ - {{\mathcal {C}}}_3^-= & {} {{\mathcal {O}}}(N_c) + {{\mathcal {O}}}(N_f) + \cdots \end{aligned}$$

(18)

They are therefore fully anticorrelated in the ± correlators. Importantly, the anticorrelated terms include the leading fermion loop corrections, ${{\mathcal {O}}}(N_f)$. These relations also imply the following scaling of the renormalization factors:

$$\begin{aligned} {Z_Q^++Z_Q^-\over 2}= & {} 1+ {{\mathcal {O}}}\left( {1\over N_c^2}\right) + {{\mathcal {O}}}\left( {N_f\over N_c^3}\right) +\cdots \nonumber \\ {Z_Q^+-Z_Q^-\over 2}= & {} {{\mathcal {O}}}\left( {1\over N_c}\right) + {{\mathcal {O}}}\left( {N_f \over N_c^2} \right) +\cdots , \end{aligned}$$

(19)

and a similar one for the Wilson coefficients, $k^\sigma$. This dependence can be explicitly checked in the perturbative coefficients known up to two loops in the $\overline{\mathrm{MS}}$ scheme [41, 42].

These results imply the following scaling of the amplitudes:

$$\begin{aligned} A^\pm = 1 \pm {\tilde{a}} {1 \over N_c}\pm {\tilde{b}} {N_f \over N^2_c}+{\tilde{c}} {1 \over N^2_c}+ {\tilde{d}} {N_f \over N^3_c}+\cdots , \end{aligned}$$

(20)

where the coefficients ${\tilde{a}}- {\tilde{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 ${{\mathcal {O}}}(1)$.

Not only the leading corrections $N_c^{-1}$ are, therefore, fully anticorrelated in the ratios, but also the leading effects of dynamical quarks, ${{\mathcal {O}}}(N_f)$. Note that this analysis does not predict the sign of the different terms, i.e., the sign of the ${\tilde{a}}-{\tilde{d}}$ coefficients, only the (anti)-correlation between the two isospin channels. This way, a negative sign of ${\tilde{a}}$ and ${\tilde{b}}$ results into an enhancement of the ratio $A^-/A^+$.

3.2 ’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 \rightarrow \pi$ 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 \rightarrow \infty$, the latter keeps the ratio $N_f/N_c$ constant. From Eq. (20), it is then clear that $\tilde{a}$ and $\tilde{b}$ have the same scaling in the Veneziano limit (the same holds for $\tilde{c}$ and $\tilde{d}$). 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 $\tilde{a}-\tilde{d}$ coefficients.

4 $\Delta S=1$ amplitudes in Chiral Perturbation Theory

4.1 Chiral Dependence of the $K \rightarrow \pi$ amplitudes

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. [28, 46]. Here we just summarize the required formulæ, 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 $\bar{Q}^+$ and $\bar{Q}^-$ 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^\pm$, that need to be determined non-perturbatively:

$$\begin{aligned} \mathcal {H}_W^{ChPT} = g^+ \mathcal {O}^+ + g^- \mathcal {O}^-, \end{aligned}$$

(21)

with

$$\begin{aligned} \mathcal {O}^\sigma = \sum _{ijkl} c^\sigma _{ijkl} F^4 (U \partial _\mu U^\dagger )_{ij} (U \partial ^\mu U^\dagger )_{kl}, \end{aligned}$$

(22)

where U is the chiral meson field, i, j, k, l are flavour indices, and $c^\sigma _{ijkl}$ are Clebsch-Gordan coefficients (see Appendix A in Ref. [28]).

By means of the chiral weak Hamiltonian in Eq. (21) and the standard NLO ChPT Lagrangian, the chiral predictions for the normalized amplitudes in Eq. (16) are found to be:

$$\begin{aligned} A^\pm = g^\pm \left[ 1\mp 3 \left( \frac{M_\pi }{4 \pi F_\pi } \right) ^2 \left( \log \frac{M_\pi ^2}{\mu ^2} + L_\pm ^r(\mu ) \right) \right] , \end{aligned}$$

(23)

where $L^r_\pm$ 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:

$$\begin{aligned} g^\pm = \lim _{M_\pi \rightarrow 0} A^\pm . \end{aligned}$$

(24)

The extracted values of $g^\pm$ can then be used to make predictions of other observables, such as the $K \rightarrow \pi \pi$ 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:

$$\begin{aligned} \begin{aligned} g^\pm = 1 \pm a_\chi {1 \over N_c} \pm b_\chi {N_f \over N^2_c} + c_\chi {1 \over N^2_c}+ d_\chi {N_f \over N^3_c} + \cdots \end{aligned} \end{aligned}$$

(25)

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_c^0)$, and identical at this order. The next term in the $1/N_c$ expansion for $L^r_\pm$ could in principle differ:

$$\begin{aligned} L^r_\pm = L^{(0)} + \frac{1}{N_c} L^{(1)}_{\pm } + \cdots . \end{aligned}$$

(26)

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

$$\begin{aligned} \begin{aligned} A^+ A^-&= g^+ g^- \left[ 1 +3 \left( \frac{M_\pi }{4 \pi F_\pi } \right) ^2 {( L^r_- - L^r_+ )} \right] , \end{aligned} \end{aligned}$$

(27)

with

$$\begin{aligned} g^+ g^-= & {} 1 + \alpha {1 \over N_c^2} + \beta {1 \over N_c^3}+ \ldots , \end{aligned}$$

(28)

$$\begin{aligned} L^r_- - L^r_+= & {} \frac{ L^{(1)}_- - L^{(1)}_+ }{N_c} + \ldots , \end{aligned}$$

(29)

where $\alpha$ and $\beta$ depend on the coefficients $a_\chi - d_\chi$.

4.2 Relation to $K \rightarrow \pi \pi$ amplitudes

Once the effective couplings $g^\pm$ have been extracted from the chiral extrapolations of the ratios $A^\pm$, they can be used to compute the $K \rightarrow \pi \pi$ weak decay amplitudes. The two pions in the final state can be in a state with total isospin $I=0$ or 2:

$$\begin{aligned} iA_Ie^{i\delta _I} = {\langle {\left( \pi \pi \right) _I | \mathcal {H}_W^{ChPT} |K^0 }\rangle }, \end{aligned}$$

(30)

where $\delta _I$ is the two-pion scattering phase. The ratio of the two amplitudes can be calculated at leading order in ChPT using the Hamiltonian in Eq. (21) [28, 48]:

$$\begin{aligned} \frac{A_0}{A_2} = \frac{1}{2\sqrt{2}} \left( 1 + 3 \, \frac{g^-}{g^+} \right) . \end{aligned}$$

(31)

The measured hierarchy of $\sim 22$ between $A_0$ and $A_2$ must then be translated into a large ratio of the couplings $g^\pm$. Note that for $g^+ = g^- = 1$, the expected large-$N_c$ result is recovered, $A_0/A_2 = \sqrt{2}$. Large $1/N_c$ corrections in the $g^-/g^+$ ratio could therefore be the origin of the $\Delta 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⁷. 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:

$$\begin{aligned} \begin{aligned}&\text {Re } \frac{A_0}{A_2} \Big |_{M_\pi ,M_D \rightarrow 0, M^{\text {phys}}_K}= {\frac{1}{2\sqrt{2}} \left( 1 + 3 \, \frac{g^-}{g^+} \right) }\\&\quad {+ \frac{17}{12\sqrt{2}} \left( 1+ \frac{1}{17}\frac{g^-}{g^+} \right) \frac{M_K^2}{(4 \pi F_K)^2} \log \frac{\Lambda _{\mathrm{eff}}^2}{M_K^2}}, \end{aligned} \end{aligned}$$

(32)

where $\Lambda _{\mathrm{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 $\Lambda _{\mathrm{eff}} > M_K$ ($\Lambda _{\mathrm{eff}} < M_K$).

Table 3 Summary of the simulation parameters of the various ensembles used in this work

5 Lattice setup

5.1 Simulation and matching of sea and valence sectors

Our lattice setup is the same as the one presented in Ref. [27], 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 [50, 51]. We have chosen the Iwasaki gauge action (following previous experience with 2+1+1 simulations [52]) and clover Wilson fermions for the sea quarks, with the plaquette-boosted one-loop value of $c_{\mathrm{\scriptscriptstyle sw}}$. The simulation parameters are shown in Table 3. We find that a separation of $\ge 10$ units of Montecarlo time produces no autocorrelation in the ratios. The lattice spacing is found to be $a\sim 0.075$ fm for all values of $N_c$ (see also Ref. [27]). In addition, we have produced two ensembles with a finer lattice spacing, $a \sim 0.065$ fm, to estimate discretization effects.

Table 4 Summary of results for our ensembles with Iwasaki gauge action and O(a)-improved Wilson fermions with $c_\mathrm{\scriptscriptstyle sw}=0$ in the valence sector throughout. The value of the lattice spacing is $a \simeq 0.075\text { fm}$ for the “A” ensembles (see Ref. [27]), whereas it is $a \simeq 0.065\text { fm}$ for “B” ensembles. We provide the pion mass in the valence sector, $aM^\text {v}_\pi$, and the PCAC mass, $am^\text {v}_{\text {pcac}}$. We also include the results for the ratios in Eq. (13), and in the last column, the chiral parameter $\xi \equiv {M_\pi ^2}/{(4\pi F_\pi )^2}$. Moreover, $\xi _L$ labels $\xi$ corrected by finite-volume effects as explained in the main text

In order to achieve automatic O(a) improvement⁸ [55] and avoid the mixing of different-chirality operators for weak decays, we employ maximally twisted valence quarks [56], i.e., the mixed-action setup [57] previously used in Refs. [53, 54]. Working in twisted quark field variables, maximal twist is ensured by tuning the untwisted bare valence mass $m^\mathrm{v}$ to the critical value for which the valence PCAC mass is zero:

$$\begin{aligned} \lim _{m^\mathrm{v} \rightarrow m_{\mathrm{cr}}} m^\mathrm{v}_{\mathrm{pcac}} \equiv \lim _{m^\mathrm{v} \rightarrow m_{\mathrm{cr}}} \frac{\partial _0{\langle {A_0^{ij}(x) P^{ji}(y)}\rangle }}{2 {\langle {P^{ij}(x) P^{ji}(y)}\rangle }}= 0. \end{aligned}$$

(33)

The bare twisted mass parameter $\mu _0$ is tuned such that the pion mass in the sea and valence sectors coincide, $M_\pi ^\mathrm{v}= M_\pi ^s$.

Since twisted mass already provides O(a) improvement, the clover improvement parameter $c_{\mathrm{\scriptscriptstyle sw}}$ can be chosen to be an arbitrary value in the valence sector. We choose $c_{\mathrm{\scriptscriptstyle sw}}=0$ in the valence sector⁹ 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 non-perturbative determination in Ref. [58] for $N_c=3$ (see below). Finally, we also observe that $c_\mathrm{\scriptscriptstyle sw}=0$ leads to smaller statistical errors.

In Table 4 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$\sigma$. In addition, the valence and sea pion masses are matched also within 1 or 2$\sigma$. The bare results for the ratios are also presented in the same table, together with the chiral parameter $\xi = {M_\pi ^2}/{(4\pi F_\pi )^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. [6, 7]. Those results used quenched simulations, with plaquette gauge action and twisted mass fermions. The lattice spacing was $a\sim 0.093$ fm and the the pion mass was fixed at around $M_\pi = 550-590$ MeV for $N_c =3-8$ and 17. In this work, we perform a reanalysis of these quenched data.

5.2 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_\pi 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 $\xi$ they can be of $O(1\%)$ and thus we correct for them, as explained in Ref. [27], following Refs. [59, 60].

Since $B_K$ and $\bar{R}^+$ differ by a volume-independent proportionality factor, we can use the results in Ref. [61], where the finite-volume effects of $B_K$ have been calculated. In addition, it is known that the finite-volume and chiral corrections of $\bar{R}^+$ and $\bar{R}^-$ are fully anticorrelated [46]. Thus, we find:

$$\begin{aligned} \bar{R}^\pm (L) = \bar{R}^\pm \left[ 1 \pm 6 \sqrt{2 \pi } \xi \frac{e^{- M_\pi L}}{(M_\pi L)^{3/2}}(M_\pi L - 4) \right] . \end{aligned}$$

(34)

The correction for these quantities is numerically negligible for our ensembles. While additional finite-volume effects could be present (see Ref. [60]) 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 $\sim 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. [58] to the one-loop perturbation theory results in their setup (they used $c_{sw} =0$). The difference is roughly $\sim 5\%$ for $N_c=3$. On the other hand, we have computed the ratios using $c_{\mathrm{\scriptscriptstyle sw}}=1.69$ in the valence sector for the 3A10 ensemble. Using the perturbative renormalization constants for this new value of $c_{\mathrm{\scriptscriptstyle sw}}$ we get a result that differs from our $c_{\mathrm{\scriptscriptstyle 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_{\mathrm{\scriptscriptstyle sw}}$ (the perturbative one-loop corrections are also significantly larger for the larger value of $c_{\mathrm{\scriptscriptstyle 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.

6 Results

6.1 $N_c$ scaling of $K\rightarrow \pi$ amplitudes

The physical amplitudes $A^\pm$ can be obtained, as explained in Eq. (16), from the bare ratios in Table 4, and the renormalization coefficients in Tables 1 and 2. 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:

$$\begin{aligned} \begin{aligned} \frac{{A}^- + {A}^+}{2}&= 1 +\tilde{c} \, \frac{1}{N_c^2} + \tilde{d} \, \frac{N_f}{N_c^3} + \ldots ,\\ \frac{{A}^- - {A}^+}{2}&= -\tilde{a} \, \frac{1}{N_c}- \tilde{b} \, \frac{N_f}{N_c^2} + \ldots . \end{aligned} \end{aligned}$$

(35)

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 $\sim 570$ MeV.

• 2. Dynamical results ($N_f=4$) at a heavy pion mass $\sim 560$ MeV (ensembles A10).

• 3. Dynamical results ($N_f=4$) at a lighter pion mass $\sim 360$ MeV (ensembles A40).

The results for the coefficients $\tilde{a}-\tilde{d}$ for the three scenarios are presented in Table 5 and Fig. 4. The coefficients are all of ${{\mathcal {O}}}(1)$ and therefore of natural size. Importantly the sign of the $\tilde{a}$ and $\tilde{b}$ 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 that $\tilde{b}, \tilde{d} \sim {{\mathcal {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 to $\tilde{b}$ and $\tilde{d}$ being absent for $N_f=0$. The other two coefficients, $\tilde{a}$ and $\tilde{c}$, 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 $\tilde{a}$, 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.

Fig. 4

Half-sum and half-difference of the amplitudes $A^\pm$ as a function of $N_c^{-1}$ for three different cases: (i) quenched results from Ref. [6] in blue, (ii) new dynamical results at a pion similar to the quenched case (red), and (iii) dynamical results at a lighter pion mass (orange). The fit results are shown in Table 5. Error bars include only statistical errors

Table 5 Summary of results for the $1/N_c$ fits to the half-sum and half-difference of the amplitudes $A^\pm$. Errors are only statistical

6.2 Kaon B-parameter ($B_K$)

The kaon B-parameter, $B_K$, is defined from the matrix element of the $\Delta S=2$ operator that mediates neutral kaon oscillations at physical kinematics:

$$\begin{aligned} {\langle {\bar{K}^0| O^{\Delta S=2}(\mu ) | K^0}\rangle } = \frac{8}{3} f_K^2 M_K^2 {\bar{B}}_K(\mu ). \end{aligned}$$

(36)

It is customary to quote the renormalization group independent (RGI) version, labelled as $\hat{B}_K$. Its value at the physical point has been computed accurately in $N_f=2$, $2+1$, and $2+1+1$ simulations [58, 62,63,64,65,66] (see Ref. [67] for a review).

In our setup, $\hat{B}_K$ coincides with the renormalized ratio $\bar{R}^+$ up to a normalization. Specifically, we have

$$\begin{aligned} \hat{B}_K = \frac{3}{4} \hat{c}^+(a^{-1}) \bar{R}^+ \end{aligned}$$

(37)

where $\hat{c}^+$ can be read off Table 1. There are two essential differences in our setup: all meson masses are degenerate, in particular $M_K = M_\pi$, and we have an active light charm quark. Both can significantly affect the value of $\hat{B}_K$.

We show our results in Fig. 5. We observe a very significant $N_c$ dependence of $\hat{B}_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 of $\hat{B}_K$, while the $N_f=4$ result is about 25% smaller. We have included as bands the Buras-Bardeen-Gerard (BBG) Dual QCD prediction from Ref. [20], using inputs on meson masses from our own simulations in both cases — quenched and dynamical. We find that our results are reasonably compatible with the BBG prediction, in particular regarding the suppression of $\hat{B}_K$ in the presence of a light charm.

Fig. 5

Lattice results for $\hat{B}_K$, defined in Eq. (37), in the case of $N_f=0$ (see Refs. [6, 7]), and $N_f=4$ (this work). Error bars are only statistical errors. We also include the predictions from Ref. [20], where the band indicates the values obtained when varying the involved matching scale M from 600 to 1000 MeV

To conclude this subsection, we can use the scaling in $N_c$ to infer a value of ${\hat{B}}_K$ with three active flavours and quasi-physical kinematics. For this, we use the coefficients $\tilde{a}-\tilde{d}$ in Table 5 for the case of $N_f=4$ and $M_\pi =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 value $\hat{B}_K$ as in Eq. (37), extracting ${\bar{R}}_+$ and using the ${{\hat{c}}}^+(a^{-1})$ for three-flavour QCD¹⁰. We find

$$\begin{aligned} \hat{B}_K\big |_{M_K = M_\pi } = 0.67(2)_{\text {stat}} (6)_{Z^+} (3)_{\text {fit}} \, , \end{aligned}$$

(38)

including statistical error, and a $\sim 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 of ${\bar{R}}^\pm$ instead of $A^\pm$.

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 of $\hat{B}_K$ in the degenerate case, to the quasi-physical (QP) situation with $M_\pi =0$ and $M_K$ at its physical value:

$$\begin{aligned} \hat{B}^{QP}_K = \hat{B}_K\big |_{M_K = M_\pi } \left[ 1 + \frac{2}{3} \left( \frac{M_K}{4 \pi F_K}\right) ^2 \log \frac{\Lambda _{\mathrm{eff}}^{B_K}}{M_K} \right] , \end{aligned}$$

(39)

where $\Lambda _{\mathrm{eff}}^{B_K}$ labels an unknown scale that parametrizes the effect of the unknown LECs. For $\Lambda ^{B_K}_\mathrm{eff}>M_K$, ${\hat{B}}^{QP}_K$ is larger than $\hat{B}_K$ and could be compatible with the existing results at the physical point from $N_f=2+1,N_c=3$ simulations [58, 62,63,64,65,66].

6.3 Extraction of the effective couplings $g^\pm$

The main goal of this work is to compute the ratio $g^-/g^+$ by extrapolating $A^\pm$ to the chiral limit. For the required chiral extrapolation, we follow the same strategy as in Ref. [48]. 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

$$\begin{aligned} \frac{g^-}{g^+} \equiv \left( g^- g^+ \right) \times \frac{1}{(g^+)^2}. \end{aligned}$$

(40)

This approach results in a milder chiral extrapolation, that will hopefully introduce a smaller systematic error.

Table 6 Results for Fit 1: the simultaneous chiral and $N_c$ fits for $A^+$ and $A^+ A^-$. Errors are only statistical

Table 7 Results for Fit 2: the chiral fit at $N_c=3$ for $A^+$ and $A^+ A^-$. Errors are only statistical

Fig. 6

Chiral extrapolation of $A^+$ and the product $A^+ A^-$. The data points are also shown in Table 4. Empty squares for $N_c=3$ indicate a finer lattice spacing. Solid lines indicate a simultaneous chiral and $N_c$ fit as in Eq. (23). Dashed lines represent the chiral extrapolation of the data points for $N_c=3$ following Eqs. (23) and (27). Errors are only statistical

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:

$$\begin{aligned} \begin{aligned} \text {Fit 1: }&g^+ = 0.187(21),\ \ \ g^+ g^- = 0.91(4), \\ \text {Fit 2: }&g^+ = 0.190(27),\ \ \ g^+ g^- = 0.80(6). \end{aligned} \end{aligned}$$

(41)

The complete results of these fits are shown in Tables 6, and 7, and also in Fig. 6.

From these results, we obtain for the ratio of couplings at $N_c=3$:

$$\begin{aligned}&\frac{g^-}{g^+} \Bigg |_{\text {fit 1}} = 26(6), ~\frac{g^-}{g^+} \Bigg |_{\text {fit 2}} = 22(5), \end{aligned}$$

(42)

where errors are only statistical, but correlations are taken into account.

6.4 $K \rightarrow \pi \pi$ amplitudes in ChPT

Using the result for the ratio of couplings in Eq. (42), and the NLO ChPT prediction in Eq. (32), we can obtain an indirect result for the ratio of isospin amplitudes in the $K \rightarrow \pi \pi$ decay for $N_c=3$. In Fig. 7, we show this prediction as a function of an unknown effective scale $\Lambda _{\mathrm{eff}}$. This prediction, valid for $M_\pi =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:

$$\begin{aligned} \text {Re }\frac{A_0}{A_2} \Bigg |_{N_f=4} = 24(5)_{\text {stat}}(4)_{\text {fit}}(5)_{Z^\pm }(3)_{\text {NLO}}, \end{aligned}$$

(43)

where the central value comes from the fit 2 result in Eq. (42). 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 Sect. 5.2 —, and (iv) a $10\%$ error from the NLO effects — see Fig. 7. Combining all error sources in quadrature results in a $\sim 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 $\Delta I=1/2$ enhancement.

Fig. 7

NLO ChPT prediction (in red) for the ratio of $K \rightarrow \pi \pi$ isospin amplitudes as a function of the NLO LEC, $\Lambda _{\mathrm{eff}}$. We use the input of Fit 2 in Eq. (42). This prediction is valid for $M_\pi =M_D =0$, and $M_K$ 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)

7 Conclusions

We have presented the first non-perturbative study of the scaling of $\Delta 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\simeq 0.075$ fm and $a\simeq 0.065$ fm and pion masses in the range $360-570$ MeV. We have analysed the $K\rightarrow \pi$ amplitudes $A^\pm$, mediated by the two current-current operators $Q_\pm$ of the $\Delta S=1$ weak Hamiltonian in Eq. (1).

The diagrammatic analysis of the large-$N_c$ scaling of these observables presented in Sect. 3 allows to classify the subleading $N_c$ corrections, and demonstrates the anticorrelation of the leading ${{\mathcal {O}}}(1/N_c)$ and ${\mathcal O}(N_f/N^2_c)$ contributions in the $A^\pm$ 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 ${{\mathcal {O}}}(1)$. They can nevertheless explain the large enhancement of the ratio $A^-/A^+$ for $N_c=3$ with respect to the $N_c\rightarrow \infty$ limit. This involves an unprecedentedly large unquenching effect in this ratio, that is nevertheless compatible with natural size ${{\mathcal {O}}}(N_f/N_c^2)$ corrections.

The amplitudes $A^\pm$ in the chiral limit can be matched to their ChPT counterparts, which depend on the leading low-energy couplings, $g^\pm$, of the chiral effective weak Hamiltonian. From a chiral extrapolation of the combinations $A^+$ and $A^+ A^-$, we have then extracted the couplings $g^\pm$, which are finally used to predict in ChPT the ratio of $K \rightarrow (\pi \pi )_{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 “$\Delta I=1/2$ rule”. Our estimate for this ratio in the theory with a light charm is

$$\begin{aligned} \text {Re }\frac{A_0}{A_2} \Bigg |_{N_f=4} = 24(5)_{\mathrm{\scriptscriptstyle stat}}(7)_{\mathrm{\scriptscriptstyle sys}}, \end{aligned}$$

(44)

which suggests that the enhancement may indeed be largely dominated by intrinsic QCD effects.