Kaon Mixing Beyond the SM from Nf=2 tmQCD and model independent constraints from the UTA

We present the first unquenched, continuum limit, lattice QCD results for the matrix elements of the operators describing neutral kaon oscillations in extensions of the Standard Model. Owing to the accuracy of our calculation on \Delta S=2 weak Hamiltonian matrix elements, we are able to provide a refined Unitarity Triangle analysis improving the bounds coming from model independent constraints on New Physics. In our non-perturbative computation we use a combination of Nf=2 maximally twisted sea quarks and Osterwalder-Seiler valence quarks in order to achieve both O(a)-improvement and continuum-like renormalization properties for the relevant four-fermion operators. The calculation of the renormalization constants has been performed non-perturbatively in the RI-MOM scheme. Based on simulations at four values of the lattice spacing and a number of quark masses we have extrapolated/interpolated our results to the continuum limit and physical light/strange quark masses.


Introduction
The fundamental target of present-day research activity in Particle Physics is the search for New Physics (NP) effects beyond the Standard Model (SM) predictions. Two main routes are followed, one based on the study of processes in which the direct production of NP particles at high energy colliders like the LHC takes place, and a second one based on the indirect investigation of NP effects coming from the exchange of virtual NP particles.
In the so-called indirect approach a crucial role is played by Flavor Physics processes that are sensitive to NP through loop effects. These processes vanish at tree level in the SM, and some of them are theoretically very clean, despite the fact that they are loop mediated, and in some cases also CKM or helicity suppressed. Among them, ∆F = 2 transitions have always provided some of the most stringent constraints on NP. For instance, the constraints from K 0 −K 0 oscillations are particularly stringent for NP models that generate transitions between quarks of different chiralities [1]- [4]. Therefore, an accurate determination of the ∆S = 2 bag parameters (B-parameters) is crucial to the improvement of NP constraints.
In the present work, we provide the first accurate lattice determination of the ∆S = 2 B-parameters relevant for physics beyond the SM, calculated in the continuum limit and using data from unquenched, N f = 2, dynamical quark simulations 2 . Our results represent a significant improvement with respect to the (quenched) input values of Refs. [6] and [7] used so far in phenomenological analyses. The calculation of the B K parameter that is relevant for the K 0 −K 0 mixing in the SM has been presented in [8] using three values of the lattice spacing. In the present work we update that value by adding a fourth (finer) lattice spacing. We note that the difference between the two results is about half standard deviation.
The outline of the paper is as follows. Section 2 contains a brief description of the ∆S = 2 matrix elements of the effective weak Hamiltonian describing the most general pattern of K 0 −K 0 oscillations. In Section 3, based on the results of this work for the ∆S = 2 B-parameters, we discuss the implications for NP of our updated Unitarity Triangle (UT) analysis [9]. In Section 4 we illustrate the main theoretical features of the lattice setup employed in our simulations (twisted mass lattice QCD [10], [11]) and we describe the strategy for obtaining accurate numerical estimates of the B-parameters as well as ratios of kaon four-fermion matrix elements. In Section 5 we collect our numerical results. In Section 6 we give our estimates of the various B-parameters and matrix elements ratios and compare the present results with the previous determinations existing in the literature. In five Appendices we discuss a number of technicalities: i) the renormalization properties of the four-fermion operators in our "mixed action" setup [12]; ii-iii) the RI-MOM computa-tion of renormalization constants (RCs) and corresponding results respectively; iv) tables of lattice data on pseudoscalar meson masses, decay constants and bare four-fermion matrix elements; v) complete results on renormalized four-fermion matrix elements obtained by using various formulae for the chiral extrapolation and two alternative procedures for the RI-MOM determination of RCs.

∆S = effective weak Hamiltonian
The general form of the ∆S = 2 effective weak Hamiltonian is where in the so-called SUSY basis ( [1], [13]) the four-fermion operators O i andÕ i have the form with α and β denoting color indices. Spin indices are implicitly contracted within square brackets. The Wilson coefficients C i andC i have an implicit renormalization scale dependence which is compensated by the scale dependence of the renormalization constants of the corresponding operators 3 . Notice that the parity-even parts of the operators O i andÕ i are equal. From now on and for notational simplicity we will denote by O i (i = 1, . . . , 5) the parity-even components of the operators (2.2). Due to parity conservation in strong interactions, in the study of K 0 − K 0 oscillations it is then sufficient to consider only the matrix elements K 0 |O i |K 0 . We recall that in the SM only the kaon matrix element of the operator O 1 comes into play.
The bag parameters, B i (i = 1, . . . , 5), provide the value of four-fermion matrix elements in units of the magnitude of their vacuum saturation approximation. More explicitly they are defined by the equations [14] for i = 2, . . . , 5, (2.5) with ξ i = (8/3, −5/3, 1/3, 2, 2/3). In the above relations one recognizes B 1 as the familiar B K . We recall that, as suggested by the parametrization adopted in the r.h.s. of Eqs. (2.4) and (2.5), the K 0 |O 1 (µ)|K 0 matrix element is expected to vanish in the chiral limit unlike the other four ones. An alternative way which has the merit of allowing a more accurate evaluation of matrix elements, is to consider the matrix elements ratios R i = K 0 |O i |K 0 / K 0 |O 1 |K 0 , i = 2, . . . , 5, as first proposed in Ref. [6]. For details on our lattice implementation see Section 4 and in particular Eq. (4.19).
For the reader's convenience we here anticipate our final continuum results for B i and R i in the MS scheme of Buras et al., defined in Ref. [15], and the RI-MOM scheme 4 at 2 GeV, see Tables 1 and 2 respectively. Details on the calculation and uncertainty estimates are given in Section 5. In Appendix E we also give the final continuum results for B i and R i in the MS and the RI-MOM scheme at 3 GeV.

