$K \to \pi$ matrix elements of the chromomagnetic operator on the lattice

We present the results of the first lattice QCD calculation of the $K \to \pi$ matrix elements of the chromomagnetic operator $O_{CM} = g\, \bar s\, \sigma_{\mu\nu} G_{\mu\nu} d$, which appears in the effective Hamiltonian describing $\Delta S = 1$ transitions in and beyond the Standard Model. Having dimension 5, the chromomagnetic operator is characterized by a rich pattern of mixing with operators of equal and lower dimensionality. The multiplicative renormalization factor as well as the mixing coefficients with the operators of equal dimension have been computed at one loop in perturbation theory. The power divergent coefficients controlling the mixing with operators of lower dimension have been determined non-perturbatively, by imposing suitable subtraction conditions. The numerical simulations have been carried out using the gauge field configurations produced by the European Twisted Mass Collaboration with $N_f = 2+1+1$ dynamical quarks at three values of the lattice spacing. Our result for the B-parameter of the chromomagnetic operator at the physical pion and kaon point is $B_{CMO}^{K \pi} = 0.273 ~ (70)$, while in the SU(3) chiral limit we obtain $B_{CMO} = 0.072 ~ (22)$. Our findings are significantly smaller than the model-dependent estimate $B_{CMO} \sim 1 - 4$, currently used in phenomenological analyses, and improve the uncertainty on this important phenomenological quantity.

where Q + γ,g (Q − γ,g ) are the parity-even (-odd) EMO and CMO, respectively, defined as: with q R,L = 1 2 (1 ± γ 5 ) q (for q = s, d). In Fig. 1   Note that at least one mass insertion is required in the diagrams, both in the SM and beyond, in order to induce the chirality flip described by the magnetic operators.
A quick inspection of the diagrams of Fig. 1 shows that the Wilson coefficients of the magnetic operators in the SM and NP model are proportional to where M N P represents the typical NP scale, e.g. the gluino mass in the SUSY case, and the factors m s /M W and δ LR are generated in the diagrams by the mass insertion. In the SUSY case, for instance, δ LR represents the off-diagonal matrix element of the squark mass matrix normalized to the average squark mass. The transition rate is controlled in the SM by the weak coupling α W (M W ). This is not generally the case for NP models. In the SUSY transition shown in Fig. 1, for example, the process is mediated by the strong interactions. Therefore, the proportionality of C N P γ,g in Eq.
(3) to the strong coupling constant α s (M N P ), rather than to the weak coupling as in the SM, compensates in part for the stronger high-energy scale suppression (M N P > M W ) in the NP model. Thus, the magnetic interactions receive potentially large contributions from physics beyond the SM.
It is also worth noting that the chirality flipping factor m s /M W , which appears in the Wilson coefficients of the magnetic operators in the SM, is of the same size as Λ QCD /M W , which represents the additional suppression factor of the coefficients of dimension-6 operators in the effective Hamiltonian. For this reason, the role of the magnetic operators tends to be marginal in the SM, while it is potentially more relevant for the searches of NP.
The K → π matrix element of the EMO Q + γ , which is relevant for instance for the CP violating part of the rare K L → π 0 + − decays [1], has been computed on the lattice both in the quenched [2] and unquenched N f = 2 [3] case. Since the electromagnetic field strength tensor F µν factorizes out of the hadronic matrix element, the lattice computation only involves the quark bilinear operator s σ µν d, and it is relatively straightforward.
Computing the hadronic matrix elements of the CMO Q ± g is, instead, by far more challenging. The main difficulty is represented by the complicated renormalization pattern of the operator, which also involves power divergent mixing with operators of lower dimensionality (see Ref. [4] and Section III). The relevant matrix elements with an initial kaon involve one, two or three pions in the final states, and are of great phenomenological interest for various processes: the long distance contribution to K 0 −K 0 mixing [5], ∆I = 1/2, K → ππ transitions and ε /ε [1], CP violation in K → 3 π decays [6]. These matrix elements are parameterized in terms of suitably defined , where the low-energy constant B CM O is estimated to be of order 1 in the chiral quark model of Ref. [7]. Therefore, the three B-parameters of Eqs. (4)(5)(6) are related by SU(3) chiral symmetry, which predicts at LO their equality: Such an equality is expected to be broken at higher orders in ChPT.
In this work we evaluate the B-parameter appearing in Eq. (4) from the lattice QCD computation of the π|Q + g |K matrix element. We perform numerical simulations by employing the gauge configurations generated by the European Twisted Mass Collaboration (ETMC) with N f = 2+1+1 dynamical quarks, which include in the sea, besides two light mass-degenerate quarks, also the strange and charm quarks with masses close to their physical values [8][9][10]. The same gauge ensembles have been used in Ref. [4] to determine non-perturbatively the power divergent mixing coefficients controlling the mixing of the CMO with operators of lower dimension. As for the mixing coefficients with operators of the same dimensionality and for the multiplicative renormalization constant (RC), we adopt the predictions of perturbation theory at one-loop obtained in Ref. [4].
Preliminary results for the B-parameter of the CMO have been presented in Ref. [11] and our final result at the physical pion and kaon point is where the first error comes from the numerical lattice simulations, while the second error accounts for the perturbative uncertainty in the one-loop determination of the multiplicative renormalization constant. Our result (8) represents the first lattice QCD determination of a matrix element of the CMO. In the SU(3) chiral limit we get B CM O = 0.076 (23). Our findings are significantly smaller than the model-dependent estimate B CM O ∼ 1 − 4 currently used in phenomenological analyses [1] and improve the uncertainty on this important phenomenological quantity.
The plan of the paper is as follows. In section II we describe the lattice setup and give the simulation details. In section III we recall the main results obtained in Ref. [4] on the determination of the power-divergent mixing coefficients needed for the renormalization of the CMO. The mixing subtraction is evaluated in section IV and it is shown that the renormalized CMO correlator can be determined with a remarkable level of precision. The lattice data for the matrix elements of the renormalized CMO are presented in section V and analyzed in terms of both SU(2) and SU (3) ChPT. Finally, section VI contains our conclusions.

