{'fraction_non_alphanumeric': 0.04484465297247811, 'fraction_numerical': 0.022586828403044, 'mean_word_length': 4.058538350217077, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
K^0-\bar K^0 mixing beyond the standard model and CP-violating electroweak penguins in quenched QCD with exact chiral symmetry

We present results for the \Delta S=2 matrix elements which are required to study neutral kaon mixing in the standard model (SM) and beyond (BSM). We also provide leading chiral order results for the matrix elements of the electroweak penguin operators which give the dominant \Delta I=3/2 contribution to direct CP violation in K->\pi\pi decays. Our calculations were performed with Neuberger fermions on two sets of quenched Wilson gauge configurations at inverse lattice spacings of approximately 2.2 GeV and 1.5 GeV. All renormalizations were implemented non-perturbatively in the RI/MOM scheme, where we accounted for sub-leading operator product expansion corrections and discretization errors. We find ratios of non-SM to SM matrix elements which are roughly twice as large as in the only other dedicated lattice study of these amplitudes. On the other hand, our results for the electroweak penguin matrix elements are in good agreement with two recent domain-wall fermion calculations. As a by-product of our study, we determine the strange quark mass. Our main results are summarized and discussed in Sec. VII. Within our statistics, we find no evidence for scaling violations.


I. INTRODUCTION
Flavor changing neutral currents (FCNC), and amongst them neutral kaon mixing, have been an extremely rich source of information about electroweak interactions and have been highly instrumental in designing the current standard model (SM) of particle physics. Like all FCNC in the SM, ∆S = 2 processes only appear at one loop and are highly suppressed. Therefore, they are also a sensitive probe for new physics and provide a rich supply of constraints for models of physics beyond the standard model (BSM), such as supersymmetric models, left-right symmetric models, multi-Higgs models, etc.
In the SM, neutral kaon mixing is dominated by W -exchange box diagrams whose short-distance contributions can be described by an effective Hamiltonian in which the W -bosons and up-type quarks have been integrated out. Because of the left-handed nature of the coupling of W -bosons to quarks, only one four-quark operator arises: where a and b are color indices.
On the other hand, in models of new physics, where heavy right-handed particles can also experience flavor changing interactions, the effective ∆S = 2 Hamiltonian is generically of the form: where, in addition to the SM operator of Eq. (1), there are the four operators, which we give here in the supersymmetric basis of [1]: as well as the parity transformed operators [1]: In Eq. (2), the Wilson coefficients C i andC i incorporate the short-distance physics associated with the new heavy degrees of freedom which may participate in extensions of the SM, as well as the SM contribution in the case of C 1 . These coefficients can be computed perturbatively. Moreover, they depend on the flavor mixing parameters of the new physics model considered. These parameters can therefore be constrained by experimental measurements of observables sensitive to ∆S = 2 processes. The relevant observables include the K L -K S mass difference, ∆M K , and the quantity ǫ, which parameterizes indirect CP violation in K → ππ decays. The short-distance contributions to these quantities are given by (cf. e.g. [2]): where the subscript I denotes the isospin of the two pion state. They are measured experimentally to a precision better than 1% (cf. [3]).
In addition to the short-distance contribution encoded in Eq. (5), ∆M K also receives "long-distance" contributions in the SM [4]. These arise when the quarks which participate in the box diagram involve the up quark. They can either be "very long-distance", involving the exchange of pseudo-Goldstone boson pairs between two ∆S = 1 vertices. Or they can be "slightly-less long-distance", involving the exchange of non-Goldstone, u, d, s hadrons. The former may be large, but their determination has very large uncertainties [5]. The latter have recently been studied in a model of large-N c QCD and are found to be of order 10% [6]. Such contributions are not a problem in the calculation of ǫ, since this quantity is dominated by the top and charm quark contributions for which the short distance approximation of Eq. (6) is very good. In general, assumptions have to be made about these long distance contributions, as well as the values of the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements in the presence of new physics. Nevertheless, under reasonable hypotheses, interesting constraints on extensions of the SM can be found, as shown for SUSY extensions, for instance, in [7].
Regardless of the hypotheses made, to validate or constrain models of new physics with measurements of ∆M K and ǫ, one clearly has to be able to quantify the non-perturbative, strong-interaction physics present in the matrix elements of the operators O i andÕ i between K 0 andK 0 states. This is where a non-perturbative technique such as lattice QCD enters. Now, while the operators in Eqs. (1), (3) and (4) have both parity even and odd parts, it is clear that only the parity even parts contribute to the matrix elements relevant for K 0 -K 0 mixing, K 0 |O i |K 0 and K 0 |Õ i |K 0 . Moreover, since the strong interaction preserves parity, we have K 0 |Õ i |K 0 = K 0 |O i |K 0 , i = 1, 2, 3. Thus, in the remainder, we need only consider the operators of Eqs. (1) and (3).
The transformation properties of the operators O i under the chiral group SU (3) L × SU (3) R are interesting. It is commonly known that O 1 belongs to the (27,1) representation of this group, which implies that its matrix element between pseudoscalar states vanishes like the square of the mass of this state in the chiral limit, as can be seen from its vacuum saturation value given below in Eq. (7). It is straightforward to work out that O 2 and O 3 belong to the (6,6) representation, while O 4 and O 5 belong to the (8,8). This means that the matrix elements of these operators between pseudoscalar states go to a constant in the chiral limit, as indicated once again by their vacuum saturation values given below in Eq. (8). Thus, the matrix elements of the non-SM operators are enhanced compared to the chirally suppressed SM one by the large factor M 2 K /(m s + m d ) 2 ≈ 25, as obtained for current values of (m s + m d ) at 2 GeV in the MS scheme. Clearly, a reliable determination of these enhanced matrix elements is needed.
The first objective of our work is to calculate the matrix elements K 0 |O i |K 0 , i = 1, · · · , 5, using lattice QCD. The calculations are performed using Neuberger fermions [8][9][10][11] on quenched gluon backgrounds generated with the standard Wilson gauge action. The gauge configurations were generated at two values of the bare coupling β = 2N c /g 0 = 6.0 and 5.85 on 18 3 × 64 and 14 3 × 48 lattices, respectively. With the lattice spacings determined in Sec. III, i.e. a −1 = 2.17 (9) GeV at β = 6.0 and a −1 = 1.49(4) GeV at β = 5.85, these lattices have sides of size 1.6 and 1.9 fm. This work extends to new matrix elements, larger lattices and to another value of the lattice spacing-allowing for an estimate of discretization errors-our earlier investigation of the SM matrix element K 0 |O 1 |K 0 , performed on a 16 3 × 32 lattice at β = 6.0 [12]. The great advantage of using Ginsparg-Wilson fermions [13] such as those proposed by Neuberger, is that unlike more traditional fermion formulations (e.g. staggered or Wilson fermions), they possess an exact, continuum-like, chiral flavor symmetry at finite lattice spacing [14]. 1 In particular, this means that the operator mixing pattern is the same as in the continuum-a non-negligible simplification for the operators considered here-and that discretization errors appear only at sub-leading a 2 order, which at currently used lattice spacings is certainly an advantage. Moreover, to avoid potentially large perturbative corrections, we perform all renormalizations non-perturbatively in the RI/MOM schemeà la [16]. Preliminary results of our calculations at β = 6.0 were presented in [17,18].
Since the matrix elements are quantities of dimension four, they suffer strongly from the uncertainty in the lattice spacing, which is of order 10% in quenched calculations [19]. It is therefore advantageous to give dimensionless measures of these matrix elements, by normalizing them with quantities which are similar enough that some of the systematic errors are expected to cancel in the ratio. For the SM contribution, determined by the operator O 1 , one usually computes the B-parameter B 1 or B K , which measures the deviation of the matrix element from its vacuum saturation value. We have where µ is the renormalization scale and M K and F K are, respectively, the mass and decay constant of the kaon, which are measured precisely experimentally. Our normalization of the decay constant is such that F K = 113 MeV and we take M K = 495 MeV. One can also define B-parameters for the matrix elements of O 2,···,5 , much in the same way. Neglecting sub-leading chiral terms, we have: for i = 2, · · · , 5 with our phase conventions, with N i = { 5 3 , − 1 3 , −2, − 2 3 } and where we have used the partial conservation of the axial current,sγ µ γ 5 d, to rewrite the matrix elements of the pseudoscalar densitysγ 5 d in terms of the s and d quark masses in the second equality. However, unlike B K which gives a dimensionless measure of the matrix element of O 1 in terms of well measured quantities, the relation of B 2,···,5 to the corresponding matrix elements involves matrix elements ofsγ 5 d, or equivalently light-quark masses, which are dimensionful quantities not directly accessible to experiment and whose determination requires a non-perturbative calculation. Thus, in providing these B-parameters, we are essentially shuffling some of the systematic errors into the determination of the normalization factor or worse, compounding them if the correlations between the parameters and the normalization factors are not taken into account. This point was first emphasized in [20], where alternate suggestions for normalizing the matrix elements were given, all involving the mass of the K * . Since the systematic errors for vector mesons may be quite different from those on pseudoscalars, we propose another dimensionless measure of the matrix elements, defined through the ratios: for i = 2, · · · , 5, where M and F are the mass and "decay constant" of the lattice kaon which we denote by P 0 to indicated that the mass of the strange and down quarks that compose it can differ from their physical values. The factor of M 2 is present to keep the ratios from diverging in the chiral limit. Though we are not interested in this limit here, the strong mass dependence induced by this divergence could complicate the necessary interpolation to the mass of the kaon. It is worth noting that the ratios R BSM i (µ, M 2 K ) measure directly the ratio of BSM to SM matrix elements and, as such, can be used in expressions for ∆M K and ǫ beyond the SM, in which the SM contribution is factored out.
Another dimensionless measure of the matrix elements which involves only quantities pertaining to pseudoscalar mesons are the ratios for i = 1, · · · , 5. From these ratios, the physical matrix elements can be obtained by considering 2F 4 K G i (µ, M 2 K ). Moreover, in the chiral limit, the G i for i = 2, · · · , 5, once multiplied by the correct Wilson coefficients, can be thought of as contributions to the couplings of operators in an effective, weak chiral Lagrangian. For simplicity, in the following we will abusively call the G i couplings.
As already mentioned, O 4 and O 5 belong to the (8,8) representation of SU (3) L × SU (3) R . This is the same chiral representation as the ∆I = 3/2 parts of the electroweak penguin operators Q 7 and Q 8 , which give the dominant ∆I = 3/2 contribution to ǫ ′ . In the SU (3)-flavor limit, we have: where the relative phase between the two matrix elements depends on the phase conventions chosen for pion and kaon states 2 and where we define: In this same limit we also have where B Moreover, at leading order in the chiral expansion, these K + → π + matrix elements are related to the K 0 → (ππ) I=2 matrix elements relevant for the calculation of ǫ ′ in the following way: where F is the pseudoscalar decay constant which, at this order, can be taken to be F π , F K or any other value in the chiral regime and where the constant of proportionality depends on phase and normalization conventions. Eq. (15) becomes exact, of course, in the chiral limit. Thus, the second objective of our work is to determine the matrix elements π + |Q 3/2 7,8 |K + /F , in the chiral limit. Ideally then, we would compute G 5,4 of Eq. (10) on the lattice for mesons P in the chiral regime and extrapolate the result to the chiral limit to obtain the electroweak penguin couplings G 3/2,χ 7,8 of the corresponding weak, effective chiral Lagrangian, or equivalently, the matrix element (ππ) I=2 |Q 7,8 |K 0 in the chiral limit. However, two problems stand in our way. The first is that the mesons in our simulation, except for perhaps the two to three lightest ones, are beyond the chiral regime. The second is that the chiral behavior of G 5,4 in the quenched and the N f = 3 unquenched theory differ, indicating that we cannot expect quenched results for G 5,4 in the chiral limit to have much to do with their physical values.
Due to computational limitations, we cannot do much about the first problem at present, but we can try to address the second. Indeed, as discussed in Sec. VI and first argued in [18], the matrix elements π + |Q 3/2 7,8 |K + can be normalized in such a way that the resulting quantities have the same NLO chiral logarithm in the quenched and N f = 3 theories. In the SU (3)-flavor limit, we can consider where F is the decay constant common to the pions and kaons. This similarity in the chiral expansions allows for the possibility that the quenched and unquenched results for these quantities in the chiral limit may not be too different. 2 In our conventions, the matrix elements π + |Q 3/2 7,8 |K + differ by a sign from the ones given in [21].
Because of its importance in constraining the CP violating phase of the Cabibbo-Kobayashi-Maskawa (CKM) matrix with the measurement of ǫ (cf. e.g. [22,23]), B K has received much attention, both analytically and in lattice QCD. Analytically, it has been evaluated with two-point [24,25] and three-point [26][27][28] QCD sum-rules, as well as in models of large N c QCD [29][30][31]. The most recent calculations were performed in the chiral limit, except for that of [29], where it was shown that chiral corrections represent approximatively 50% of the final number. On the lattice B K also has a long history, which has been reviewed recently in [32][33][34]. In particular, there are at present two quenched benchmark calculations [35,36] and, very recently, the first N f = 2 + 1 unquenched calculation with staggered fermions, at a single value of the lattice spacing [37]. The point of our calculation is not to establish a new quenched benchmark for B K , but rather to quantify K 0 -K 0 mixing beyond the SM in a consistent framework. In fact, our dimensionless ratios R BSM i of Eq. (9), in which a number of systematic errors are expected to cancel, may be used in the future, in conjunction with high precision unquenched determinations of B K , to yield improved calculations of K 0 -K 0 mixing beyond the SM.
Unlike the situation that we have described in the preceding paragraph for B K , there is only one other lattice calculation of the full set of ∆S = 2 matrix elements relevant for K 0 -K 0 mixing [20]. It was performed with tree-level O(a)-improved Wilson fermions, which are not chirally symmetric and suffer for potentially large O(α s a) errors. In fact we find results for the BSM matrix elements which are substantially larger than those of [20], implying stronger constraints on BSM models than previously found.
The situation regarding electroweak penguin matrix elements is somewhat intermediate. There are a number of analytical calculations [21,[38][39][40] as well as a handful of recent quenched lattice calculations [41][42][43][44][45] and a preliminary N f = 2 unquenched calculation [46]. There are also older calculations of the corresponding B-parameters [47,48] as well as exploratory calculations which attempt to go beyond the chiral limit [45,46], including one [49] which is based on the finite-volume formalism of [50]. Yet there is no consensus on the value of these matrix elements even in the chiral limit. Thus, our results, obtained in a formulation of lattice QCD with full chiral symmetry, full O(a)-improvement, at two values of the lattice spacing and using non-perturbative renormalization, are useful for clarifying the current situation.
Our paper is organized as follows. In Sec. II we briefly summarize the parameters of our simulation. In Sec. III we describe our analysis of two-point functions relevant for the determination of the matrix elements of interest. This includes a description of our extraction of the strange quark mass. In Sec. IV we use three-point functions to determine the bare B-parameters and other quantities related to the four-quark operator matrix elements, as functions of our bare quark masses. Sec. V is dedicated to the non-perturbative calculation of the renormalization constants required to subtract the logarithmic divergences present in the bare results. In Sec. VI we discuss the quark-mass dependence of the results renormalized conventionally at 2 GeV 3 in the RI/MOM scheme and discuss the necessary interpolations and extrapolations to the physical quark masses. As in the previous sections, the results obtained at the two values of the lattice spacing are discussed in parallel, and a direct comparison of these results is presented. In Sec. VII we combine the results obtained at the two values of the lattice spacing for the matrix elements, matrix element ratios and quark masses interpolated or extrapolated to their physical mass point. We give our final physical results at 2 GeV in both the RI/MOM and MS-NDR schemes, including estimates of the dominant systematic errors within the quenched approximation. We also compare our results to those obtained both on and off the lattice. The reader interested only in our final results can skip directly to this section. Finally, in Sec. VIII, we conclude with a few closing remarks.
The paper further contains a number of appendices. Appendices A-D contain tables which are referred to in earlier sections, but which we place in appendices to improve the readability of the paper. For completeness, we provide in Appendix E the non-perturbative renormalization constants in the MS-NDR scheme as well as the quark-mass dependence of the four-quark matrix element results renormalized in the MS-NDR scheme at 2 GeV in Appendix F. We also give, in Appendix G, results in the MS-NDR scheme at the individual lattice spacings for the B-parameters and other matrix element ratios at the physical point. These results are used in Sec. VII to obtain our final physical results in that scheme.

II. NUMERICAL IMPLEMENTATION
We describe gluons with the Wilson gauge action and quarks with the overlap action where m q is the bare quark mass and D[0] is the Neuberger-Dirac operator [9] In Eq. (18), X = D W − ρ/a is the Wilson-Dirac operator with a negative mass parameter ρ which separates physical from doubler modes. We generated 100 quenched gauge configurations each at a bare value of the gauge coupling, β = 6.0 and β = 5.85, corresponding to inverse lattice spacings of a −1 = 2.12 GeV, respectively a −1 = 1.61 GeV, as obtained from the Sommer scale defined by r 2 0 F (r 0 ) = 1.65 and r 0 = 0.5 fm [51,52]. 4 On each gauge configuration, we invert the Neuberger-Dirac operator for a number of values of the quark mass ranging from roughly m s /3 to 1.5m s . At β = 6.0, the bare masses in lattice units are am q = 0.030, 0.040, 0.060, 0.080, 0.100 and at β = 5.85, they are am q = 0.030, 0.040, 0.053, 0.080, 0.106, 0.132. Details on our implementation of the Neuberger-Dirac operator can be found in [53]. All statistical errors are estimated with a single-elimination jackknife, which is propagated throughout our analysis.
In our simulations, we choose the negative mass parameter ρ to vary with the lattice spacing a, namely ρ = 1.4 at β = 6.0 and ρ = 1.6 at β = 5.85, to optimize the locality of the Neuberger-Dirac operator [54]. At fixed lattice spacing, a change in ρ corresponds to an O(a 2 ) redefinition of the action and one may worry about the effect on the continuum limit of our procedure. However, as long as ρ, viewed as a function of the lattice spacing a, remains in the single particle sector and has a smooth limit in the range 0 < ρ < 2 for a → 0, changing ρ with β simply represents a modification of the approach to the continuum limit. Moreover, for a function ρ which approaches its continuum limit value with corrections of O(a 2 ), this modification only appears at O(a 4 ). We require that our two values of ρ lie on such a curve.

III. TWO-POINT FUNCTIONS, SPECTRAL QUANTITIES AND SCALE SETTING
To obtain the vacuum saturation values of the matrix elements of interest, we compute the following two-point functions: and where P ≡dγ 5ŝ is the pseudoscalar density and A 0 is the time-component of the axial current A µ ≡dγ µ γ 5ŝ . Herê s = (1 − aD[0]/2ρ)s. The introduction of the "hatted" quark fields guarantees that the bilinear and four-quark operator matrix elements have the proper chiral properties and therefore renormalize as in the continuum and are free of O(a) discretization errors. MoreoverP andĀ 0 are obtained from P and A 0 by interchanging the flavors d ↔ s. We work in the SU (3)-flavor limit where m s = m d = m q and we allow the common bare quark mass am q to take on the values described in the preceding section. While the correlation function C P A0 (x 0 ) is only used to normalize the three-point function involving operator O 1 as described below, we fit the time-symmetrized correlation function C asymptotic behavior has been reached while the final time is chosen to match the range used in three-point function fits, as described in Sec. IV. Because of the chiral symmetry of Neuberger fermions, the axial Ward identity can be used to derive the pseudoscalar decay constant F from this matrix element, without a finite renormalization. Results for aM and aF are summarized in Table XI for β = 6.0 and Table XII for β = 5.85.
Because Neuberger fermions satisfy an index theorem [14,55], one may worry about the presence of finite-volume, quenched, zero mode contamination in the C (2) P P (x 0 ) correlator. However, we showed in [53] that these contributions are not visible within our statistical accuracy at the volumes and quark masses that we consider. In [53], we compared directly results of fits to C (2) P P (x 0 ) and C (2) (P +S)(P −S) (x 0 ), where S is the scalar density, and found no appreciable difference. The two correlators are both dominated by the pseudo-Goldstone boson state for asymptotic times x 0 and C (2) (P +S)(P −S) (x 0 ) is free from zero-mode contributions by chirality. For completeness, this comparison is repeated here. Results obtained for the time-symmetrized correlation function C (2) (P +S)(P −S) (x 0 ), for the same time fitting ranges as above, are also given in Tabs. XI-XII.
As Tabs. XI-XII show, C P P (x 0 ) yields results for aM and aF which are compatible with those given by C (2) (P +S)(P −S) (x 0 ). The compatibility is also displayed graphically in Figs. 1-2, where the results for aF from the two correlation functions are plotted against the corresponding values M 2 /(4πF ) 2 . We conclude that the effects of finite-volume zero-mode contributions are not statistically significant at the volumes and quark masses that we consider, at least for two-point functions at asymptotic times. For simplicity, we will consider only C (2) P P (x 0 ) in the remainder of this paper.
To determine the lattice spacing a, we study aF as a function of M 2 /(4πF ) 2 and interpolate aF quadratically to the point M 2 /(4πF ) 2 = M 2 K /(4πF K ) 2 , with M K and F K given in the Introduction. The interpolations are shown in Fig. 1 for β = 6.0 and in Fig. 2 for β = 5.85, where the result corresponding to am q = 0.030 is omitted because M P L ≃ 3.5 for this point, which makes it more susceptible to finite volume effects. With the value of the decay constant in lattice units obtained at the kaon point, we set a −1 = F K /(aF ) and obtain a −1 = 2.17(9) GeV at β = 6.0 and a −1 = 1.49(4) GeV at β = 5.85. These values for the lattice spacing are compatible with those derived from the Sommer scale in Sec. II, namely a −1 = 2.12 GeV at β = 6.0 and a −1 = 1.61 GeV at β = 5.85, given the intrinsic O(10%) model-dependence associated with r 0 . We choose to set the scale and fix the strange quark mass with M K and F K because these quantities are closely related to the matrix elements of interest. This procedure may therefore absorb some of the quenching artefacts associated with the physics that we are studying. In any event, we will be considering mainly dimensionless quantities such as B-parameters and matrix-element ratios. The value of the lattice spacing will then only enter logarithmically through the renormalization scale dependence of the matrix elements and in fixing the kaon point in the mild interpolations that we will have to perform.
As a side product of our analysis, we determine the value of the strange quark mass. The renormalized quark mass is obtained from the bare quark mass in lattice units as m RI (2 GeV) = (am) * a −1 /Z RI S (2 GeV), where the lattice spacing is the value determined above, and Z RI S (2 GeV) is the renormalization constant of the scalar density determined non-perturbatively in Sec. V. We assume that the meson mass depends only on the sum of the masses of the quarks which compose it. In Fig. 3, we show a quadratic interpolation to the kaon point of the RI/MOM sum of quark masses at 2 GeV, (m q1 + m q2 ) RI (2 GeV), as a function of M 2 /(4πF ) 2 at β = 6.0. Also shown in this plot are the results at β = 5.85, where the interpolation, performed without the point corresponding to am q = 0.030 for the reasons discussed above, has been omitted for clarity. Only statistical errors are shown, and already the agreement is very good. The interpolation to the kaon point gives (m s +m) RI (2 GeV), wherem is the average up and down quark mass. Our results for this quantity, obtained at the two values of the lattice spacing are summarized in Table I. There is no evidence for scaling violations. Results in the MS scheme are given in Appendix G. The final result will be given in Sec. VII.
shown as a vertical dotted line. The dashed curve is an extension of the fit curve outside the fit range. Also shown are the results at β = 5.85, but the interpolation is omitted for clarity. Only statistical errors are shown.

IV. THREE-POINT FUNCTIONS, BARE B-PARAMETERS AND FOUR-QUARK OPERATOR MATRIX ELEMENTS
To obtain the bare matrix elements of the ∆S = 2 operators, O i , between K 0 andK 0 states, and the corresponding B-parameters, we construct the following ratios of three-to two-point functions: for i = 2, · · · , 5, with the N i 's given after Eq. (8). As indicated, these ratios become constant for asymptotically large initial (y 0 ) and final (x 0 ) times. These constants are the bare B-parameters, B i , i = 1, · · · , 5, whose renormalized analogs were defined in Eq. (8). In Eqs. (21)- (22), the subscripts on B i J f Ji , i = 1, · · · , 5, label the sink J f (x) and the source J i (y), which we defined below Eq. (20) for the case of P (x) and A 0 (x). In addition, the bare operators O i in Eqs. We fit the B i P P ratios of Eqs. (21)- (22) to a constant in the symmetric time intervals given by: 12 ≤ x 0 /a ≤ 19 and 45 ≤ y 0 /a ≤ 52 for i=1, · · · , 5, at β = 6.0; 10 ≤ x 0 /a ≤ 12 and 36 ≤ y 0 /a ≤ 38 for i=1 and 10 ≤ x 0 /a ≤ 14 and 34 ≤ y 0 /a ≤ 38 for i=2, · · · , 5, at β = 5.85. The points nearest the origin are are chosen so that the ratios are dominated by the ground state contribution. The other extremities of the ranges are picked so that time reversed contributions are less than 1%, even for the lightest meson (next to lightest at β = 5.85). The results for the bare B-parameters are given in Table XIII at β = 6.0 and Table XIV at β = 5.85, as functions of bare quark mass.
Though no evidence was found in the previous section for finite-volume, zero-mode contamination in two-point functions, one could worry about their effect on three-point functions. In order to verify whether these modes might be a problem, we considered an alternate ratio for determining B 1 which, unlike B 1 P P (x 0 , y 0 ) of Eq. (21), is entirely free of zero-mode contamination because of its chiral structure. This ratio is the one that we advocated in our first determination of B K with Neuberger fermions [12], performed on a lattice whose four-volume is roughly three times smaller than the one considered here: where L 0 is the time-component of the left-handed weak current, L µ =dγ µ (1 − γ 5 )ŝ. As found in [12], this ratio reaches asymptotic behavior at early times. Fits of this ratio in symmetric time ranges with x 0 /a ≥ 8 and y 0 /a ≤ 56 at β = 6.0 and with x 0 /a ≥ 6 and y 0 /a ≤ 42 at β = 5.85 gives bare values of B 1 which are entirely compatible with those obtained from B 1 P P (x 0 , y 0 ). We take this as evidence that the three-point functions which appear in Eqs. (21)- (22) are not significantly affected by zero-mode contributions for the masses that we consider in our relatively large volumes.
In computing matrix elements of the operators O 2 and O 3 , it is also possible to find sources which entirely eliminate finite-volume, zero-mode contributions. In the definition of B 2,3 P P (x 0 , y 0 ) in Eq. (22), one must replace the pseudoscalar source and sink in the numerator by R ≡ S + P , where S is the scalar density, S ≡dŝ. This eliminates zero-mode contributions in the three-point function by chirality. In the denominator, the two point function P (x)P (0) must be replaced by R(x)L(0) , and P (0)P (y) by L (0)R(y) , where L ≡ S − P . Both these two-point functions are free from zero modes. The resulting ratios, B 2,3 RR (x 0 , y 0 ) are thus also free from zero-mode contributions. As one can easily convince oneself, such a manipulation is not possible in correlation functions relevant for computing the matrix elements of the operators O 4 and O 5 . There, zero-mode contributions cannot be eliminated by simply using chirality. Nevertheless, in [18], we have checked that sources with different chiral properties, and thus different zero-mode contributions, give identical results. In any case, we do not pursue this issue further here and assume that the absence of zero-mode contributions observed in pseudoscalar two-point functions and in the three-point function of Eq. (21) is evidence that these contributions are not a problem. And we take as our final results for the B-parameters those given by the ratios of Eqs. (21)-(22), because the use of identical sources for the different operator matrix elements yields cancellations of some statistical fluctuations in the ratio of these matrix elements.
From these B-parameters and our determinations of M and F from the two-point function C (2) P P , we then derive a number of other quantities of interest. The building blocks for these quantities are the bare matrix elements of our five operators O i , i = 1, · · · , 5, which we obtain as follows: In Eq. (24), Z A is the normalization constant of the local axial current, A µ , defined after Eq. (20). It is there because the vacuum saturation value of the matrix element of O 1 , used in the definition of B 1 , is obtained in terms of the matrix element of A 0 , while in constructing the matrix element in Eq. (24), we have used the partially conserved axial current, through the axial Ward identity. Z 2 A could be omitted here and reintroduced when we renormalize our results. However, for our simulation parameters it leads to a large correction (see Table II). We have thus chosen to introduce it, because it allows for direct comparison of quantities containing P 0 |O 1 |P 0 of Eq. (24) to those that are obtained using a vacuum saturation calculated with matrix elements of the local axial current A 0 .
In fact, we have performed such a comparison using 0|A 0 |P 0 obtained from both C P A0 and C A0P , together with C (2) P P . The bare matrix elements determined using the three different estimates of its vacuum saturation value are entirely compatible. However, the matrix element's statistical error does depend quite strongly on the procedure used to obtain it from B 1 . The matrix element obtained using a vacuum saturation value from C (2) P A0 is the one with the smallest statistical error. It is followed by the one obtained which makes use of C (2) A0P and finally by the one of Eq. (24), which makes use of only C (2) P P . Nonetheless, the best measure of the physical matrix element K 0 |O 1 |K 0 is given by 16 3 M 2 K F 2 K B 1 , which does not necessitate reconstructing P 0 |O 1 |P 0 for our lattice mesons, P 0 . Moreover, obtaining P 0 |O 1 |P 0 from B 1 and a fit to C (2) P A0 gives results for the BSM ratios of Eq. (9) which have significantly larger statistical errors than those obtained using Eq. (24). The same is true of the matrix element obtained from B 1 and a fit to C (2) A0P . In fact, both these determinations of BSM ratios have very comparable errors. We thus take P 0 |O 1 |P 0 as obtained through Eq. (24) and avoid discussing fits of the two-point functions C (2) P A0 and C (2) A0P which have already been detailed in [53].
In addition, as discussed in the Introduction, we have chosen not to give the values of the matrix elements themselves, which have mass dimension four and are therefore very sensitive to the intrinsic uncertainty associated with the determination of the lattice spacing in quenched calculations. Rather we choose to present two dimensionless measures  (18) of these matrix elements, namely the BSM ratios of Eq. (9) and the couplings of Eq. (10), the latter enabling a determination of the matrix element itself for a given meson mass, using the corresponding value of the leptonic decay constant.
Thus, from the matrix elements obtained using Eqs. (24)- (25) and the values of M and F from fits to C P P (Tabs. XI-XII), we straightforwardly construct the bare equivalents of the BSM ratios defined in Eq. (9). We summarize our results for these ratios as a function of bare quark mass at β = 6.0 in Table XV and at β = 5.85 in Table XVI.
We also construct the bare equivalents of the couplings G i of Eq. (10), from the matrix elements as given by Eqs. (24)- (25), and the pseudoscalar decay constant obtained from C (2) P P and the axial Ward identity and given, again, in Tabs. XI-XII. Our results for these contributions to the bare couplings as a function of bare quark mass are summarized in Table XVII   To avoid potentially large perturbative errors, we compute all renormalization constants non-perturbatively in the RI/MOM schemeà la [16]. Thus, we fix gluon configurations to Landau gauge and numerically compute the relevant amputated forward quark Green functions with legs of four-momentum p, Λ Oi (m q , p, g 0 ), where m q is the mass of the quark. The renormalization constants are then determined by requiring that the corresponding renormalized vertex functions have their tree-level values. We therefore define the ratios: . Because of the chiral symmetry of Neuberger fermions, the parity even and odd parts of the operators renormalize in the same way. We therefore consider the parity odd parts, as CPS symmetry then guarantees that off-diagonal mixing terms vanish exactly. Chiral symmetry also guarantees that the product Z A × Tr {Λ V (m q , p)P V } is the quark wavefunction renormalization squared, Z q . Since we can obtain Z A very cleanly from the axial Ward identity [53], this method allows us to circumvent a more direct calculation of Z q , a procedure which induces large discretization errors. We choose Λ V instead of Λ A because the former has a slightly cleaner signal. In Table II, we summarize the results for Z A obtained from the axial Ward identity.
The momentum combinations that we consider in constructing the ratios of Eq. (26) are the following. We allow p 1,2,3 × (L/2π) to take on values from 0 to 2, and p 0 × (T /2π), values from 0 to 4. For larger momenta, we restrict p 1,2,3 × (L/2π) to take on values from 2 to 3 and p 0 × (T /2π), values from 7 to 10. Moreover, we use the lattice-fermion definition of momentum, p 2 ≡ (1/a 2 ) 3 µ=0 (sin ap µ ) 2 . These choices reduce O(4)-breaking discretization errors so that the R RI ij (m q , p, g 0 ) are to a good approximation functions of p 2 . At large p 2 (p 2 > ∼ 1.5 GeV 2 ), an operator product expansion (OPE) in 1/p 2 can be performed on the Green functions which appear in the definition of R RI ij (m q , p, g 0 ). In this expansion, the leading term, proportional to the unit operator, gives the desired renormalization constants in the RI/MOM scheme. Perturbative mass effects appear at sub-leading orders, suppressed by powers of p 2 , as do possible non-perturbative contributions. On the lattice, of course, O(4) symmetry is broken and there will be discretization errors. However, as already noted, for our choice of momenta, we expect the R RI ij (m q , p, g 0 ) to be functions of p 2 to good approximation. Therefore, the leading terms in the OPEs of these ratios at finite m q , including discretization error terms, are of the form (see also [12]): where the ellipses stand for higher-order terms in the OPE and sub-leading discretization errors. Here U RI ik (p 2 ) describes the running of the renormalization constants in the RI/MOM scheme, Z RI ij (p 2 ) = U RI ik (p 2 )Z RGI kj , where the superscript "RGI" (renormalization group invariant) indicates that the corresponding renormalization constant is the constant which translates lattice results into scale-and scheme-independent quantities. 5 The p 2 dependence in A (2) ij and B (2) ij is specified to accommodate the expected logarithmic corrections to the power behavior in p 2 . Note that B (2) ij is independent of quark mass, because this dependence only appears at higher orders in the OPE.
To determine whether power corrections and discretization errors are relevant in extracting the renormalization constants from the R RI ij (m q , p, g 0 ), we consider the "RGI" ratios: obtained from the R RI ij (m q , p, g 0 ) by dividing out the running. In the absence of discretization errors and sub-leading perturbative and OPE terms, these "RGI" ratios should be constant functions of p 2 . In Fig. 4, we plot R RGI ij (m q , p, g 0 ) for i = j = 1 as functions of p 2 for our smallest and largest quarks masses, am q = 0.03 and 0.10 at β = 6.0. We do the same in Fig. 5 for R RGI ij (m q , p, g 0 ) with i, j ∈ 2, 3, and in Fig. 6 for the mixing of operators O 4 and O 5 . Here the running is implemented at two loops [57,58], with N f = 0 and Λ QCD = 0.238 (19) GeV from [59]. A number of observations are in order: • to a very good approximation, R RGI ij (m q , p, g 0 ) are functions of p 2 ; • the R RGI ij (m q , p, g 0 ) are not generically constant functions of p 2 at values of p 2 where one might consider extracting the renormalization constants; when they are, this behavior seems to be the result of a turnover instead of genuine p 2 -independence; • the R RGI ij (m q , p, g 0 ) do in some cases exhibit a significant mass dependence, but that dependence tends to vanish when p 2 increases as predicted by the OPE; • discretization errors are clearly visible but remain moderate and there is no evidence for contributions of order higher than (ap) 2 .
Our conclusion is that power corrections and discretization errors may give significant contributions and they should be accounted for in extracting the renormalization constants, as already briefly discussed in [12]. This is even more true at β = 5.85 where discretization errors enter at lower physical momenta. So, instead of performing, in some cases delicate, chiral extrapolations of R RI ij (m q , p, g 0 ) to m q = 0 6 and taking Z RI ij (g 0 , p 2 ) = R RI ij (m q , p, g 0 ) for p 2 ≃ 4 GeV 2 as is usually done [16,56], we fit the R RI ij (m q , p, g 0 ) at finite m q to OPE forms such as the one given in Eq. (27). This procedure not only eliminates the mass effects through their p 2 dependence, but also other power-suppressed contributions as well as discretization errors. In performing these fits, we neglect possible logarithmic corrections in p 2 in the coefficients A (2,4) ij and B (2) ij . They are at least subdominant to the leading corrections and sub-subdominant compared to the perturbative behavior that we want to isolate. 7 Moreover, they are not known analytically and they would be difficult to identify in the data. We also neglect discretization errors in (am) 2 since they should be significantly smaller than either the (ap) 2 discretization errors or the 1/p 2 corrections in the range of p 2 considered.
We obtain our central values for the renormalization constants by fitting to all points p 2 such that 3 GeV 2 ≤ p 2 ≤ 14 GeV 2 at β = 6.0 and 3 GeV 2 ≤ p 2 ≤ 8.5 GeV 2 at β = 5.85. The reason for this rather large range of p 2 is that we want to make sure that all OPE and discretization terms, which make a significant contribution to R RI ij (m q , p, g 0 ) in regions where we can hope to extract reliable values of the renormalization constants, are accurately measured with the data. Moreover, down to 3 GeV 2 , the two-loop running remains under control. For R RI 11 (m q , p, g 0 ), we consider even lower momenta, down to p 2 = 1.5 GeV 2 . The reason for this is the presence of a large power correction whose coefficient can only be determined at these low values of p 2 . And because the running of O 1 is mild, the two-loop evolution of this operator remains sensible even down to these low values of the momentum. The terms kept in the fits for the various renormalization constants are summarized in Table III. We initially individually fit R RI ij (m q , p, g 0 ) obtained for our lightest or next to lightest and heaviest quark masses ((am min q , am max q ) = (0.030, 0.100) at β = 6.0 and (0.040, 0.132) at β = 5.85). We find that, despite a factor of over 3 between these two masses, the resulting Z RGI ij (g 0 ), as well as the discretization error coefficients B ij (g 0 ), are independent of quark mass within statistical errors, giving us confidence in the consistency of our OPE description. Therefore, for our final results we perform a combined fit to the R RI ij (m q , p, g 0 ) for the two masses, constraining Z RGI ij (g 0 ) and B ij (g 0 ) to be mass-independent, but allowing, of course, A ij to be different for the two masses. The resulting fits at β = 6.0 are shown in Figs. 4-6, together with the horizontal line describing the behavior of the perturbative term alone, as reconstructed from the fit. The distance between this line and the data points quantifies the impact of OPE and discretization corrections to the leading perturbative term. Similar plots are obtained at β = 5.85, though the stability of the fits is reduced by the presence of discretization errors at even lower values of p 2 . The resulting renormalization constants, Z RI ij (2 GeV), in the RI/MOM scheme at 2 GeV are summarized in  Fig. 4, but for the scalar density. Here, running is implemented at four loops, as described in the text.
for β = 6.0 and Table XXII for β = 5.85. Since the MS-NDR scheme is also often used, for completeness we give the renormalization constants in that scheme in the Appendix E in Tabs. XXXVIII-XXXIX.
To renormalize the B-parameters, we further need to renormalize the matrix elements' vacuum saturation values. For O 2,···,5 , these are given by matrix elements of the pseudoscalar density, so that the renormalization constant needed is Z P . This constant is actually difficult to obtain directly in the RI/MOM scheme of [16], because the relevant vertex receives chirally enhanced power corrections [16]. Fortunately, the chiral symmetry of Neuberger fermions guarantees that Z P = Z S , where Z S is the renormalization constant for the scalar density, which is free of Goldstone pole contaminations. Thus, we consider a ratio R RI S (m q , p, g 0 ), which is analogous to the ratios R RI ij (m q , p, g 0 ) defined in Eq. (26), with the operators O i and O j replaced by the scalar density and the squares replaced by first powers in the numerator. As above, we study the momentum dependence of this ratio and fit it to an OPE expression of the form given in Eq. (27), in the interval 2 GeV 2 ≤ p 2 ≤ 14 GeV 2 at β = 6.0 and 2 GeV 2 ≤ p 2 ≤ 8.5 GeV 2 at β = 5.85. We find that a good description of the data is obtained by retaining the discretization term proportional to (ap) 2 together with 1/p 2 and 1/p 4 OPE terms (see Table III). The data and the fit at β = 6.0 are shown in Fig. 7, together with the horizontal line describing the behavior of the perturbative term alone, as reconstructed from the fit. One difference here is that the running of the operator is implemented at four loops [60][61][62], with N f = 0 and the same value of Λ QCD as above. The result in the RI/MOM scheme at 2 GeV is given in Table XXIII. The value at the same scale in the MS scheme is given in Appendix E in Table XL. These constants are also required for determining the renormalized quark masses of Tabs. I-LV, from the corresponding bare quark masses.
Though it is straightforward to obtain the renormalization constants for the B-parameters from those for the operators in Tabs. XXI-XXII and for Z P in Table XXIII, we nevertheless give our results for these constants in the RI/MOM scheme at β = 6.0 in Table XXIV and at β = 5.85 in Table XXV, because they may be of use and because we can account for the correlations in the ratios of Z ij to Z 2 P . Results in the MS-NDR scheme are summarized in Tabs. XLI-XLII in Appendix E.
To estimate the systematic errors associated with our procedure for determining renormalization constants, we consider fits to a higher momentum range, where power corrections in 1/p 2 will play a less important rôle. We choose the range 7.5 GeV 2 ≤ p 2 ≤ 14 GeV 2 at β = 6.0 and 5 GeV 2 ≤ p 2 ≤ 8.5 GeV 2 at β = 5.85, where momenta belong to the high momentum cone described below Eq. (26). We also eliminate the power corrections in 1/p 2 to which our results are not sensitive in this higher range of momenta. The list of terms kept is summarized in Table IV and the resulting renormalization constants are given in Tabs. XXVI-XXVII in the RI/MOM scheme. The resulting values for Z S are given in Table XXIII and the results for the renormalization constants of the B-parameters are summarized in Tabs. XXVIII-XXIX.
As the renormalization constants obtained from the two fitting ranges indicate, our description of the RI/MOM  IV: Same as Table III, but for the restricted range 7.5 GeV 2 ≤ p 2 ≤ 14 GeV 2 at β = 6.0 and 5 GeV 2 ≤ p 2 ≤ 8.5 GeV 2 at β = 5.85.
am min q am max q 1/p 4 1/p 2 1/p 4 Green functions is very stable at β = 6.0, giving us confidence that we have isolated correctly the various OPE and discretization terms which contribute to these functions. However it is clear that for a high precision calculation of matrix elements, a non-perturbative renormalization procedure in which corrections in powers of 1/p 2 are absent and in which perturbation theory is only used at much higher scales would be highly desirable, especially on coarser lattices. Proposals in that direction have been made very recently in [63,64].

VI. MASS DEPENDENCE AND PHYSICAL MATRIX ELEMENTS
Having determined the renormalization constants, we can now renormalize the bare matrix elements and Bparameters for our various meson masses. Renormalizing these quantities before studying their behavior in mass is useful because it allows for comparisons of this behavior with results obtained in other simulations, performed with either different actions or different parameters. In our case, we undertake such a comparison between the results obtained on our finer and coarser lattices, at β = 6.0 and β = 5.85, respectively. Moreover, this renormalization is also important if one is to perform a combined chiral and continuum extrapolation using finite lattice spacing chiral perturbation theory.
In Tabs. XXX-XXXI we present the results obtained at β = 6.0 and β = 5.85, respectively, for the B-parameters in the RI/MOM scheme at 2 GeV, as functions of bare quark mass. The corresponding results for the renormalized BSM ratios, R BSM i , are given in Table XXXII for β = 6.0 and Table XXXIII for β = 5.85. In Table XXXIV for β = 6.0 and Table XXXV for β = 5.85, we summarize our results for the couplings G i of Eq. (10), again as a function of light-quark mass. Finally, Tabs. XXXVI-XXXVII display the results that we obtain for D 3/2 7,8 in the RI/MOM scheme at β = 6.0 and β = 5.85, respectively. For completeness, we present all of these results renormalized in the MS-NDR scheme at 2 GeV in Appendix F, in Tabs. XLVII-LIV.
In order to finalize our calculation of the B-parameters, BSM ratios R BSM i and couplings G i relevant for K 0 -K 0 mixing beyond the standard model, we have to interpolate our renormalized results to the kaon mass. To perform the interpolation we fit our data, with the exclusion of the single point for am q = 0.030 at β = 5.85 which we leave out for the reasons explained in Sec. III, to simple polynomial forms in M 2 /(4πF ) 2 . As we did for setting the scale in Sec. III, we fix the kaon point to be M 2 /(4πF ) 2 = M 2 K /(4πF K ) 2 . More complicated chiral perturbation theory fits could be tried, but their applicability in our mass range is questionable and for the simple interpolations that we have to perform, any reasonable description of the data is sufficient.
We begin with the interpolation of the B-parameters. We find that B 1 , and B 5 are consistent with linear behavior in M 2 /(4πF ) 2 while B 2 , B 3 and B 4 are well fit by parabolas in our range of masses. We show these fits for illustrative purposes in Fig. 8, but do not reproduce the fit parameters here since they are not of physical interest. Our results for the B-parameters at the kaon mass in the RI/MOM scheme at 2 GeV are summarized in Table V     We further interpolate the couplings G i . As Fig. 10 confirms, G 5 is well fit by a constant, G 4 by a straight line, while G 1 , G 2 and G 3 require quadratic terms in M 2 /(4πF ) 2 . Again, we do not reproduce the fit parameters here. We only note that in the chiral limit G 1 is consistent with zero, which is the expected behavior. The results for the couplings G i at the kaon mass in the RI/MOM scheme at 2 GeV are given in Table V.
From these results for B i , G i and R BSM i at the kaon mass in the RI/MOM scheme, it is straightforward to obtain the corresponding results in the MS-NDR scheme at 2 GeV. This is done in Appendix G and the results are presented there in Table LVI. Finally, we turn to the most delicate part of this discussion: the chiral extrapolation of the electroweak penguin matrix elements. As stated in the Introduction, the quantities D  Fig. 8, but for the "couplings" Gi, i = 1, · · · , 5. Also shown, for indicative purposes, is an extrapolation of G1 to the chiral limit, which is denoted by a dashed vertical line at M 2 /(4πF ) 2 = 0.  NLO in the quenched and N f = 3 theories, while the corresponding matrix elements or parameters G 5,4 do not. Indeed, at NLO in the N f = 3 theory [65], and [66] In the quenched theory, we have at NLO [67], and [68] π + |Q 3/2 where it is understood that the low energy constants ρ F , . . . may differ in the quenched and N f = 3 theories. Thus, for D 3/2 7,8 , we find in both cases, the following chiral behavior: where again, the low energy constants are distinct in the two theories. However, the fact that the functional behavior is the same and the fact that the coefficient of the logarithm is small suggests that D 3/2,χ 7,8 may be comparable in the two theories, especially since our lattice has been tuned to reproduce results at the kaon mass and the extrapolation should be mild. This mildness of the mass behavior is actually seen in our results for D 3/2 7,8 renormalized in the RI/MOM scheme at 2 GeV, which are plotted as functions of M 2 /(4πF ) 2 in Fig. 11.
Instead of attempting a hazardous chiral extrapolation of D 3/2 7,8 , we choose to interpolate these quantities to the kaon point. Because F K is used to fix the lattice cutoff, we take D , with G 5,4 (µ, M 2 K ) obtained above in the RI/MOM scheme at 2 GeV. We then use the fact that any reasonable chiral extrapolation (linear, quadratic, log form of Eq. (33), log form plus a quadratic term in M 2 /(4πF ) 2 , . . .) gives values which are compatible with the results at the kaon point. This is shown in Fig. 11 (16) interpolated to the kaon point and our best estimate for (ππ)I=2|Q7,8|K 0 in the chiral limit written as a K + → π + matrix element according to Eq. (15), to avoid ambiguities due to normalization conventions for the two pion state. All results are given in the RI/MOM scheme at 2 GeV. Please see text for a discussion of how the central value and systematic error on the matrix elements is obtained.   [20] This work [20] This work [20]

VII. FINAL PHYSICAL RESULTS AND DISCUSSION
As seen the preceding Section, results at β = 6.0 and 5.85 are in excellent agreement. We could therefore assume that we are in the scaling regime and fit the a 2 dependence of our results to a constant. However, as discussed briefly in Sec. V, the determination of the renormalization constants is difficult at β = 5.85. We therefore choose to take the results obtained on our finer lattice as our final central values and use the results at β = 5.85 to estimate possible discretization errors. The lattice spacing squared at β = 6.0 is essentially half way between that at β = 5.85 and the continuum limit. Assuming that a linear extrapolation of the results in a 2 to the continuum limit gives a worst case scenario, we estimate our symmetric discretization errors to be the absolute value of the difference of the central values of the results obtained at the two lattice spacings. We then add this error in quadrature with the error coming from the systematic uncertainty in the renormalization constants. The procedure assumes that O(a 2 ) errors are dominant over higher order discretization errors and that their coefficient is the same at the two values of β, as discussed in Sec. II.
Let us first consider the BSM matrix elements. We give our final results for the B-parameters, couplings G i and BSM ratios R i in Table VII, in the RI/MOM scheme at 2 GeV. In this table we also provide the results derived from the only other calculation of these matrix elements [20]. The latter were obtained in the quenched approximation, using tree-level, O(a)-improved Wilson fermions. The first point that we wish to make is that our result for B K = B 1 is in excellent agreement with the recent world averages of lattice results given in [32][33][34]. The second point is that our BSM matrix elements are enhanced at M K by a factor of five to twenty compared to the chirally suppressed standard model one, as shown by the values that we obtain for the BSM ratios R i . The third point is that this enhancement is significantly larger than the one found in [20]. This difference can only partially be explained by our smaller value of B 1 , which accounts for about a 23% enhancement in our results for R i over those of [20]. The remainder of the enhancement must therefore come from elsewhere. It does not come from the B-parameters B i for i = 2, · · · , 5, which agree within errors. Since the B-parameters agree much better than do the matrix elements themselves, we must already disagree on the vacuum saturation values, that is on the matrix element of the pseudoscalar density, 0|sγ 5 d|K 0 . Since the axial ward identity requires that 0|sγ 5 d|K 0 = √ 2F K M 2 K /(m s +m), we must in fact find a smaller (m s +m). Using the results of Tabs. I-LV and the method described above for obtaining final results on the matrix elements, we find, for the quark masses: which is in excellent agreement with the earlier Neuberger fermion calculations of [53,69]. Our result for this sum  of quark masses is almost 28% lower than the one found by the authors of [20], using tree-level, O(a)-improved Wilson fermions at the couplings which they use to obtain the BSM matrix elements. Indeed, in [70], they obtain (m s +m) MS (2 GeV) = 130(2)(18) MeV with the same gauge configurations as in [20]. Since this mass appears squared in the relation of the BSM matrix elements to their respective B-parameters, this difference accounts for the large enhancement of our BSM matrix elements over those of [20] and, combined with our lower value of B 1 , for the factor of roughly two enhancement of our BSM ratios R i . Note that the result of Eq. (34) is compatible with the rigorous lower bounds calculated in [71,72]. For completeness, we determine the value of the strange quark mass. Using m s /m = 24.4 ± 1.5 as obtained from chiral perturbation theory [73], the results in Eq.
where the uncertainty in the chiral perturbation theory mass ratio does not contribute significantly to the error. Our result is fully compatible with the quenched, continuum limit, benchmark result of the ALPHA collaboration, m MS s (2 GeV) = 97.(4) MeV, obtained using F K , as we do, to set the lattice spacing [74]. We therefore believe that the results of [20] for the BSM matrix elements suffer from large systematic errors which are absent in our calculation.
Because the MS-NDR scheme is frequently used, we also provide our results for the B-parameters, couplings G i and BSM ratios R i in this scheme at 2 GeV in Table VIII.
We turn now to the matrix elements of the electroweak penguin operators, (ππ) I=2 |Q 7,8 |K 0 . We compute them in the chiral limit, as discussed in Sec. VI. We obtain final results, as for the BSM matrix elements and quark masses, by combining our results at β = 6.0 and 5.85 in the manner discussed above. These results are summarized in Table IX, both in the RI/MOM and MS-NDR schemes at 2 GeV.
These matrix elements have already been obtained, both in quenched lattice QCD and with a variety of analytical techniques. In Table X, we compare our results with those obtained by others in the NDR scheme at 2 GeV, also in the chiral limit. The results quoted for CP-PACS [42] were obtained by multiplying their results for (ππ) I=2 |Q 7,8 |K 0 at 1.3 GeV in the MS-NDR scheme by − 2/3(F χ + F K )/(2F π ), to match our definition of the matrix elements. They were then run up from 1.3 to 2 GeV using N f = 0 two loop running with the parameters discussed in Appendix E.
Those of RBC [43] were obtained from the expression π + |Q 3/2 8,8 (7,8),0 /f 2 ) × a −4 with f = 137(10) MeV [75]. These bare results were then renormalized in the RI/MOM scheme with the renormalization constants given in the paper at 2.13 GeV, matched onto the MS-NDR scheme and finally run down to 2 GeV with our two-loop running formulae. Finally, the results from SPQcdR [45] are obtained by extrapolating, to the chiral limit, quenched results for K → ππ matrix elements computed for a variety of initial and final state kinematics, using NLO chiral perturbation theory. X: Comparison of our results for the matrix elements of the electroweak penguin operators in the chiral limit, with those obtained in other lattice calculations (above the double line) and with a variety of analytical techniques (below the double line). The "??" indicate that we have not estimated systematic errors in converting the results to our conventions. Moreover, the factors of 2/3 in front of the SPQcdR results [45] are required, we believe, to translate their results to our normalization. a Results are given in the MS-NDR scheme at 2 GeV.  [76] 0.12(2) 2.00(36) Knecht et al '01 [21] 0.11(3) 3.5(1.1) Narison '00 [38] 0.21 (5) 1.40 (35) a We thank C.-J. David Lin for correspondence on this subject.

Ref. Action
Agreement of our results with lattice calculations performed using domain-wall fermions is excellent for both matrix elements. Only the Wilson fermion results of [45] differ significantly. This may be due to the absence of chiral symmetry in the Wilson discretization which significantly complicates the renormalization procedure and allows potentially large O(a) discretization errors. In regards to the results obtained using analytical techniques, agreement is found with the results of [38,39] and to a lesser extent [40]. The authors of [21,76], however, find results for the matrix element of Q 3/2 7 which are a factor of two lower than ours, and for the matrix element of Q 3/2 8 , a factor of two to three higher. It should be noted that a quenched investigation of electroweak penguin matrix elements with modified overlap fermions, and perturbative renormalization, has been performed in [44]. However, the functional form for the chiral extrapolation of the matrix elements is incorrect and at the kaon mass, the results obtained at two different lattice spacings do not appear to be compatible. In addition, preliminary results from a quenched calculation using HYP staggered fermions were presented in [77]. Finally, SPQcdR has also attempted to determine the (ππ) I=2 |Q 7,8 |K 0 matrix elements at NLO in chiral perturbation theory by fitting the required low energy constants to quenched lattice results for a variety of initial and final state kinematics [45].

VIII. CONCLUSION
We have presented results for the full set of ∆S = 2 matrix elements which are required to study K 0 -K 0 mixing in the standard model and possible extensions. We have also presented results, at leading chiral order, for the matrix elements of the electroweak penguin operators Q 3/2 7,8 , which give the dominant ∆I = 3/2 contribution to ǫ ′ . As a by-product of our calculation, the strange quark mass was obtained. Quarks were simulated with Neuberger fermions, which possess an exact chiral flavor symmetry at finite lattice spacing. This not only greatly simplifies the complicated mixing pattern under renormalization of the operators that we consider, but also guarantees that our results are free of the leading O(a) discretization errors. To avoid potentially large perturbative errors, we implemented all renormalizations non-perturbatively in the RI/MOM scheme, modifying the techniques of [16] to accommodate non-perturbative power corrections and high-momentum discretization errors. For most renormalization constants, we find that these effects are important. The calculations were performed on two sets of quenched configurations, generated with the Wilson gauge action on an 18 3 × 64 lattice at β = 6.0 and on a 14 3 × 48 lattice at β = 5.85. The two lattice spacings allowed us to investigate scaling violations. Within our statistics, we find no evidence for such violations and use the results on our coarser lattice to estimate possible discretization errors.
Our main conclusion is that the non-SM, ∆S = 2 matrix elements are significantly larger than found in the only other dedicated lattice study of these amplitudes [20]. In tracing the source of this difference, we found that we already disagree on the much simpler matrix element of the pseudoscalar density between a kaon state and the vacuum, which is the building block for the vacuum saturation values of the BSM ∆S = 2 amplitudes. Through the axial Ward identity, the matrix element of the pseudoscalar density is related to the sum of the strange and down quark masses, which we find to be roughly 30% smaller than the value obtained in [70], with the same tree-level improved Wilson fermion action and gauge configurations as used in [20]. Since our result for this sum of masses is in agreement with the continuum limit, benchmark result of [74], we are convinced that the stronger enhancement of non-SM ∆S = 2 matrix elements that we observe is correct.
Regarding B K and the matrix elements of the electroweak penguin operators Q 3/2 7,8 , we find good agreement with recent quenched lattice calculations, in particular those which make use of domain-wall fermions [35,42,43,78], which have an approximate chiral symmetry at finite lattice spacing. Agreement with analytical results depends on the particular calculation.
Beyond this, our investigation validates the use of Ginsparg-Wilson fermions for the calculation of weak matrix elements. The simplified operator mixing and non-perturbative O(a) improvement guaranteed by the full chiral flavor symmetry of these fermions has proven essential for obtaining reliable results in this delicate calculation. We hope to have convinced the reader that these benefits outweigh the numerical overhead that accompanies the use of these fermions. Of course, our results suffer from the shortcomings of the quenched approximation, though we have tried to minimize their impact by considering dimensionless measures of the various matrix elements. Nevertheless, because of the numerous advantages already discussed, we believe that the way to remove this approximation is by studying Ginsparg-Wilson quarks on unquenched backgrounds generated with numerically cheaper fermion formulations, such as O(a)-improved Wilson fermions. Work in that direction is beginning. In this appendix, we have regrouped the results for the pseudoscalar meson masses and decay constants, in lattice units, which are discussed in Sec. III. Please refer to that section for explanations.  In this appendix, we have regrouped the results for the bare B-parameters, BSM ratios, "couplings" G i and ratios D 3/2 7,8 as a function of bare quark mass at β = 6.0 and 5.85. These results are discussed in Sec. IV. Please refer to that section for explanations.      In this appendix, we have regrouped the tables of renormalization constants discussed in Sec. V. Please refer to that section for explanations.            In this appendix, we have regrouped the results for the B-parameters, BSM ratios, "couplings" G i and ratios D 3/2 7,8 , renormalized in the RI/MOM scheme at 2 GeV, and as a function of bare quark mass at β = 6.0 and 5.85. These results are discussed in Sec. VI. Please refer to that section for explanations.       Because they may be of use in future calculations, we give non-perturbative results for the renormalization constants in the MS-NDR scheme. These are obtained from the RGI renormalization constants obtained with the fits described in Sec. V, through the relation Z NDR ij (2 GeV) = U NDR ik (2 GeV)Z RGI kj . Here, U NDR ik (p 2 ) describes the running of the renormalization constants in the MS-NDR scheme which, again, we implement at two loops [58] with N f = 0 and Λ QCD from [59]. The matching of the renormalization constants for the pseudoscalar density Z P = Z S is implemented at four loops [60][61][62], with N f = 0 and the same value of Λ QCD . Our results for all of these renormalization constants are given in Tabs. XXXVIII-XLVI, including our estimate of the systematic error associated with the analysis performed in Sec. V.          Because MS-NDR is a popular scheme, we give in this appendix results for the B-parameters, BSM ratios, "couplings" G i and ratios D 3/2 7,8 , in that scheme at 2 GeV, and as a function of bare quark mass at β = 6.0 and 5.85. While we do not use these results in our paper, they should be useful for comparisons with future lattice calculations done in the MS-NDR scheme. The results are obtained by combining the bare results of Sec. IV, which are tabulated in Appendix B, and the renormalization constants obtained in Appendix E, using the extended p 2 range.     Because the MS-NDR scheme is commonly used, we provide here results in the MS-NDR scheme at the individual lattice spacings for the B-parameters and other matrix element ratios at the physical point. These results are used in Sec. VII to obtain our final physical results in that scheme. They are obtained by multiplying the RI/MOM results obtained in Sec. VI by U NDR (2 GeV)[U RI (2 GeV)] −1 . Here, U NDR (p 2 ) and U RI (p 2 ) are the matrices which describe the running of the matrix elements of interest in the MS-NDR and RI/MOM scheme, respectively. In the case of the four-quark operators, we implement them at two loops [57,58] with N f = 0 and Λ QCD from [59], as in Sec. V. The running for the pseudoscalar density is implemented at four loops [60][61][62], with N f = 0 and the same value of Λ QCD . Our results for all quantities of interest in the MS-NDR scheme at 2 GeV and at the two values of the lattice spacing are summarized in Table LV. The systematic errors quoted in these tables are obtained in that same way as in Sec. VI.    (16) interpolated to the kaon point and our best estimate for (ππ)I=2|Q7,8|K 0 in the chiral limit written as a K + → π + matrix element according to Eq. (15), so as to avoid ambiguities due to normalization conventions for the two pion state. All results are given in the MS-NDR scheme at 2 GeV. Please see text for a discussion of how the central value and systematic error on the matrix elements is obtained.