MS (2 GeV)
0.52 (2) 0.70(2) 1.22 (7) 1.00(4) 0.69(5) 1 -12.9(4) 4.5(2) 21.2(7) 4.7(3) 3 Model-independent constraints on ∆S = 2 operators and New Physics scale from the Unitarity Triangle analysis ∆F = 2 processes provide some of the most stringent constraints on NP generalizations of the SM. Several phenomenological analyses of ∆F = 2 processes have been performed in the last years, both for specific NP models and in model-independent frameworks. A generalization of the UT analysis, which allows for NP effects by including the most significant flavour constraints on NP available at the time was performed in Ref. [9]. The result was a simultaneous determination of the CKM parameters and the size of NP contributions to ∆F = 2 processes in the neutral kaon and B d,s meson sectors. The NP generalization of the UT analysis consists in including in the theoretical parametrization of the various observables the matrix elements of operators which, though absent in the SM, may appear in some of its extensions. The analysis shows that the constraints coming from K 0 −K 0 matrix elements are the most stringent ones, in particular for models that generate transitions between quarks of different chiralities (see Refs. [1] - [4]). Thus an accurate determination of the ∆S = 2 B-parameters is crucial to the improvement of the NP constraints.
The results for the ∆S = 2 B-parameters obtained in the present work come from unquenched N f = 2 lattice QCD data carefully extrapolated to the continuum limit. They hence represent a significant progress with respect to the input values used in the UT analysis performed in Ref. [9], where quenched lattice numbers (without a systematic continuum limit extrapolation analysis) computed more than five years ago in Refs. [6] and [7] were employed. For this reason, we present here an update of the analysis of Ref. [9] based on our new values of the ∆S = 2 B-parameters. The new ingredients entering the analysis are collected in Tables 1 and 2. For all the other input data we use the numbers quoted in Ref. [19] in the Summer 2012 analysis.
In the present NP-oriented analysis, the relations among experimental observables and the CKM matrix elements are extended by taking into consideration the most general form of the ∆S = 2 effective weak Hamilonian (see Eq. (2.1)). The effective weak Hamilonian is parameterized by Wilson coefficients of the form where F i is the (generally complex) relevant NP flavor coupling, L i is a (loop) factor which depends on the interactions that generate C i (Λ), and Λ is the scale of NP, i.e. the typical mass of new particles mediating ∆S = 2 transitions. For a generic strongly interacting theory with an unconstrained flavor structure, one expects F i ∼ L i ∼ 1, so that the phenomenologically allowed range for each of the Wilson coefficients can be immediately translated into a lower bound on Λ. Specific assumptions on the flavor structure of NP correspond to special choices of the F i functions. For example Minimal Flavor Violation (MFV) models [20]- [25] correspond to F 1 = F SM and F i =1 = 0. Following Ref. [9], in deriving the lower bounds on the NP scale Λ, we assume L i = 1, that corresponds to strongly-interacting and/or tree-level coupled NP. Two other interesting possibilities are given by loop-mediated NP contributions proportional to either α 2 s or α 2 W . The first case corresponds for example to gluino exchange in the minimal supersymmetric SM. The second case applies to all models with SM-like loop-mediated weak interactions. To obtain the lower bound on Λ entailed by loop-mediated contributions, one simply has to multiply the bounds we quote in the following by α s (Λ) ∼ 0.1 or α W ∼ 0.03.
In agreement with Ref. [9], we find that in the K 0 sector, due to the non-vanishing chiral limit (chiral enhancement) of their matrix elements, all bounds coming from the contributions of non-standard operators (i.e. from the operators O i with i = 1) are more than one order of magnitude stronger than the bound from the SM O 1 operator.
The results for the upper bounds on the Im C K i coefficients and the corresponding lower bounds on the NP scale Λ are collected in Table 3 where they are compared to the previous results of Ref. [9]. The superscript K is to recall that we are reporting the bounds coming from the kaon sector we are here analyzing. Although several input parameters have been updated with respect to Ref. [9] (see Ref. [19]), the more stringent constraints on the Wilson coefficients of the non-standard operators and, consequently, on the NP scale, mainly come from the improved accuracy achieved in the values of the ∆S = 2 B-parameters obtained in the present work. This can be realized by comparing the small improvement of the bound coming from Im C K 1 , obtained using a value of the B 1 -parameter very close to the one taken in Ref. [9], with those coming from the other coefficients using the new B-parameters. We observe that the analysis is performed (as in [9]) by switching on one coefficient at the time in each sector, thus excluding the possibility of having accidental cancellations among the contributions of different operators. Therefore, the reader should 95% allowed range Lower limit on Λ (GeV −2 ) (TeV)  Table 3: 95% probability range for the Im C i K coefficients and the corresponding lower bounds on the NP scale, Λ, for a generic strongly interacting NP with generic flavor structure (L i = F i = 1). In the lower panel the results of [9] are displayed for comparison. keep in mind that the bounds may be weakened if, instead, some accidental cancellation occurs.
In Fig. 1 we show the comparison between the lower bounds on the NP scale obtained for the case of a generic strongly interacting NP with generic flavor structure by the constraints on the Im C K i coefficients coming from the present generalized UT analysis, and the previous results of Ref. [9].
As a specific example of NP models we consider the warped five-dimensional extensions of the SM discussed in Ref. [26], where the origin of hierarchies in quark masses and mixings is explained via the localization properties of quark wave functions in the fifth dimension. In particular, in the Randall-Sundrum (RS) scenario one has where M G and g * s ∼ 6 are the mass and coupling of Kaluza-Klein excitations of the gluon, Y * ∼ 3 is the five-dimensional Yukawa coupling (whose flavour structure is assumed to be anarchic), m d ∼ 3 MeV and m s ∼ 50 MeV are MS quark masses at the high scale and v = 246 GeV is the Higgs vev. Running from a reference scale of 5 TeV, we obtain at 95% probability Im C K 4 ∈ [−4.7, 10.6] · 10 −18 , from which we get The lower bounds on the NP scale, provided by the constraints on Im C K i (i = 1, . . . , 5) for generic NP flavor structure, are shown as brown bars. For comparison, we plot the bounds of Ref. [9] as yellow bars.
Considering instead gauge-Higgs unification (GHU) models, from Ref. [26] we have 4) where in this case g * ∼ 4 is the five-dimensional gauge coupling in units of the radius of the compact dimension. We obtain the bound