II. SIMULATION DETAILS
Our lattice setup is based on the gauge configurations generated by ETMC with N f = 2 + 1 + 1 dynamical quarks [8,9], adopted in Ref. [10] for the determination of the up, down, strange and charm quark masses, using the experimental value of the pion decay constant f π to set the lattice scale 1 .
The gauge fields are simulated using the Iwasaki gluon action [12], while sea quarks are implemented with the Wilson Twisted Mass Action at maximal twist [13][14][15]. In order to avoid the mixing of strange and charm quarks induced by lattice artifacts in the unitary twisted mass formulation, we have adopted the mixed action setup described in Ref. [16], where the valence strange quarks are regularized as Osterwalder-Seiler (OS) fermions [17], while the valence up and down quarks have the same action as the sea. The use of different lattice regularisations for the valence and sea quarks of the second generation preserves unitarity in the continuum limit and brings no complications for the operator renormalization pattern in mass-independent schemes, while producing only a modification of discretization effects. Therefore, the uncertainty related to the use of a non-unitary action at finite lattice spacing is incorporated directly in the error due to discretization effects, which will be addressed in Section V. Moreover, since we work with both valence and sea quarks at maximal twist, physical observables are guaranteed to be automatically O(a)-improved [15,16].
The details of the ETMC gauge ensembles with N f = 2+1+1 dynamical quarks are collected in     Quark propagators are obtained using the multiple mass solver method [18,19], which allows to invert the Dirac operator for several quark masses at a relatively low computational cost. The statistical accuracy of the meson correlators is significantly improved by using the "one-end" stochas-tic method [20], which includes spatial stochastic sources at a single time slice chosen randomly.
Statistical errors are evaluated using the jackknife procedure.