Non-perturbative lattice computation of the ∆S = 2 matrix elements
Lattice QCD provides an ideal first principle framework in which non-perturbative computation of hadronic matrix elements can be performed with controlled systematic uncertainties. In the last few years a number of lattice determinations of B K of increasing precision have appeared in the literature based on a variety of lattice regularizations with N f = 2 or N f = 3 dynamical fermions [8], [27]- [31]. For recent reviews see Refs.
[32]- [34]. Very little has been done in the literature concerning the calculation of the physical matrix elements of the full operator basis of Eq. (2.2). Calculations of the whole set of ∆S = 2 renormalized operators have been carried out, in the quenched approximation, using improved Wilson fermions (see Refs. [6] and [14]) or the chirality conserving overlap and domain-wall fermions (see Refs. [7] and [36], respectively). Very recently a N f = 2+1 dynamical quark calculation appeared [5] that uses domain-wall fermions at one value of the lattice spacing 5 .
The first calculation of the kaon matrix elements of the whole operator basis was performed employing Clover improved Wilson fermions. Since the Clover term coefficient was set to its tree-level value, matrix elements were affected by O(g 2 0 a) discretization errors. Simulations were carried out at two values of the gauge coupling corresponding to lattice spacings ∼ 0.07 and ∼ 0.09 fm ( [6], [14]). The major source of systematic errors was, however, the uncertainty related to the construction of the multiplicatively renormalizable lattice operators O i . In fact, owing to the breaking of the chiral symmetry intrinsic in the Wilson fermion action, the bare counterparts of each of them mix with all the other operators of equal dimension including those with "wrong" chiral transformation properties [37]. All the mixing coefficients and the overall RC were computed in the non-perturbative RI-MOM scheme [16].
Similar quenched computations were carried out using overlap and domain-wall fermions (see Refs. [7] and [36], respectively). Though performed at pretty coarse lattice spacings (namely a ∼ 0.09 and ∼ 0.13 fm in the case of overlap fermions and a ∼ 0.1 fm in the case of domain-wall fermions), these simulations have the advantage that the renormalization properties of the operators entering the four-fermion basis are as in the continuum, and lattice artifacts are O(a 2 ). Also in this case the non-perturbative RI-MOM scheme was used in the computation of the various RCs.
In the following sections we present a new unquenched computation of the kaon matrix elements of the full ∆S = 2 four-fermion operator basis employing lattice data from simulations with N f = 2 dynamical fermions, performed at four rather fine values of the lattice spacing in the interval [0.05, 0.1] fm. We are thus able to safely extrapolate the lattice estimators of all the relevant matrix elements to the continuum limit (CL).
We use a mixed fermion action setup where we adopt different regularizations for sea and valence quarks. In particular we introduce maximally twisted (Mtm) sea quarks [11] that we take in combination with Osterwalder-Seiler (OS) [38] valence quarks. This strategy has been suggested in Ref. [12] as a way of setting up a computational framework allowing for a calculation of ∆S = 2 four-fermion matrix elements that is both automatically O(a) improved and free of wrong chirality mixing effects. A proof of the latter point is given in Appendix A and its validity is numerically verified in Appendix B, while O(a) improvement of physical quantities is a genuine property of the setup of Ref. [12]. As a consequence unitarity violations due to different sea and valence quark regularization yield only O(a 2 ) artifacts, provided renormalized sea and valence quark masses are matched. In our case, the matching of the renormalized quark masses is obtained by simply taking identical values for the corresponding sea and valence bare mass parameters.
The interesting lattice setup briefly described above has already been successfully tested in B K computations both in the quenched approximation [39] and on ensembles with N f = 2 [8] and N f = 2 + 1 + 1 dynamical quarks [40], as well as in unquenched (N f = 2 + 1 + 1) studies of meson masses and decay constants [41] and nucleon sigma terms [42].

Sea and valence quark regularization
The Mtm-LQCD action of the light quark flavor doublet can be written in the so-called "physical basis" in the form [11] (4.1) The subscript sea is to remind us that this action will be used to generate unquenched gauge configurations. The field ψ describes a mass degenerate up and down doublet with bare (twisted) mass µ sea . The parameter M cr is the critical mass that one has to fix nonperturbatively at its optimal value (as proposed in Refs. [43]- [45] and implemented in Refs. [46] and [47]) to guarantee the O(a)-improvement of physical observables and get rid of all the unwanted leading chirally enhanced cutoff effects. In the gauge sector the tree-level improved action proposed in Ref. [48] has been used.
For valence quarks we use the OS regularization [38]. The full valence action is given by the sum of the contributions of each individual valence flavour q f and reads [12] S OS val = a 4 x,fq where the index f labels the valence flavors and M cr is the same critical mass parameter which appears in Eq. (4.1). We denote by r f and µ f the values of the Wilson parameter and the twisted quark mass of each valence flavor.

Lattice operators and correlation functions
In the strategy proposed in Ref. [12], which we follow here, four species of OS valence quark flavors (q f , f = 1, . . . , 4) are introduced, two of which (q 1 and q 3 ) will represent the valence strange quark with masses µ 1 = µ 3 ≡ µ "s" , while the other two (q 2 and q 4 ) will be identified with the light up/down quarks having masses µ 2 = µ 4 ≡ µ ℓ . The corresponding r f Wilson parameters must obey the relation In the numerical computations reported in the present work we have averaged over the two cases r 1 = ±1, holding r 2 , r 3 and r 4 related to r 1 as in Eq. (4.3).
As in the case of the computation of B K (≡ B 1 ), in the calculation of the bag parameters B i (i = 2, . . . , 5), we need to consider the axial currents and the pseudoscalar quark densities In addition, we need to consider the following set of four-fermion operators where square parentheses denote spin invariants and α and β are color indices. In the mixed action (MA) approach defined above the following properties can be proved (see Appendix A and Ref. [12]).
The key statement (iii) follows by noting that the set of fermionic Wick contractions contributing to the three-point correlation functions (see Eq. • the relevant renormalization constants in the RI-MOM scheme can also be evaluated with no O(a) artefacts.
The first of these properties was derived in Ref. [12] for generic hadron masses and matrix elements 6 . The second property above follows from the remark that the one-particleirreducible vertices entering the RI-MOM renormalization conditions are O(a) improved in our MA setup, just because such vertices turn out to be invariant under parity transformations of their external momenta. The argument here is closely analogous to the one presented for the renormalization constants of quark bilinear operators in the Appendix of Ref. [18].
In the construction of correlation functions we follow the general procedure outlined in Ref. [8]. We use periodic boundary conditions in every direction for all fields, except for the quark fields on which we impose anti-periodic boundary conditions in the time direction. At time slices y 0 and y 0 +T /2 "wall" operators with K 0 -meson quantum numbers are inserted.
The first operator is constructed in terms ofq 2 and q 1 quark fields and the second in terms ofq 4 and q 3 quark fields. Explicitly they have the expressions yq 4 ( y, y 0 + T /2)γ 5 q 3 ( y, y 0 + T /2) (4.10) The correlators we then need to compute are 7 To improve the signal-to-noise ratio a sum has been performed over the spatial position of the four-fermion operator, and for each gauge configuration the time slice y 0 is randomly chosen. An important contribution to the reduction of statistical fluctuations comes also from summing over the spatial position of both kaon interpolating fields in Eq. (4.11). After summing over the the spatial position of the four-fermion operator, which is known from experience to be crucial for the signal, to project on zero 3-momentum states, just one further spatial sum is needed. The second spatial sum, which we are able to do, gives a further signal improvement. These spatial sums were implemented and carried out at a reasonably low computational price by means of the stochastic technique discussed in sect. 2.2 of Ref. [8].
For large time separation y 0 ≪ x 0 ≪ y 0 + T /2 the (plateau of the) following ratio estimate , i = 2, . . . , 5 (4.14) provides an estimate of the B  8 In the following the superscript (b) denotes bare quantities.

Estimates of renormalized quantities
Recalling Eqs. (4.7) and (4.9), the renormalized values of the bag parameters will be given by the formula where a sum over the repeated index j is understood. As for the renormalized expression of the four-fermion operator ratios of Eq. (4.17), we choose to evaluate the rescaled quantitỹ where M lk and F lk are the mass and decay constant of the pseudoscalar meson made of the (q l q k ) quark pair. In order to compensate the chiral vanishing of the K 0 |O 1 |K 0 matrix element we have multiplied the ratios R i by the factor M 12 M 34 /F 12 F 34 ; its form has been chosen in such a way to partially compensate the lattice artifacts affecting the different lattice discretizations of kaon mesons (resulting from different choices of the OS r f -parameters) we use. Furthermore, we have remultiplied the quantity in the square bracket by the ratio of the experimental values of the kaon decay constant (f exp K = 156.1 MeV) and its mass (m exp K = 494.4 MeV). The definition we get in this way is in line with the one proposed in Ref. [7]. Based on the discussion of Section 4.2 (in particular item iii), we note that in the continuum limit and at the physical values of the u/d and s quark masses the quantityR i of Eq. (4.19) provides the right estimate for the ratio of the renormalized matrix elements of interest, i.e.
We end this section by recalling that in the twisted mass mixed action setup of Ref. [12] we are using, lattice estimators of physical quantities are only affected by O(a 2 ) lattice artifacts. This is true in particular for the kaon mass and decay constant as well as for the kaon four-fermion matrix elements.

Simulations, data analysis and results
The ETM Collaboration has generated N f = 2 configuration ensembles at four values of the inverse bare gauge coupling, β and at a number of light quark masses, µ sea . The values of the simulated lattice spacings lie in the interval [0.05, 0.1] fm. Bare quark mass parameters are chosen so as to have light pseudoscalar mesons ("pions") in the range 280 ≤ m PS ≤ 500 MeV and heavy-light pseudoscalar mesons ("kaons") in the range 450 ≤ m P S ≤ 650 MeV. Simulation details are given in Table 4.
The value of the light u/d quark mass parameter, aµ ℓ , is common to sea and valence quarks, while the heavier quark (the would-be strange quark that we denote by "s", see Table 4) is quenched. As discussed in Section 5.2, we will get to the physical kaon mass by suitably interpolating (extrapolating) data in µ "s" (µ ℓ ) to the "physical" value µ s (µ u/d ), while simultaneously taking the continuum limit. The "physical" values µ u/d and µ s of the quark masses are known and can be found in Ref. [49]. The quark bilinear RCs, Z P and Z S , have been computed in the non-perturbative RI-MOM scheme in Ref. [18]. A RC computation for the full basis of the four-fermion operators using RI-MOM techniques is presented in Appendix B. In Appendix C we collect the values of the four-fermion RCs that are used in this work, as well as Z P and Z S .

Extracting bare estimates from lattice data
Bare results for the ratio of the four-fermion matrix elements R

Computation at the physical point
Extracting physical quantities from lattice data requires performing extrapolations and/or interpolations of renormalized lattice estimators to the physical point (continuum limit and "physical" value of quark masses). β = 3.80, a ∼ 0.10 fm

RCs computation and combined continuum-chiral extrapolation
We have computed the full matrix of the four-fermion operator RCs in a mass independent scheme. We carry out the non-perturbative calculation adopting the RI-MOM approach. The implementation of the RI-MOM setup has been presented in Refs. [8] and [18]. We should mention that in our RC estimators cutoff effects, though parametrically of O(a 2 ), are numerically reduced owing to the subtraction of perturbatively evaluated O(a 2 g 2 ) contributions. After that, two different, but by now standard [18], procedures are employed to deal with O(a 2 p 2 ) discretization effects. The first, called M1, consists in linearly extrapolating to zero the residual (after the perturbative subtraction) O(a 2 p 2 ) terms. The second one (socalled p 2 -window method, or M2 for short) leads to RC estimates obtained by averaging data over a fixed (in physical units) and very narrow momentum interval. We carry out continuum and chiral extrapolations in a combined way. For all bag parameters, B i , and ratios,R i (see Eq. (4.19)), we have tried out a fit ansatz of the following 2τ /T 2τ /T 2τ /T where we have made explicit the dependence of the fit parameters A (n) Y and D Y on the renormalized strange quark mass 9 in units of r 0 (r 0μs ). We studied separately the cases 2τ /T 2τ /T 2τ /T 2τ /T of linear and polynomial ansatz. We have also considered NLO ChPT fit functions for B i based on the formulae given in [50] in the case of SU(3). Those formulae transformed to NLO SU(2) ChPT read: withB 0 = 2.84(11) GeV (renormalized in MS at 2 GeV) and f 0 = 121.0(1) MeV, as we used in [8]. The sign before the logarithmic term is minus (-) for i = 1, 2, 3 and plus (+) for i = 4, 5. As forR i and i = 2, 3 the ChPT fit formula at NLO coincides with the linear fit ansatz, while for the cases i = 4, 5 we usẽ (7) [49]. Their values in the MS scheme at 2 GeV are In the four panels of Fig. 4 we show the combined chiral and continuum fit (see Eq. (5.1) and Eq. (5.3)) for the ratiosR i against the renormalized light quark mass for i = 2, . . . , 5, respectively. The RCs used in these plots are the ones computed with the M1method and are expressed in the MS scheme of Ref. [15] at 2 GeV. Lattice data correspond to points taken at the pair of quark masses (r 0μℓ , r 0μs ).
In panels (a) and (b) (corresponding to cases with i = 2, 3 respectively) we display the curves that correspond to the polynomial fit function (5.1) at the four β values we are considering in this paper. The black solid line represents the continuum limit curve. The dashed black line represents the continuum limit curve that is obtained if a linear fit ansatz in µ ℓ is used. Black open circles and triangles stand for the results at the physical quark mass point from the polynomial and the linear fit ansatz, respectively. Recall that forR i with i = 2, 3 ChPT fit formula at NLO coincides with a linear fit ansatz. In panels (c) and (d) (corresponding to cases with i = 4, 5 respectively) we display the curves corresponding to the ChPT fit formula and the linear fit function. In this case black open circles and triangles stand for the results at the physical quark mass point from the ChPT fit and the linear fit ansatz, respectively.
Similarly, in Fig. 5 we present the combined chiral and continuum fit for the B iparameters, again renormalized in the MS scheme of Ref. [15] at 2 GeV. In all four panels we display the curves corresponding to the ChPT fit formula at the four β values. The black solid line represents the continuum limit curve. The dashed black line represents the continuum limit curve that is obtained if a linear fit ansatz inμ ℓ is used.
Note the nice agreement (within one standard deviation) of the two fit ansätze for both theR i ratios and the B i parameters. show the behaviour vs. the renormalized light quark mass of the combined chiral and continuum fits (according to the polynomial formula (5.1) with n = 2) of theR i ratios, with i = 2 and i = 3 respectively, renormalized in the MS scheme of Ref. [15] at 2 GeV with the M1-type RCs. The full black line is the continuum limit curve. In panels (c) and (d), solid lines, instead, show the combined chiral and continuum fit described by NLO-ChPT, Eq. (5.3) for i = 4 and i = 5, respectively. The full black line is the continuum limit curve. The dashed black line represents the continuum limit curve in the case of the linear fit ansatz. Black open circles and triangles stand for the results at the physical point corresponding to the polynomial (panels (a) and (b)) and ChPT fit (panels (c) and (d)), and linear fit ansatz, respectively.
In Tables 1 and 2 of Section 2 we have gathered our final continuum results for R i and B i in the MS of Ref. [15] and RI-MOM scheme at 2 GeV respectively. The final value of B i for i = 2, . . . , 5 has been computed by averaging the estimates obtained from the three 2)) for the B i parameters with i = 2, . . . , 5 respectively, renormalized in the MS scheme of Ref. [15] at 2 GeV with the M1-type RCs. The full black line is the continuum limit curve (5.1). The dashed black line represents the continuum limit curve in the case of the linear fit ansatz. Black open circles and triangles stand for the results at the physical point corresponding to the ChPT fit and linear fit ansatz, respectively. kinds of fit ansatz discussed above, and using bootstrap error analysis. The half difference between the two more distant results has been taken as an estimate of the systematic error associated to the extrapolation procedure. The total uncertainty is obtained by adding in quadrature the statistical and the systematic error. For i = 1 we update the result for B K published in Ref. [8]; note that the difference between the two results is about half standard deviation.
In Appendix E we provide more detailed results obtained from the various fitting procedure we have investigated. We also show that the continuum extrapolated quantities that are eventually obtained by employing M1-type and M2-type RCs turn out to be perfectly consistent between each other within statistical errors. Also in Appendix E, see Tables 23  and 24, we quote our continuum results for B i and R i in the MS and the RI-MOM scheme respectively at 3 GeV.
As already stated, an alternative (indirect) way to compute the ratio of the kaon matrix elements of the renormalized operators O i , i = 2, . . . , 5 to that of O 1 is based on the formula (see eqs. (2.4) and (2.5)) This of course requires knowledge of the B i parameters and the renormalized quark masses. We find that the two evaluations (indirect and direct, based on Eq. (5.6) and Eq. (4.19), respectively) lead to compatible results within errors. However, the indirect estimates suffer from much larger final uncertainties. This is due to several reasons. One is related to the quadratic dependence on the quark mass, which makes the relative error on the mass to give a significant contribution to the final error. Furthermore in the indirect method one has to consider extra uncertainties due to the error of the bilinear operators' RCs that are used to compute the B-parameters. A comparison of the direct and indirect results obtained for the ratios R i is provided in Appendix E.