III. RENORMALIZATION OF THE CHROMOMAGNETIC OPERATOR
In this Section we briefly review the main results obtained in Ref. [4] for the renormalization of the CMO, whose specific renormalization pattern depends on the details of the lattice regularization, i.e. on the choice of the lattice action. A detailed analysis of the implications of the discrete symmetries of the twisted-mass action was carried out in Ref. [4], showing that the renormalization of the CMO involves the mixing among 13 operators of equal or lower dimensionality 2 , including also non gauge invariant operators vanishing by the equation of motion.
In the case of on-shell matrix elements the mixing simplifies, and the renormalized parity-even CMO can be written as [4] O where O CM ≡ 16π 2 Q + g = gsσ µν G µν d, S =sd, P = isγ 5 d and O 4 = (sd) are bare operators, valence quarks are taken with the same value of the Wilson r-parameter, i.e. r s = r d (see Ref. [4]), and µ s (µ d ) denotes the bare strange (light) quark mass.
Note that the quadratically divergent mixing of the CMO with the scalar density S is common to any regularization, whereas the mixing with the pseudoscalar density P (softened by the proportionality to the quark masses) is peculiar of twisted mass fermions, and it is a consequence of the non conservation of parity. Moreover, in Eq. (9) the power divergent mixing coefficients c 12 and c 13 are scheme and renormalization scale independent [22], while the multiplicative RC Z CM and the coefficients c i with i = 2, 3, 4 depend on both the scheme and the renormalization scale.
It is well known [23] that the determination of power divergent coefficients, controlling the mixing with operators of lower dimension, cannot rely on perturbation theory. The reason is that potential non-analytic (in g 2 ) contributions to these coefficients, like those proportional to powers of (1/a) exp(−1/(β 0 g 2 )) ∼ Λ QCD , do not appear in the perturbative expansion. Therefore, while for the present study the multiplicative renormalization factor Z CM and the coefficients c i with i = 2, 3, 4 have been evaluated in perturbation theory at one-loop, the coefficients c 12 and c 13 in Eq. (9) have been determined in a non-perturbative way by imposing two suitable subtraction 2 The operator mixing pattern considered in Ref. [4] includes one operator (O6), which is not independent from the other ones, and misses one five-dimensional operator [21], which mixes with the CMO only at two loops in perturbation theory. The results presented in this work are therefore not affected.
conditions [4]. The first one is that the matrix element of the CMO between external kaon and pion at rest must vanish in the SU(3) chiral limit, namely from which the coefficient c 13 can be determined. The second requirement is the vanishing of the parity violating matrix elements of the CMO up to terms of O(a), specifically from which the coefficient c 12 can be calculated once the coefficient c 13 is determined from Eq. (10).
In Table II we present the numerical results for the various mixing coefficients, obtained in Ref. [4] at the three values of the inverse coupling β given in Table I. The central values and the errors of the coefficient c 13 shown in the last column correspond to the averages and the spread of the two non-perturbative determinations corresponding to the choice "LP" given in Table IV of Ref. [4]. For the mixing coefficient c 12 the uncertainty of the non-perturbative results has been found at the level of 60% [4], and therefore the g 2 -dependence of the non-perturbative determination of c 12 , shown in the penultimate column, can be safely neglected. GeV, and of the mixing coefficients c i in Eq. (9), obtained in Ref. [4]. The results correspond to the three values of the inverse coupling β given in Table I using one-loop perturbation theory, except for c 12 and c 13 in the last two columns, which have been obtained non-perturbatively (see text). The perturbative results have been evaluated using the bare coupling g 2 0 = 6/β for the power divergent coefficients c 12 and c 13 and the boosted coupling g 2 It can be seen that: i) the power divergent coefficient c 13 has been determined nonperturbatively with a very high level of precision; ii) the coefficients c 3 and c 4 vanish at one loop; iii) the coefficient c 2 which starts at O(g 2 ) is rather small; and iv) the multiplicative renormalization factor Z CM receives at one loop a sizable correction (∼ 70%).
In Table II we also provide the one-loop results for the power divergent coefficients c 12 and c 13 .
For the latter the difference between the one-loop and the non-perturbative results is less than 10%. The bulk of the difference is compatible with being a correction of O(g 4 ). Thus, genuine non-perturbative contributions to c 13 are likely to be small, even though a firmer conclusion in this sense would require the calculation of c 13 at two loops at least.
As far as the coefficient c 12 is concerned, its size is smaller by (at least) one order of magnitude with respect to c 13 both perturbatively and non-perturbatively (see Table II). In addition, the corresponding operator is proportional to the first power of the quark masses, with a(m s + m d ) ∼ 0.02 in our simulations. For these reasons, the subtraction of the linear divergence in Eq. (9) has a numerically negligible impact on the determination of the CMO matrix elements (see next Section).