Conclusions
Accurate measurements of the K 0 −K 0 mixing amplitudes can yield useful hints on New Physics if theory can provide comparatively accurate calculations of quantities parametrizing beyond the SM effects. This requires a precise, first principle, evaluation of the kaon matrix elements of the full basis of four-fermion operators entering the most general effective ∆S = 2 weak Hamiltonian.
In this paper we have presented the first unquenched lattice QCD determination in the continuum limit of the matrix elements of the full ∆S = 2 four-fermion operator basis.
We have used N f = 2 unquenched tm-LQCD gauge configurations produced by the ETM Collaboration in combination with maximally twisted valence quarks of the OS type.
The mixed action setup proposed in Ref. [12] offers the possibility of obtaining automatically O(a) improved results and an operator renormalization pattern identical to that of a chirally invariant regularization at the rather cheap price of mere O(a 2 ) unitarity violations. Using data at four lattice spacings (with a in the interval [0.05, 0.1] fm) and a number of pseudoscalar masses ("pions") in the range [280, 500] MeV, we are able to safely carry out the continuum and the light quark mass limit of the observables of interest. All results are non-perturbatively renormalized in the RI/MOM scheme.
We get in this way the most accurate estimates to date of ∆S = 2 effective weak Hamiltonian matrix elements. The total error on the R i ratios is between 4% and 6% and on the bag parameters, B i , between 3% and 7%. Tables 5 and 6 show a comparison between our results for R i and B i (in RI/MOM at 2 GeV) and the data at fixed lattice spacings coming from the two old quenched calculations of Refs. [6] and [7] 10 .
For the B-parameters, one finds large differences between the central values of our results and those of Refs. [6] and [7], which vary between 5% and 25% (though the errors are typically comparably large). With respect to Ref. [6], the differences are even larger when the results are compared in terms of the ratios R i , presumably due to a combined effect, in this case, of having overestimated the values for both B 1 and the strange quark mass in the computation of [6]. We emphasize that, with respect to the old quenched calculations, having performed in the present study simulations at four values of the lattice spacing and quite smaller values of the pion masses provides us with a much better control over the main sources of systematic uncertainties, besides the quenched approximation. Current experience suggests that the possible systematic errors related to the quenching of the strange and charm quarks, which still affect our calculation, are negligible within the present uncertainties. Forthcoming results from simulations in the continuum limit with N f = 2 + 1 and N f = 2 + 1 + 1 dynamical flavours will provide a check of this expectation. We should add that our continuum limit results for R i and B i (i = 2, . . . , 5) are in the same ballpark with the numbers given at one lattice spacing in Ref. [5] where N f = 2 + 1 dynamical quarks are employed.
As an interesting phenomenological application of the results obtained in this paper we have carried out a new UT analysis along the lines of the work of Ref. [9]. Thanks to the improved accuracy of the present determination of the ∆S = 2 B-parameters, we could substantially strengthen the existing upper bounds on the Wilson coefficients of the operators of the non-standard sector of the effective weak Hamiltonian, and consequently increase the lower bound on the New Physics scale.  (2) 0.87 (7) 0.90(10) 0.67 (7) 0.72(9) 3 1.22 (7) 1.41 (12 Table 5 for the R i ratios.

A Renormalization properties of ∆S = 2 four-fermion operators
In this appendix we want to spell out the renormalization properties of the four-fermion operators of interest for the description ofK 0 − K 0 oscillations in the mixed action (MA) lattice setup of Section 4. We will do this by exploiting the results of Ref. [51]. In particular, we show that the operators in Eq. (4.6), that we report here for the reader convenience, are taken between the pseudoscalar operators P 12 =q 1 γ 5 q 2 and P 43 = q 4 γ 5 q 3 , one gets the same Wick contraction multiplicities one would obtain in QCD upon evaluating the kaon matrix elements of the operators (2.2). Naturally, apart from the issue of renormalization, the physical matrix elements will be obtained (in the continuum limit) by finally setting in our MA setup µ 1 = µ 3 = µ s and µ 2 = µ 4 = µ ℓ , with µ s (resp. µ ℓ ) corresponding to the bare strange (resp. degenerate up-down) quark mass.
A key result of Ref. [51] (stated there in the quark basis that is most natural for un- where σ µν = [γ µ , γ ν ]/2, one gets In Eq. (A.2) we have omitted color indices as they are always contracted within each square parenthesis.
In order to make direct contact with the formulae of Ref. [51] we must pass from the q f -basis, in which the valence quark action (4.2) was written and where the Wilson term is (maximally) twisted, to the χ f -basis, where the Wilson term takes its standard form. This is achieved by the chiral transformation x,fχ and assuming r 4 = ±1 (as used in the present work), According to Ref. [51], for the renormalized operatorsQ M A i one gets Since at µ f = 0 the fermion actionS OS val is indistinguishable from a standard massless Wilson fermion action, in any mass independent renormalization scheme, the operators Q , From the renormalizability of (correlation functions evaluated in) the MA lattice setup of Section 4 and the exact conservation of the individual valence flavours it immediately follows that the operator renormalization pattern of eqs. (A.8)-(A.10) is independent (up to cutoff effects, as usual) from the values of sea and valence quark masses. It is hence possible to determine the relevant renormalization constants in any mass-indipendent renormalization scheme by extrapolating to the chiral limit suitable renormalization constant estimators evaluated at non-vanishing quark masses. Following this strategy we computed non-perturbatively in the RI-MOM scheme the renormalization matrix Z Q , see Eq. (A.10), as detailed in Appendix B and summarized in Appendix C.

B RI/MOM computation of renormalization constants of four-fermion operators
In order to convert our lattice results for the bag parameters B i to their physical continuum counterparts, in the same renormalization scheme and at the same scale as the corresponding perturbative Wilson coefficients used in the phenomenological analysis, we need the renormalization constants (RCs) of the operators O M A i[+] , i = 1, 2, . . . , 5, see Eq. (4.6), or equivalently Eq. (A.1). As discussed in Section 4, in our mixed action (MA) setup for lattice correlation functions these operators represent the analogs of the parity-even parts of the ∆S = 2 four-fermion operators (2.2) that are relevant in the formal continuum theory.
In this Appendix, we give details on the non-perturbative computation of the RCs performed using the RI'-MOM scheme ( [17], [18] , i = 1, 2, . . . , 5, defined in Eq. (A.2). To lighten our notation, in the following we will drop the superscript and the sign subscript, denote these operators simply by Q i and assemble them in the array Q. The generic renormalization pattern of the bare operators Q (b) is of the form where the scale-dependent renormalization matrix Z is block-diagonal, with a continuumlike block structure (the same as for e.g. the matrix in Eq. (A.7)), while ∆ ∆ ∆ is a sparse off-diagonal and scale-independent matrix of the form However, as shown in Appendix A, using the MA lattice setup of Section 4, the "wrong chirality mixing" terms ∆ ij , are reduced to mere O(a 2 ) effects, the renormalization matrix Z coincides with the matrix Z Q of Eq. (A.10) and we recover a continuum-like renormalization pattern. This is a very important advantage of our approach, which we will implement in practice using the following strategy: compute the quark propagators in the q f -basis (also called physical basis of tmQCD at maximal twist, in which the critical Wilson term is twisted, see Eq. (4.2)), impose RI-MOM renormalization conditions on the operators Q i and extract the renormalization matrix (Z) and, for check purposes, the mixing matrix (∆ ∆ ∆).

B.1 Procedure for extracting the RCs
To determine the matrices Z and ∆ ∆ ∆ in Eq. (B.1) we proceed as follows. We start by computing the lattice quark propagator and the four-point Green functions with an insertion of the operator Q i , namely a 16 The lower(upper) case Greek (Latin) symbols denote uncontracted spin (color) indices. The corresponding amputated Green functions are given by For the sake of clarity, we will use matrix notation, denoting the matrices by boldface symbols and omitting color and spin indices. The amputated Green functions will be collected in the 1 × 5 row vector Λ Λ Λ(p) = ( Λ 1 , Λ 2 , Λ 3 , Λ 4 , Λ 5 ) (p, p, p, p) . In Eq. (B.8), Z q is the quark field RC and ∆ ∆ ∆ is, as we said before, the mixing matrix. We have also introduced the 5 × 1 column vector of spin projectors (see Eq. (37) of Ref. [51] for the explicit form of these projectors) which act on the amputed Green functions by (i, j = 1 · · · , 5), where the trace is taken over spin and colour, and obey the orthogonality relations j the tree level amputated Green function of the operator Q j . It is convenient to express Λ Λ Λ in terms of a "dynamics" matrix D, defined by , we see that, once the dynamics matrix is known, we can determine both the renormalization and the mixing matrices from the relation This matrix equation can be solved for Z and ∆ ∆ ∆ by exploiting the block diagonal structure of the Z matrix. In fact, it is easy to see that the three diagonal blocks of the renormalization matrix (i, j = 1; i, j = 2, 3 and i, j = 4, 5) are given by whereas the mixing coefficients are easily obtained from the equations, We can now summarize our procedure to determine the renormalization matrix of the parity-even part of the four-fermion operators of the SUSY basis of Eq. (2.2).
Step 1 The Green functions (B.4) and (B.5), and from them the dynamics matrix D, are evaluated in the Landau gauge for a sequence of sea, µ sea , and valence, µ val , quark mass values at each of the four lattice spacings we consider here. The bare parameters and the statistics of this computation are detailed in Table 2 of Ref. [18]. One can also find there (see Eqs. (3.6) and (3.7)) the set of discrete lattice momenta, p ν (p 1,2,3 = (2π/L) n 1,2,3 , p 4 = (2π/T ) (n 4 + 1/2)), that we include in the present calculation.
To minimize the contributions of Lorentz non-invariant discretization artifacts, we take into consideration only momenta satisfying the constraint In the following, we shall often use the short-handp 2 = νp 2 ν .
Step 2 For each β and each choice of the scalep 2 , the renormalization relation (B.14) is enforced at all values of µ sea and µ val given in Table 2 of Ref. [18]. By doing so, we obtain at nonzero quark masses the estimators Z RI ′ ij (p 2 ; a 2p2 ; µ val ; µ sea ), which are then extrapolated to µ val = 0 (see Section B.2) and µ sea = 0 (see Section B.3).
Step 4 Using the NLO continuum QCD evolution of the renormalization matrix Z calculated in Refs. [52,15], the first argument of Z RI ′ ij (p 2 ; a 2p2 ; 0; 0) is brought to a reference scale µ 2 0 . In this step, we assume that the scalesp 2 and µ 2 0 are large enough to make NLO perturbation theory accurate. This is the same level of accuracy achieved in the determination of the Wilson coefficients.
Step 6 In order to reduce the statistical error, the lattice RC estimators are averaged over two equivalent patterns of Wilson parameters (r 1 , r 2 , r 3 , r 4 ), namely (1, 1, 1, −1) and (−1, −1, −1, 1), as well as over different lattice momenta corresponding to the samẽ p 2 . We have checked that performing these averages before or after taking the chiral limit leads to consistent results.

B.2 Valence chiral limit
In view of the relation (B.14) and since the extraction of Z q (see Ref. [18]) poses no particular problems, our discussion will be mainly focused here on the quark mass dependence of the dynamics matrix D. At fixed values of β, a 2p2 and aµ sea , we fit the dynamics matrix elements D ij to the ansatz Here we have introduced a term with a pole in µ val ∼ m 2 P S to cope with the expected Goldstone boson (GB) pole contribution to the elements of the D matrix. The existence of such a GB-pole term can be understood as follows. At asymptotically large p 2 , nonperturbative effects giving contributions potentially divergent in the chiral limit to the Green functions (B.4) do vanish and the latter turn out to be polynomial in the quark mass parameters [16]. At finite values of p 2 , however, the contributions to (the spectral decomposition of) these Green functions from one-GB intermediate state with momentum q and mass m P S , give rise to terms proportional to (q 2 + m 2 P S ) −1 and suppressed by some power of 1/p 2 . If several one-GB intermediate states contribute to the spectral representation of the Green functions (B.4) several terms, each behaving as (q 2 + m 2 P S ) −1 and suppressed by some power of 1/p 2 , will show up. These results follow straightforwardly from the "polology' study of the Green functions (see e.g. the discussion in the book [53]) or from the well known Lehmann-Symanzik-Zimmermann (LSZ) reduction formalism. Now, since in the Green functions (B.4) the four-fermion operator is inserted at zero four-momentum transfer (q = 0), one expects, from the time orderings where two quark fields can create from the vacuum a pseudoscalar (i.e. GB) one-particle state, a contribution proportional to 1/m 2 P S , suppressed by some power of 1/p 2 . Similarly, from those time orderings where two quark fields create a GB-state and two further quark fields destroy another GB-state, contributions do arise that behave as (1/m 2 P S ) 2 and are twice more strongly suppressed at large p 2 . In conclusion, by exploiting (along the lines of Appendix A of Ref. [16]) the largep 2 behaviour of the matrix element 11 0|q f (p)q f ′ (−p)|P f ′ f and taking also into consideration the four factors of S q (p) −1 that stem from the relations (B.5), one finds that the dynamics matrix D ij contains GB-pole contributions of the following kinds: We recall that the kinematics of our Green functions corresponds to an exceptional momentum configuration where p 1 = p 3 = −p 2 = −p 4 = p, and thus q = 0. As in the chiral limit S −1 q (p) ∼ γ µ p µ , the result (B.18) implies that single and double GB-pole terms are suppressed by 1/p 2 and 1/(p 2 ) 2 factors, respectively. A second observation is 11 Here q ′ f (p) (q f (p)) denotes the four-dimensional Fourier transform of the quark field q ′ f (x) (q f (x)), while |P f ′ f is the pseudoscalar meson state with valence quarks of flavour f and f ′ . that thanks to the choice r 4 = −r 3 in our MA setup the lattice axial currentq 4 γ µ γ 5 q 3 is conserved (only broken by soft mass terms) and hence the matrix elements of the operator because, owing to r 2 = r 1 , the lattice axial currentq 2 γ µ γ 5 q 1 is broken by discretization effects (see Ref. [11] and Appendix A of Ref. [8]).
For the case j = 1, when no similar GB-pole simplifications can occur, double GBpole terms strongly suppressed (like 1/(p 2 ) 2 ) at large p 2 are to be expected in D ij (p). However, precisely owing to this strong suppression in practice, within our statistical errors and in the ranges of quark masses and p 2 we explore (see Table 2 of Ref. [18] and sect. B.6), we hardly see in our lattice data any effects that can reliably be ascribed to double GB-pole contributions. On the contrary, we do find clear numerical evidence for single GB-pole contributions, which indeed at high p 2 are only suppressed as 1/p 2 . We thus decided to ignore double GB-pole terms in our valence mass chiral extrapolations.
This choice is also justified a posteriori by the results of the valence chiral fits based on the ansatz (B.17). A subset of these results is illustrated in Fig. 6. There we display typical examples of the effect of GB-pole subtractions in the matrix elements of the dynamics matrix at two values of β. As can be seen, after the subtraction, a smooth dependence upon µ val (or equivalently on M 2 12 ) is observed. Combining the valence chiral limit lattice estimator of D ij and Z q , we are able to get reliable estimates of the intermediate quantities Z lat ij (p 2 ; a 2p2 ; 0, aµ sea ).

B.4 Removal of O(a 2 g 2 ) cutoff effects
We will obtain improved chiral limit RC estimators, Z RI ′ −impr ij (p 2 ; a 2p2 ), by removing from our Z RI ′ ij (p 2 ; a 2p2 ; 0, 0) lattice data perturbative discretization errors. This can be done up to O(a 2 g 2 ) exploiting the one-loop perturbative results obtained [54,55] in the massless lattice theory for the quark propagator form factor Σ 1 , related to the quark-field RC by Z q (p) = Σ 1 (p), and the dynamics matrix elements, v.i.z.

B.5 Absence of wrong chirality mixings
In Fig. 9, one can clearly see that for all the operators of interest the mixing coefficients ∆ ij are very small (in fact vanishing within errors in the range ofp 2 that we eventually use for extracting RCs). We also find that this is systematically more and more so as β increases, well in line with our expectation that in our lattice setup wrong chirality mixing effects are reduced to mere O(a 2 ) artifacts. For these reasons the effects of ∆ ∆ ∆ have been neglected in our final RC analysis, where we have assumed a fully continuum-like relation between renormalized and bare operators. In addition we checked that repeating the whole analysis with the tiny effects of ∆ ∆ ∆ on the relation (B.1) properly taken into account leads to no significant changes in the values of RCs.

B.6 Final RC estimates from M1 and M2 method
Having extrapolated the (improved) RCs estimators to the valence and sea chiral limit at each value of the momenta, we evolve Z RI ′ −impr ij (p 2 ; a 2p2 ) from the scalep 2 to a common scale µ 2 0 by using the known matrix formula for the NLO running of the operators Q i [52,15], obtaining Z RI ′ −impr ij (µ 2 0 ; a 2p2 ). This step is necessary in order to disentangle the O(a 2p2 ) cutoff effects from the genuine continuum p 2 dependence. Notice also that the actual value of µ 0 has no impact on the RGI results of the RC's. As is customary, we take µ 0 = a −1 (β) for each β, with a −1 (3. Of course we still allow for a residual dependence on a 2p2 . In order to deal with these cutoff effects, following Ref. [18], we use two methods. Method M1 consists in fitting Z RI ′ −impr ij (µ 2 0 ; a 2p2 ) to the linear ansatz   , are finally used to evaluate via the NLO running matrix formula of Ref. [15], the quantities Z MS ij (M 1) and Z RGI ij (M 1). Therefore, the MS scheme we use here is the one defined by Buras et al. in Ref. [15]. This definition of the MS scheme, which has become standard, differs from the one of Ref. [52] proposed by Ciuchini et al. in the treatment of the four-fermion evanescent operators appearing in the calculation of the two-loop anomalous dimensions.
In Fig. 10 the simultaneous best linear fits in a 2p2 at our four β's of Z ij are shown. We recall that both in the analysis and in the figures of this Appendix, only data points corresponding to the momentap satisfying the constraint (B.16) are used and shown.
The idea of the M2 method is instead to separately average at each β the values of Z RI ′ −impr ij (µ 2 0 ; a 2p2 ) over a narrow interval of momenta (ideally just one point), which has to be kept fixed in physical units for all β's. We have chosen this interval to bep 2 ∈ [ 8.0, 9.5 ] GeV 2 . In this way, at the price of giving up the reduction of cutoff effects implied by the M1 method, no assumptions are introduced in the RC analysis about the detailed form of lattice artifacts and/or the adequacy of NLO anomalous dimensions to describe the RC-evolution at scales belowp 2 ∼ 9 GeV 2 .
The a 2 -scaling of renormalized quantities (in this work operator matrix elements) constructed using RCs determined with the M2 method will of course be in general different from the one of their M1 method counterparts, but the continuum limit results for these quantities, if attainable from both methods with controlled errors, should be consistent with each other (see e.g. Appendix E).

C Renormalization Constant results
In Tables 7 and 8 we collect values of Z P and Z S calculated in the RI-MOM scheme. Results are obtained with methods M1 and M2 [18] at each value of the gauge coupling in MS and RI-MOM at 2 GeV. We have used the three-loop conversion formula from RI-MOM to MS [56].

E Results for R i and B i
In this appendix we present in detail our results in the MS scheme of Ref. [15] at 2 GeV for the quantities R i and B i (c.f. Eqs. (4.19, 4.20) and Eq. (4.18) respectively). We also give the R i results computed in the indirect way of Eq. (5.6). In Table 21 we gather results obtained employing M1-type RCs and using ChPT (NLO) fit formula, polynomial and linear fit functions (see Eqs. (5.2-5.3), and n = 2 and n = 1 of Eq. (5.1) respectively). In Table 22 we show the respective results when using M2-type RCs. Instead of using the definition of Eq. (4.19), we have employed a slightly different but equivalent one which reads where indices M1 and M2 refer to the use of the respective type of renormalisation constants and we define G . We find that the quantity defined in Eq. (E.1) has smaller O(a 2 ) effects.
In Figs. 11 and 12 we show the combined fit for the ratios,R ′ i and bag parameters B i (i = 2, . . . , 5) against the light quark mass when M2-type RCs are used.
We remark that a good agreement between the continuum limit results for the bag parameters B i and the matrix elements ratios R i obtained using M2-type RCs and their counterparts based on M1-type RCs, as we observe for i = 1, 2, 3, 4 (for i = 1 see also Ref. [8]), provides a valuable check of the smallness of residual systematic errors in the evaluation of RI-MOM RCs with the M1-method. In particular possible systematic errors stemming in the M1-method from the inadequacy at non-high momenta (p 2 ) of the perturbative operator anomalous dimensions used in the analysis or from the removal of the leading cutoff effects via a linear fit inp 2 are strongly reduced or absent when using the M2-method for RCs. This is so because in this latter approach (see Ref. [18] and Appendix B) the RCs are extracted from Landau gauge correlators at a rather highp 2 -value (fixed to ∼ 9 GeV 2 for all β's) but comes at the price of generically larger lattice artifacts on the RCs, which we partly suppress by removing the perturbatively known O(a 2 g 2 ) contributions. For the case of O i with i = 1, 2, 3, 4 the resulting cutoff effects on R i and B i (see Figs 11 and 12) appear to be under control and the continuum extrapolation is reliable. On the contrary, the case of B 5 and R 5 (see panel d) of the figures above) is a typical one where too large cutoff effects affecting the M2-type RCs make unreliable the results appearing in the i = 5-lines of Table 22.    Table 22: R i (direct computation through Eq. (E.1) and indirect computation through Eq. (5.6)) and B i results using M2-type RCs for three kinds of fit function, namely a ChPT (NLO) fit, a polynomial and a linear fit with respect to the light quark mass. For i = 2, 3 the ChPT (NLO) fit formula for R i coincides with the linear one (we refer to results of the 3rd column). The results in the lines corresponding to i = 5 here are not reliable, due to very large cutoff effects resulting in this case from the use of M2-type RCs (see text).
Finally we give our continuum results for B i and R i in the MS scheme of Buras et al., defined in Ref. [15], and the RI-MOM scheme at 3 GeV, see Tables 23 and 24 1 -14.6(5) 5.0(3) 25.6(9) 6.2(4)   Figure 12: Solid lines in panels (a) to (d) show the behaviour vs. the renormalized light quark mass of the combined chiral and continuum fits (according to the ChPT fit formula (5.2)) for the B i parameters with i = 2, . . . , 5 respectively, renormalized in the MS scheme of Ref. [15] at 2 GeV with the M2-type RCs. The full black line is the continuum limit curve (5.1). The dashed black line represents the continuum limit curve in the case of the linear fit ansatz. Black open circles and triangles stand for the results at the physical point corresponding to the ChPT fit and linear fit ansatz, respectively.