IV. LATTICE QCD CORRELATORS
In order to evaluate the matrix elements of the renormalized CMO (9), using the values of the mixing coefficients given in Table II,  operator O CM the gluon tensor G µν is replaced by its clover discretization P µν , namely [24] O CM = g 0 ψ s σ µν P µν ψ d , where and the sum is over the four plaquettes P j (x) in the µ-ν plane stemming from x and taken in the counterclockwise sense.
The K → π matrix elements of the bare local operators O i (i = {CM, S, P }) are extracted from the large (Euclidean) time distance behavior of a convenient combination of 2-and 3-point correlation functions, which for both initial and final mesons at rest are defined as where t is the time distance between the source and the sink, t is the time distance between the insertion of the operator O i and the source, while P K (x) = is(x)γ 5 u(x) and P π (x) = id(x)γ 5 u(x) are the local interpolating fields of the K and π mesons, respectively. The Wilson parameters r of the two valence quarks in both initial and final mesons are always chosen to have opposite values, i.e. r s = r d = −r u , so that the squared meson masses differ from their continuum counterpart only by terms of order O(a 2 mΛ QCD ) [15].
The statistical accuracy of the correlators (14-15) is significantly improved by using the all-to-all quark propagators evaluated with the so-called "one-end" stochastic method [20], which includes spatial stochastic sources at a single time slice chosen randomly. Statistical errors are evaluated using the jackknife procedure.
At large time distances 2-and 3-point correlation functions behave as where Z π(K) is the matrix element 0|P π(K) (0)|π(K) and M π(K) is the mass of the π(K) meson.
Both quantities are determined adopting the fitting function (16)  The matrix elements K|O i |π can be extracted from the time behavior of the following ratios where s i (t, t ) is the sign of correlator C Kπ i (t, t ) and the correlation function C π(K) (t) is given by which at large time distances behave as i.e. without the backward signal. At large time distances one has so that the bare matrix elements K|O i |π can be calculated from the plateau of R i independently of the matrix elements Z π and Z K of the interpolating fields. In order to minimize excited state effects, in what follows the source-sink separation is fixed to t = T /2. Therefore, the region of time distances, where both the initial and final ground states dominate leading to the plateau (21), corresponds to [t min , T /2 − t min ]. Such a time interval will be adopted to extract the CMO matrix elements.

The size of the various terms involved in the renormalization of the CMO (see Eq. (22) below)
and expressed in lattice units, namely a 3 R CM (t, T /2), c 13 aR S (t, T /2), c 12 a(µ s +µ d )aR P (t, T /2) and c 2 a 2 (µ 2 s +µ 2 d )aR S (t, T /2), can be inferred from Fig. 2, where the results refer to β = 2.10 and quark masses aµ d = 0.0020 and aµ s = 0.0151 (corresponding to M π 255 MeV and M K 520 MeV).
The values adopted for the mixing coefficients c 12 and c 13 are the non-perturbative ones, while the value of the coupling entering c 2 is taken from boosted perturbation theory (see Table II). It can be seen that the ratio a 3 R CM (t, T /2) and the one corresponding to the leading power divergence, i.e. c 12 a(µ s + µ d )aR P (t, T /2) and c 2 a 2 (µ 2 s + µ 2 d )aR S (t, T /2), are smaller by almost five orders of magnitude. The operators µ s µ d S and O 4 = (sd) do not contribute at one loop, since their mixing coefficients vanish (see Table II).
In the lower panel of Fig. 2 we show the renormalized CM ratio R CM (t, T /2), given by Despite being the outcome of a large numerical subtraction, due almost totally to the power divergent term (c 13 /a 2 )R S , the results for R CM (t, T /2) are clearly different from zero and exhibit plateaux, from which the renormalized CMO matrix element K| O CM |π can be determined quite precisely (see next Section). In this respect we stress that the high level of precision achieved in the non-perturbative determination of c 13 (see the last column of Table II) plays a crucial role.
Notice that the renormalized CMO ratio a 3 R CM (t, T /2) is two orders of magnitude larger than the mixing term c 2 a 2 (µ 2 s + µ 2 d )aR S (t, T /2). Since c 2 (as well as c 3 and c 4 in Eq. (9)) is known only at one loop in perturbation theory, higher-order corrections and non-perturbative effects might contribute to the the renormalized CMO ratio away from the SU(3) chiral limit. However, the smallness of the mixing term c 2 a 2 (µ 2 s + µ 2 d )aR S (t, T /2) indicates that higher orders and nonperturbative effects are not expected to play a significant role in the determination of the matrix elements of the renormalized CMO within the present statistical uncertainties.
Note that the evaluation of the B-parameter (23) does not require the knowledge of the lattice spacing, while it involves the mass RC Z m , which in our maximally twisted-mass setup is given by where Z P is the RC of the pseudoscalar density. For the latter we adopt the RI'-MOM results obtained in Ref. [10] using the two methods M1 and M2, which differ by O(a 2 ) effects.
Besides B Kπ CM O we have calculated two further B-parameters. They correspond to the transitions induced by the CMO in which either the mass m s of the strange valence quark is taken to be equal to the light-quark mass m d or the mass m d of the light valence quark is taken to be equal to the strange quark mass m s . In both cases the spectator valence quark, which does not participate in the transition, is always a up-quark with mass m ud (see Eqs. (14)(15)). We will refer to the above two transitions as the ππ and KK channels, respectively. Explicitly one has  Table I in the MS(2 GeV) scheme.
Since we have simulated three values of the strange quark mass around its physical value (see Table I channels, the results for the B-parameters exhibit controllable discretization effects, which are the beneficial consequence of the subtraction of the power-divergent mixings (thanks to the use of maximally twisted Wilson quarks the subtracted CMO matrix elements vanish in the chiral limit even at finite lattice spacing). Moreover, the impact of finite volume effects can be estimated by comparing the results corresponding to the ensembles A40.24 and A40.32, which share common values of the pion mass and the lattice spacing, and differ only by lattice size L. No significant effects are visible within the statistical errors.
As described in section I, the ChPT prediction at LO is that all three B-parameters should coincide and be independent of the light-quark mass. The latter feature is approximately fulfilled only in the case of π → π channel (see upper panel in Fig. 3). Chiral corrections beyond the LO are clearly visible for both the Kπ and KK channels (see middle and bottom panels in Fig. 3), since they exhibit a remarkable dependence on the light-quark mass and deviate strongly from the results corresponding to the ππ channel.
In order to extrapolate to the physical pion mass we need to take into account at least NLO effects, which however are not known analytically. We observe, in this respect, that the chiral corrections for the CMO matrix element contain also powers of p K · p π . This means that for the B-parameter B Kπ CM O NLO terms proportional not only to M 2 K and M 2 π , but also to M K M π should be considered. Note that the quantity M K M π ∝ m 1/2 ud is non-analytic in the light-quark mass. In the ππ and KK channels such a non-analytic term is not expected.
Thus, we introduce the variables where B 0 and f 0 are LO low-energy constants (LECs) of the SU (3) ChPT [25], and we consider for the generic channel i → j with i(j) = K, π the following SU(3)-inspired Ansatz where B CM O is the LO LEC appearing in Eq. (7) (i.e., the SU(3) chiral limit of the B-parameters), while the parameters α 1 , β 1 and γ 1 play the role of NLO LECs, and α 2 , β 2 , β 2 and γ 2 are NNLO LECs. In Eq. (27) the two terms proportional to ξ π and ξ 2 π are due to the dependence on the mass of the u-quark (not involved in the transition) in common for all the channels. Note that the NNLO term proportional to ( where • () stat+f it indicates the uncertainty induced by both the statistical errors and the fitting procedure itself; 4 In Eq. (27) an additional NNLO term proportional to (ξi + ξj)( √ ξi − ξj) 2 can be considered. We have checked that the impact of its inclusion is almost negligible. Moreover, the inclusion of additional discretization terms proportional either to a 2 ξπ or to a 2 ( √ ξi − ξj) 2 produces no significant effect within the errors.
• () chir corresponds to the uncertainty related to the chiral extrapolation, obtained using the results corresponding to the inclusion (γ 2 = 0) or the exclusion (γ 2 = 0) of the NNLO non-analytic term in Eq. (27); • () disc is the uncertainty related to discretization effects estimated by adding a term proportional to a 4 without any prior; • () Z P is the error induced by the use of the two methods M1 and M2 to obtain the mass RC Z m = 1/Z P in Ref. [10].
As a further check, we have also analyzed separately the data for the Kπ channel (see also Ref. [11]) obtained after a smooth interpolation at the physical value of the strange quark mass m phys s = 99.6(4.3) MeV [10]. This allows us to adopt the following SU(2)-inspired Ansatz where the parameters α, β and γ play the role of SU (2) where the second error accounts for the perturbative uncertainty in the one-loop determination of the multiplicative RC Z CM (see Table II) and it has been estimated to be 25% relying on the Our result (31) represents the first lattice QCD determination of a matrix element of the CMO.
In the SU(3) chiral limit we obtain B CM O = 0.076 (14) (18) P T = 0.076 (23). Our findings are significantly smaller than the model-dependent estimate B CM O ∼ 1 − 4 currently adopted in phenomenological analyses [1]. The comparison indicates also that the uncertainty on this important phenomenological quantity has been significantly reduced by lattice QCD. A further drastic improvement in the precision can be achieved by removing the uncertainty due to the use of one-loop perturbation theory for estimating the multiplicative RC Z CM , which is by far the dominating source of uncertainty in our results.

VI. CONCLUSIONS
We have presented the results of the first lattice QCD calculation of the K → π matrix elements of the chromomagnetic operator O CM = gs σ µν G µν d, which appears in the effective Hamiltonian describing ∆S = 1 transitions in and beyond the Standard Model.
Having dimension 5, the chromomagnetic operator is characterized by a rich pattern of mixing with operators of equal and lower dimensionality. The power divergent coefficients controlling the mixing with operators of lower dimension were determined non-perturbatively in Ref. [4], while the multiplicative renormalization factor as well as the mixing coefficients with the operators of equal dimension have been computed at one-loop in perturbation theory [4].