Neutral Kaon Mixing Beyond the Standard Model with $n_f=2+1$ Chiral Fermions Part 1: Bare Matrix Elements and Physical Results

We compute the hadronic matrix elements of the four-quark operators relevant for $K^0-{\bar K^0}$ mixing beyond the Standard Model. Our results are from lattice QCD simulations with $n_f=2+1$ flavours of domain-wall fermion, which exhibit continuum-like chiral-flavour symmetry. The simulations are performed at two different values of the lattice spacing ($a\sim0.08$ and $a\sim 0.11 \, \fm $) and with lightest unitary pion mass $\sim 300\, \MeV$. For the first time, the full set of relevant four-quark operators is renormalised non-perturbatively through RI-SMOM schemes; a detailed description of the renormalisation procedure is presented in a companion paper. We argue that the intermediate renormalisation scheme is responsible for the discrepancies found by different collaborations. We also study different normalisations and determine the matrix elements of the relevant four-quark operators with a precision of $\sim 5\%$ or better.


I. INTRODUCTION
The investigation of neutral Kaon mixing has been an important area for our understanding of the Standard Model (SM) of particle physics.CP-violation was first observed in K S regeneration experiments [1] and the small value of the K L -K S mass difference led to the prediction of the charm quark at the GeV scale [2,3].Neutral Kaon mixing within the SM is dominated by W -exchange box diagrams as illustrated in Fig. 1.By performing an operator product expansion, one can factorise the long-distance effects into the matrix element K0 |O 1 |K 0 of the four quark operator.
where a and b are colour indices and the summation over Dirac indices is implicit.In the SM, the only Dirac structure which contributes is "(Vector-Axial) × (Vector-Axial) " arising from the W-vertices.The four-quark operator given in Eq. ( 1) is invariant under Fierz re-arrangement, therefore gluonic exchanges do not introduce new four-quark operators.
In a massless renormalisation scheme which preserves chiral symmetry the four-quark operator O 1 does not mix with other four-quark operators, nor with lower dimensional operators.The importance of the matrix element given in Eq. (1) has motivated many lattice studies of the SM kaon bag parameter (defined in some renormalisation scheme at some scale µ) which have now achieved accuracies at the few-percent level [4][5][6][7].(Our convention for the decay constant is such that f K = 156.1 MeV.) Combined with the value of the Wilson coefficient C(µ), computed in perturbation theory, and experimental observables, such as the mass difference ∆ M K = m K L − m K S and ε K , the determination of B K (µ) provides important constraints on the Cabibbo-Kobayashi-Maskawa (CKM) matrix.Schematically, one obtains where F is a known function of the CKM factors and of well-measured quantities.In the framework of the SM, the experimental value of ε K (which parametrizes indirect CP violation) together with the theoretical determination of B K provides an important constraint on the apex of one the CKM unitary trianglesin the (η, ρ) plane -and on the overall consistency of the CKM picture.ε K is also a powerful probe of potential new physics, with sensitivity to energies well beyond those being explored directly at the LHC (see for example [8][9][10][11]).
Beyond the SM, both left-handed and right-handed currents may contribute in the K 0 -K0 mixing process and the CP-violation parameter ε K is sensitive to new CP violating phases generically predicted by these models.Here we assume that the new-physics effects occur at energy scales much higher than the interaction scale of QCD and that QCD remains a valid description of the strong interaction in the non-perturbative regime.In addition to the SM operator O 1 given in Eq. ( 1), seven four-quark operators appear in a generic effective ∆S = 2 Hamiltonian [12] 1 where and Õi=1,2,3 are obtained from the O i=1,2,3 by swapping chirality (1−γ 5 ) → (1+γ 5 ).The Wilson coefficients C i (µ) and Ci (µ) depend on the details of the new-physics model under consideration but the matrix elements K0 |O i |K 0 are model independent.(In our framework parity is conserved, therefore the operators Õi=1,2,3 are redundant).In terms of representation of SU L (3) × SU R (3), it is straightforward to show that in the chiral limit O 2 and O 3 transform like (6,6) while O 4 and O 5 belong to (8,8).Therefore these new operators mix pair-wise under renormalisation: O 2 with O 3 ; and O 4 with O 5 .
In contrast to B K (µ), studies of the extended set of matrix elements are relatively few.The first computation performed with dynamical fermions was reported by our collaboration in [15] and was done with n f = 2 + 1 DW fermions at a single lattice spacing.It was followed by a n f = 2 computation by the European Twisted Mass (ETM) collaboration using twisted-mass Wilson fermions with several lattice spacings [16].These two computations reported results in decent agreement (the matrix elements of O 2,3,4 agree within errors, O 5 only within ∼ 2 σ), suggesting that these quantities are not very sensitive to the number of flavours.However, another study by the Staggered Weak Matrix Element (SWME) collaboration using n f = 2 + 1 flavours of improved staggered fermions [17] found a noticeable disagreement for two of these matrix elements (O 4 and O 5 ).The ETM collaboration has since repeated their computation with n f = 2 + 1 + 1 flavours and found bag parameters compatible with their n f = 2 results (only within ∼ 2 σ for O 5 ) [18].The SWME collaboration has extended their previous study by adding more ensembles and improving extrapolations to the physical point [19], they confirmed their disagreement with the other studies.Since the results have been extrapolated to the continuum limit, one does not expect the fermion discretisation used (Domain-Wall, Twisted-Mass, or Staggered) to be responsible for the discrepancy.
Central to this work is an explanation for this disagreement, our arguments and preliminary results have been presented in [20] and discussed with the authors of [19].We improve upon our earlier DWF result [15] in two important ways: by adding a second lattice spacing, allowing us to take the continuum limit (with a resonable handle on the lattice artefacts) and by renormalizing the four-quark operators through nonexceptional momentum schemes.
As we will show, the second point is of great importance and is often overlooked.Some systematic errors in the original RI-MOM schemes which are very hard to control are absent in the RI-SMOM schemes we present here.
In the next section, we give an overview of our strategy and make explicit our choice of conventions (choice of basis, normalisation).Sections II and III contain our global fit procedure and the method for determining the bare hadronic matrix elements K0 |O i (µ)K 0 .In section IV we present our final results and compare with previous works.

II. EXTRAPOLATIONS TO THE PHYSICAL POINT
In this work we have considered data with pion masses in the range of m P ∼ 300−430 MeV and performed a chiral extrapolation to the physical value of m π = 140 MeV (we take the mass of the charged pions).The spatial extent our the simulated lattice is L ∼ 2.66 fm, so within this range of pion masses Lm P > 4, therefore the finite volume effects are expected to be negligible compared to our systematic errors.We work in the isospin limit, m u = m d ≡ m ud and for the same reason we do not consider isospin corrections.
Furthermore, we also require a continuum extrapolation to reach the physical point (a = 0, m π = 140 MeV).
Since we work with Domain-Wall fermions, we expect the dominant lattice artefacts to be linear in a2 (we remind the reader that a 3 corrections of the fermionic action are forbidden by chiral symmetry 2 ).Before the continuum extrapolation can be performed, a renormalisation step is also necessary: we employ the non-perturbative Rome-Southampton method [21], as explained in detail in a companion paper [22].Below we list our strategy to extract the physical quantities of interest from our lattice simulations: 1. Compute the bare matrix elements, at two values of the lattice spacing and several values of the quark masses (on already existing RBC-UKQCD ensembles).
3. Interpolate/extrapolate to the physical value of the strange quark mass.
4. Extrapolate to the physical point (Continuum/Chiral extrapolation in the light quark sector).
Central to this work is an investigation of the extrapolations to the physical point (details can be found in section IV).In particular we have studied several parametrisations of the four-quark operator matrix elements.Ideally, one would like to find a dimensionless quantity which can smoothly be extrapolated to the physical point and be free of large systematic errors.For the SM matrix element one usually defines the bag parameter B K as in Eq. ( 2): The matrix element of the four-quark operator is normalised by its Vacuum Saturation Approximation (VSA).This normalisation is widely accepted for the SM contribution, however this is not the case for the BSM matrix elements, for which different possibilities have been proposed (see for example [17,[23][24][25]).We investigate several strategies which differ by the choice of normalisation and global fit procedure, allowing us to estimate the systematic uncertainties of our work.

A. The ratios R i
A possible parameterisation of the matrix elements has been proposed in [25].Denoting by P the simulated strange-light pseudo-scalar particle (kaon) of mass m P and decay constant f P , the ratios R i are defined by such that at the physical point is the ratio of the BSM matrix element to the SM one.Previous studies have shown that these ratios are large ( ∼ O(10) ) as the BSM matrix elements are enhanced compared to the SM one [16,25,26] (this is expected from Chiral Perturbation Theory: the SM matrix element vanishes in the chiral limit whereas the BSM matrix elements remain finite).An advantage of this method compared to the bag parameters is that the denominators do not depend on the quark masses.The BSM matrix elements can be reconstructed from the ratios R i , the SM bag parameter B K , the kaon mass and decay constant only.Moreover, since the numerator and the denominator are very similar, one expects some cancellations of the statistical and systematic errors to occur in the ratio.

B. The Bag parameters B i
The renormalised bag parameters are defined as the ratio of the weak matrix elements normalised by their VSA values: For the SM bag parameter B 1 (µ) = B K (µ) with our conventions, and for the BSM ones3 , The factors N i>1 depend on the basis, as we work in the SUSY basis we have N i>1 = − 5 3 , 1 3 , 2, 2 3 .For SM bag parameter B K , the denominator consists of the precisely known quantities f K and m K .This contrasts with the BSM B i , for which the denominator is not uniquely defined, it depends on the scheme and the renormalisation scale.

C. The Combinations G ij
Another possibility, advocated for example in [17,24] is to define products and ratios of bag parameters such that the leading chiral logarithms cancel out.For some of these quantities (called "golden combinations"), this cancellation actually occurs at every order of the chiral expansion.For the other ones ("silver combinations"), only the leading logarithms cancel.Such quantities were introduced in [24] for SU (3) chiral perturbation theory and later in the context of SU (2) staggered chiral perturbation theory in [27].The relevant NLO continuum SU (2) chiral expansions can be found in Appendix V C. We follow [17] and define4 As can be seen in the Appendix V C, the quantities G 23 and G 45 have no chiral logarithms, whereas in G 21 , G 24 the cancellation only occurs for the leading logarithms.

D. Continuum and chiral fitting strategies
We start by adjusting our (renormalised) results to the physical strange mass.On the coarse lattice we perform a linear interpolation whereas a tiny extrapolation is necessary on the fine one (the numerical values are given in the next section).Then we perform a combined chiral-continuum extrapolation to the physical point.In order obtain a reliable estimate of our systematic error we follow three different strategies: • Method A. We perform a global fit according to NLO SU (2) chiral perturbation theory (see Appendix V C).The general form of the fit function we use is (we drop the renormalisation scale dependence µ for clarity), Where in this expression m P is the mass of the pseudoscalar meson made of two light quarks.The values Y i , α i and β i are free parameters and fit simultaneously between ensembles of different lattice spacings.The values for C i are listed in Table I below.We have checked that for f , using the chiral value, the physical value or the simulated value f P give compatible results.We apply the procedure to the ratios R i and to the bag parameters B i .
TABLE I: Chiral logarithm factors C i for R i and the B i .
• Method B. We perform a continuum/chiral extrapolation of R i and B i using a global fit procedure according to the following ansätz (κ i and δ i are free parameters simultaneously fit between ensembles) • Method C. We first extrapolate the combinations G ij according Method B (linearly in the pion mass squared), and then reconstruct the bag parameters.
Methods A and B are equivalent up to the chiral logarithm terms, the difference allows us to estimate how strong the chiral effects from being at non-physical pion mass are.The corresponding analysis is presented in great detail in section IV.Method C allows us to determine the bag parameters with no leading chiral logarithm, except from the standard model one, whose effect is benign (as explained in section IV).
Furthermore, the quantities G ij have different statistical and systematic errors.Performing the analysis using different quantities and extrapolation methods allows us to check the consistency of our final results and ensure our systematics are understood.The results for Method C are presented in Appendix V F.

III. LATTICE IMPLEMENTATION
Our measurements are performed on n f = 2 + 1 gauge ensembles generated by RBC-UKQCD using the Iwasaki gauge action [28,29] and the Shamir DWF formulation [30].These ensembles have been described extensively in [31] and references therein.
The finer of the two lattices used in this study has a lattice volume of 32 3 × 64 × 16 with inverse lattice spacing a −1 = 2.383(9) GeV.There are three values of light sea quark masses am sea ud = 0.004, 0.006, and 0.008, corresponding to unquenched pion masses of approximately 300, 360, and 410 MeV respectively.
For the light valence quarks we use only unquenched data, am val ud = am sea ud .The simulated strange quark mass for this ensemble is am sea s = 0.03.To reach the physical kaon mass we extrapolate using unitary (am val s = am sea s = 0.03) and partially quenched (am val s = 0.025) data, which is close to its physical value of 0.02477 (18) [5].
The coarser lattice has an extent of 24 3 × 64 × 16, and inverse lattice spacing a −1 = 1.785(5)GeV.There are three values of light sea quark mass used in the simulations, am sea ud = 0.005, 0.01 and 0.02 (we drop the heaviest of these in the chiral extrapolations).We again use only unquenched light valence quarks, corresponding to pion masses of approximately 340 and 430 MeV.The simulated strange quark mass for the ensemble is am sea s = 0.04, while the physical value has been determined to be am phys s = 0.03224 (18).
As with the fine ensemble, we interpolate between unitary (am val s = am sea s = 0.04) and partially-quenched

A. Correlation functions
We have used Coulomb gauge fixed wall-source propagators, which allow for much greater statistical resolution at similar cost to a point-source propagator inversion and should have better overlap of the ground state.The fine ensemble results were generated as part of the calculation of B K in [4].The coarse ensemble configurations were first Coulomb gauge fixed using the time-slice by time-slice FASD algorithm of [32] (to an accuracy of Θ < 10 −14 ).
Working in Euclidean space, we define the two-point functions, where O i represents a bilinear operator.For the present analysis we only consider non-flavour singlet operators with two different Dirac structures: either P the pseudo-scalar density, or A 0 the temporal component of the local axial current.The superscripts (s i ) denote the source type, either (L)ocal or (W)all source.The two-point functions are fit to their asymptotic form (T is the temporal extent of the lattice): Our conventions are such that and P = ψ1 γ 5 ψ 2 (and therefore P = ψ2 γ 5 ψ 1 ) denotes a (flavour non-singlet) pseudo-scalar sate of mass m P .
The corresponding decay constant f P is defined (at finite lattice spacing and zero momentum) by and can be extracted from an appropriate ratio of two-point functions.The superscript R denotes the fact that a finite (re)-normalisation factor is required to connect the local axial current We prefer to renormalise the axial current with Z V rather than Z A for numerical reasons, (Z A and Z V should be identical if chiral symetry is exact, however Z V is numerically easier to extract).In a similar way, the bare matrix elements P|O i |P are determined from three-point correlation functions where the operator is inserted between two well separated wall sources, In order to have a better handle on our systematics, we extract the quantities of interest in different ways (which are in principle equivalent up to lattice artifacts).Our key results are obtained through the ratio of three-point functions (k = 2, . . ., 5) which we fit to a constant in the asymptotic region: We also define the ratios of three-point over two-point functions, which at large times allows us to obtain the bare BSM bag parameters: We show some examples of plateaux in Figs 2 and 3.The simulated time extent is T /a = 64 on both lattices, but for the fine lattice we implement the Periodic ± Anti-periodic trick which is designed to reduce the round the world artifacts.Effectively this trick doubles the number of accessible points [33] (see also the discussion in [31]).Although the signal obtained from the coarse lattice time slice per time slice is different from the one of the fine lattice, the precision obtained on the ratio R Lat k (by a correlated fit) is of the same order.

B. Non-Perturbative Renormalisation (NPR)
Once the bare matrix elements have been obtained, they need to be renormalised in order to have a well-defined continuum limit.We opt for the framework which is now standard within the RBC-UKQCD FIG. 2: Example of the plateau for R Lat i (T, t, 0) as a function of the operator insertion time t/a.We show our results for the lightest kaon mass on our coarse lattice.collaboration: the non-perturbative Rome-Southampton renormalisation method [21], with non-exceptional kinematics (we use the symmetric RI-SMOM schemes) [34], momentum sources [35] and twisted boundary conditions [36][37][38].Similarly to what was done for B K and K → ππ, we define two schemes: the RI-SMOM-(γ µ , γ µ ) and RI-SMOM-( / q, / q) schemes5 (we drop the "RI" in the following).We refer to these schemes as "intermediate schemes".Our final results are the ones given in these SMOM schemes; however the matrix elements of interest are conventionally given in a MS scheme at a reference scale of 2 or 3 GeV.Although the computation of the bare matrix elements and of the renormalisation factors is done non-perturbatively, FIG. 3: Same as the previous figures but for our fine lattice.
this matching step involves (continuum) perturbation theory.MS results obtained via different intermediate schemes should be consistent, up to higher-order PT matching corrections (and lattice artifacts if the resutls are given at finite lattice spacing).The use of multiple intermediate schemes allows one to gain a better handle on these uncertainties6 . .
We also implement the original RI-MOM scheme [21], however we find that the results are not consistent with the SMOM ones.As shown in detail in the companion paper [22], we find that the RI-MOM Z-matrices exhibit large violations of the block diagonal structure expected from the chiral-flavour properties of the four-quark operators.This seems to be due to important infrared artefact which go as inverse powers of the quark mass.These pole "contamination" require a hard subtraction, and render results significantly more unreliable than their SMOM counterparts.This is indicated in Table VI by the discrepancies of the RI-MOM scheme results with the SMOM ones and with the ones obtained by the SWME collaboration, whose (1-loop) perturbative matching is free from IR contamination (see [19,20]).We do not advocate the use of these RI-MOM results, indeed we show that this choice of intermediate scheme is probably the cause of the disagreement observed between different collaborations.
Another advantage of the SMOM schemes is that the perturbative matching factors connecting them to the MS scheme are much closer to the identity matrix.This suggests a better behaved perturbative series with less matching uncertainty than for the MOM case, which would demand a higher matching scale.
Referring again to Table VI, the close compatibility of the SMOM-(γ µ , γ µ ) and ( / q, / q) results provides strong evidence that the matching uncertainty for the SMOM schemes is negligible within our error budget.

IV. RESULTS AT THE PHYSICAL POINT AND DISCUSSIONS
We report here our main results for the ratios R i , the bag parameters B i and the combinations G ij .We consider the main results of this work to be the ratios R i , because at the physical point they directly provide the ratio of the BSM matrix element to the SM one.They do not depend on the quark masses, nor on our ability to renormalise the pseudo-scalar density as the bag parameters and some of the combinations G ij are.The results for the bag parameters extracted from the combinations G ij (method C) are reported in Appendix V F. We also compute the matrix elements K0 |O i |K 0 using the different strategies.The quality of the fits can be judged from the χ 2 reported in Appendix V A, Table IX A. The ratios R i In Fig. 4, we show the results using the combined continuum-chiral fits discussed in section II, both Method A and Method B in the non-exceptional SMOM-(γ µ , γ µ ) scheme.We show all of our results in this scheme, however we note that the SMOM-( / q, / q) scheme gives very similar results.The RI-MOM results have already been presented in [15,40], they are just reported for comparison with previous work.In the figures, the dashed line represents the chiral extrapolation performed linearly in m 2 P (the pion mass squared) at fixed lattice spacing and the a 2 → 0 extrapolation is shown as a solid black line.The magenta line represents the Method A fit, in which we take the leading chiral logarithms into account.Our physical results obtained by this chiral-continuum extrapolation is the filled circle.
We note that the fit quality is very good with chi-square per degree-of-freedom (χ FIG.4: Continuum/chiral extrapolation Method A and B of the ratio R i in the SUSY basis and renormalised in the (γ µ , γ µ )-scheme.The conventions here and in the following plots are: red squares are the fine lattice data points, the blue squares the coarse ones.Open symbols represent a point which was omitted in the fit procedure.All the points have been interpolated/extrapolated to the corresponding physical strange quark mass.The magenta curve is the chiral fit and the solid point is its chiral-continuum value.The black line is the linear fit at a 2 = 0. We keep the relative scale constant for the vertical axis (around fifty percent of the extrapolated value).
or less as shown in Table IX of Appendix V A. We also note that although the ratios R i have the largest coefficients for the chiral logarithms, the effect of these terms is mild and the difference between the linear fit in m 2 P and the chiral one is at most of the order of a few per cent.The fits for Method A and Method B show similar quality as indicated by by their χ 2 /d.o.f , hence we do not see significant curvature.We take the fact that the fit quality for Method A is good as an indication that NLO Chiral Perturbation Theory is a decent description of the mass dependence of our data, this is our choice for our central values.We use the difference of the results obtained from Methods A and B to estimate the effects of the chiral logarithms.
As shown in the plots, the two methods give very close results.This might be because the ensembles we have used are at relatively heavy pion mass.However we give another argument below based on Method C, to justify that the chiral extrapolations to the physical quark masses are well under control, and that the chiral extrapolation effect is one of the most benign compared to the other systematics in this calculation.
For some of these quantities we see significant cut-off effects, especially R 5 which requires an extrapolation of the order of 15% from the fine ensemble's data to reach a 2 = 0. We observe that this is largely due to the 3-GeV renormalisation factors (for this quantity our estimate for the discretisation error is almost a factor two smaller at 2 GeV).From Fig. 4 it is interesting to note that as we approach the continuum limit R 2 , R 4 and R 5 of our BSM matrix elements are larger (in magnitude) than we previously determined just from our fine ensemble's data in [15].As other previous studies have noted, the BSM matrix elements are an order of magnitude larger than the SM one.

B. The bag parameters B i
The combined chiral-continuum plots for the B i are shown in Fig. 5 using the same conventions as in the previous section.We show our results again for the (γ µ , γ µ ) scheme.We observe that the fit quality is a bit worse for the B i compared to the R i with χ 2 /d.o.f ranging between 0.4 to 1.9 (Table IX).We also note that while the effect of chiral logarithms is almost invisible, the discretisation effects are larger than anticipated for two of these quantities B 3 and B 5 : we observe a deviation of more than 10% between the fine ensemble and the a 2 → 0 extrapolation.

C. The combinations G ij
Fig. 6 shows the results obtained in the (γ µ , γ µ ) scheme, using the same conventions as in the previous figures.Firstly, we see that there is no noticeable chiral curvature which is unsurprising as these quantities were designed for this purpose.We observe that the combinations G ij can be numerically very different.
For G 23 and G 45 , we find a rather good χ 2 /d.o.f., a linear behaviour in m 2 P (with a very small slope), however the lattice artefacts for G 45 are clearly visible (with again a difference of order 10% between the fine lattice and the extrapolated value).We have also computed an alternative combination, G23 , in order to compare our results with the SWME collaboration.Similarly to G 23 , it is defined as the ratio of the two bag parameters B 2 and B 3 , but computed in a different basis, the one introduced by Buras, Misiak, and Urban in [41].We call this basis the "BMU basis" in the following.This is also the choice of the SWME collaboration, therefore what we call G23 here is called G 23 in [17] and [19].Only B 3 differs between the two sets of operators.Within our convention the operator O 3 is defined as the colour partner of O 2 , whereas in the BMU basis, it is purely a "tensor-tensor" operator.Although in principle the two definitions are equivalent (thanks to Fierz theorem), the cutoff effects can be very different.Indeed we observe that the sign of the a 2 coefficient of two-colour partner operators are identical : positive for B 2 , B 3 and negative for B 4 , B 5 .This results in some cancellation of these artefact in the ratio G 23 (whereas in G 45 , the cutoff effects are completely dominated by B 5 and taking the ratio does improve very much from that point of view).We now turn to G23 , which reads in terms of bag parameters where B 2 and B 3 refer to the SUSY basis.In this peculiar combination, the cutoff-effects do not cancel, but on the contrary they add up, as illustrated in Fig. 7.We note that the authors of [19] also found this combination difficult to fit.
The fit of the product G 24 is very reasonable with a χ 2 /d.o.f. of around 1.4, the pion mass dependence is very mild, and there is clearly an important cancellation of the lattice artefacts in the product as the a 2 coefficients have a different sign.However, we believe that this cancellation is purely accidental.
We find that the ratio G 21 is much more difficult to fit, with a χ 2 /d.o.f of order seven.The difficulty comes mainly from the coarse ensemble, where the results seem to fluctuate around a constant value of the mass.
This effect could be due to some unfortunate statistical fluctuation or lattice artefact and need to investigated further in the future.This is rather unfortunate because the quantity this section).In the same Appendix V F, we compare the results for the bag parameters extrapolated directly (Methods A and B), to the ones extracted from the combinations G ij .We find that the combinations G do not provide more precise results (within our sytematic error budget) except for one quantity, B 3 (if G 23 is computed in the SUSY basis).
Finally, we point out that one could also first perform a continuum extrapolation of the bag parameters in the range of simulated pion mass, then compute the combinations G ij and finally perform the chiral extrapolation.We leave this for future investigations.

D. Error budget
Our central results are the BSM quantities non-perturbatively renormalised through the SMOM-(γ µ , γ µ ) and ( / q, / q) schemes, given in Tables III and IV.For these quantities, we have identified two main sources of systematic error: discretisation effects and chiral extrapolation to the physical pion mass.We have illustrated that some of our results have larger than expected O(a 2 ) lattice artefacts; since we have only two lattice spacings, we take half the difference between the fine ensemble's result (extrapolated to the physical pion mass) and the continuum extrapolation's result as an estimate of a potential curvature due to O(a 4 ) artefacts 7 In the future, it will be crucial to include a third lattice spacing to reduce (or eliminate) this 7 Exact chiral symmetry would guarantee the absence of O(a) and O(a 3 ) artefacts.Strictly speaking with Domain-Wall fermions there could be O(amres) and O((amres) 3 ) terms, however all our numerical studies show that these terms are numerically irrelevant, if not absent, as expected from naive power counting.
error and check that these quantities approach their continuum values linearly in a 2 .
Our chiral extrapolations are well under control, as illustrated in Figs. 4 and 5.We find that both a chiral perturbation theory prediction (Method A) and a linear Ansatz (Method B) in m 2 P give very good χ 2 per degree-of-freedom.We take half the difference between these to estimate our chiral extrapolation error.
We also observe that the results of the bag parameters extrapolated with a chiral fit give are very similar to those obtained from the combinations G ij (Method C), see Appendix V F. Since the combinations G ij are free from leading chiral logarithms we conclude that the chiral extrapolation to the physical quark masses are well under control.In the future, we plan to perform the computation at physical values of the quark mass [5] and therefore eliminate this error.
In Tables III and IV, we give the breakdown of our error budget.For our main results, the ratios R i renormalised in SMOM-(γ µ , γ µ ) and ( / q, / q) schemes at µ = 3 GeV, we give the statistical errors together with our estimate of the discretisation and chiral errors.We emphasise that these quantities are completely non-perturbative.We determine R 2 and R 3 with a precision better than 5%, whereas R 4 and R 5 have an error of 5% and 8% respectively.The latter are largely dominated by the discretisation errors, therefore we expect an important improvement with the future inclusion of a third lattice spacing in our analysis.
We have also converted our results to MS; since this matching is done in perturbation theory, there is an uncertainty due to the truncation of the perturbative series, in this case of order O(α 2 s ).We estimate this error by taking the difference: In Tables III and IV, this error refers to as "PT" (Perturbation Theory).Although the conversion can be done in the continuum limit, we checked that the applying the conversion to MS on the data before continuum/chiral extrapolation give the same results as if we apply it in the continuum (the difference is smaller than our statistical errors).In Table III, these results are denoted by (MS ← SMOM).We observe that the matching has very little effect on the central values and on the error budget (except of course that there is a perturbative error in addition).For the central value and the errors given in Tables III and IV, we quote the results obtained using SMOM-(γ µ , γ µ ) as an intermediate scheme 8 .The effect of the intermediate SMOM scheme is less than 3% for µ = 3 GeV and 4 − 5% for µ = 2 GeV.Regarding the total error, we find that all together, after conversion to MS, the µ = 2 GeV results are of the same size as the 3 GeV ones.(Although we also note that in general if we lower the scale, the perturbative errors increase and the discretisation errors decrease, as expected).
Scheme   We also give the error budget for the bag parameters B and their combinations G.Not surprisingly, we also find that the discretisation effects are larger than anticipated.In particular for the quantities B 3 and B 5 we quote an error of ∼ 8% and ∼ 5% at µ = 3 GeV.Clearly these errors come mainly from the NPR procedure as we observe a reduction of a factor two when we lower the scale to µ = 2 GeV.However, as for the ratios R i , the perturbative errors increase if we lower the scale and -apart from B 3 -we observe that the µ = 2 and µ = 3 results have similar total uncertainty, after conversion to MS.We expect the systematic uncertainty associated with the discretisation effects to drop drastically in the future with the inclusion of a third (finer) lattice spacing.The µ = 3 GeV results should then have have significantly reduced systematic errors in comparison to results renormalised at µ = 2 GeV

E. Final results and comparison with previous works
We report our final results for the ratios R, the bag parameters B and the combinations G in Table V.
The first error is statistical and the second combines the various systematic errors.Our main results are those given in the intermediate SMOM-(γ µ , γ µ ) and ( / q, / q) schemes.The RI-MOM results are only given for comparison with previous work.All these results are purely non-perturbative.The corresponding correlation matrices are given in Appendix V G.
For completeness, we also give our results after conversion to MS; in order to keep track of the intermediate scheme dependence, we denote them by MS ← scheme, where scheme can be one of the three intermediate schemes.We remind the reader that this conversion is done in perturbation theory, therefore the systematic errors also include an estimate of the perturbative error (except for the RI-MOM scheme as we do not find these results to conversion to MS at µ = 3 GeV, one expects the results to be independent from the intermediate scheme, up to small perturbative corrections.Table V shows that upon matching to MS the conversion has very little effect on the ratios for the non-exceptional schemes.Furthermore we the MS ← (γ µ , γ µ ) and MS ← ( / q, / q) are compatible within statistical fluctuations (in the worst case within ∼ 1.5 standard deviations).This is highly suggestive that the perturbative series for these schemes are well-behaved at this matching scale.
However, as shown in Tables III, IV and V, we observe that our new results using the non-exceptional schemes differ significantly from the ones renormalised though the RI-MOM scheme.This could be due to large higher order terms in the perturbative series for the matching of RI-MOM to MS that we neglect, although at the high matching scale we use this seems unlikely, which leaves this discrepancy to being some systematic inherent to the exceptional scheme renormalisation technique itself, such as the subtraction of the Goldstone pole (absent in the SMOM schemes).We argue below that the non-perturbative renormalisation procedure is the cause of the disagreement between the different collaborations and that it is due to systematic errors inherent in the RI-MOM scheme.
We finalise this section with a comparison of our results with previous measurements shown in Table VI.
We report the two most recent results of the ETM collaboration, who renormalised their results nonperturbatively using the intermediate, exceptional, RI-MOM scheme.We also compare our results to those of the SWME collaboration, who used 1-loop continuum perturbation theory.We choose to compare the bag parameters because the ratios R i are in general not reported by these collaborations.First, looking MS ← ( / q, / q) 0.536(9)(6)(11) 0.492 (7)  When present, the third one is the perturbative error coming from the matching to MS.Our best results are the ones obtained through the SMOM-(γ µ , γ µ ) and ( / q, / q) schemes.The RI-MOM results are presented here only for illustration and comparison purposes, we did not attempt to estimate the perturbative error for the MS ← RI − MOM case.
at the first three columns, ETM 12, ETM 15, and RBC-UKQCD 12, we see that the n f = 2 results are compatible with the n f = 2 + 1 and n f = 2 + 1 + 1 ones (only within ∼ 2.8σ for B 5 ), suggesting that these quantities do not depend strongly on the number of flavours.However the values of B 4 and B 5 quoted by the SWME collaboration differ significantly from the other determinations.In this work we show that our  When only one error is quoted, it means that the errors have been already combined.If not, the first errors are statistical and the second systematic.We argue that the renormalisation procedure is the cause of the disagreement observed for B 4 and B 5 between the different collaborations and that it is due to some underestimated systematic errors present in the RI-MOM scheme.For the RI-SMOM results, we choose the (γ µ , γ µ ) scheme.
values of B 4 and B 5 are compatible with those of the ETM collaboration if we use the RI-MOM intermediate scheme.However, if we use an SMOM scheme (as we strongly advocate in this work) our results are then compatible with the SWME collaboration .The fact that we are compatible with ETM whilst using the same renormalisation scheme suggests that the scheme dependence we see is legitimate.

F. Matrix elements of the BSM four-quark operators
We end this section with the matrix elements of interest K0 |O i |K 0 .They can obtained from the ratios R i , the bag B i or the combinations G ij with different source of systematic errors.We find that the methods give consistent results, but the error can be very different.We find that the (6,6) operators are more precise when computed from the R i whereas the bag B i give smaller systematic errors for the the (8, 8) operators.Our best estimates are given in Table VII, where we also convert to MS for the reader's convenience.The corresponding correlation matrix is given in Table VIII.As expected, there are important correlations between operators of same chirality which have to be taken into account in phenomenological application.The non-perturbative results are obtained with a precision of 5% or better, this is the most precise computation of these matrix elements.The details of this computation are given in Appendix V E.  Since the combinations G ij are designed to cancel (at least) the leading chiral logarithms, we did not perform a chiral fit on these quantities.The results presented here are for the fits performed on quantities renormalised in the RI-SMOM and RI-MOM schemes.In Table IX we give the χ 2 per degree-of-freedom of the global fit used in Method A, B and C (see Section II).Method A corresponds to fitting the ratios R i or the bag parameters B i using Chiral Perturbation theory (χPT).Method B uses a linear fit in m 2 P , where m P is the simulated pion mass.We find that the fit of the ratios R i are of very good quality with a χ 2 per-degree-of freedom of order 0.5.The fits for B 2 and B 3 are a bit worse, although the χ 2 are still reasonable (of order 1.5).It is important to stress that our data do not seem to prefer either of the method, the effects of the chiral logs are not statistically significant.
We also show the χ 2 for the linear fits of the combinations G ij , (Method C).There we find that G 21 is very hard to fit, with a χ 2 per degree of freedom of order 6.We can see from Fig. 6 that the problem seems to come from the coarse ensemble and could be due to some lattice artefacts.The other G ij have a much more reasonable χ 2 .

B. Renormalisation factors
We give the Z matrices obtained though the SMOM-(γ µ , γ µ ) and SMOM-( / q, / q) schemes, together with their conversion to MS in Tables X-XV.

C. Chiral extrapolations
We only consider physical (unitary) quarks, so m val = m sea .We use the following notation     (/ q, / q) TABLE XV: Z/Z 2 V matrices for the (8,8) operators at µ = 3 GeV on the fine lattice, a = a 32 .
such that at leading order (LO) The parameter B related to the chiral condensate should not be confused with the bag parameter (noted B in this appendix).We consider kaon SU (2) L × SU (2) R χPT, ie m u = m d m s , Λ QCD .At next to leading order (NLO) we have Denoting the matrix element K0 |O i |K 0 by O i , we have thus for the Standard Model matrix element, we find We now turn to the BSM operators (O i>1 ) in the SUSY basis.They read Rewriting Eq. ( 26) The expansions for the Bag parameters read where s 2,3 = −1 and s 4,5 = 1.
It is then clear that the combinations, have no leading order chiral logarithm.
For the matrix elements of the operators, we obtain the following expansions: Finally we consider the ratios R i : We note that the chiral logarithms in R 2 and R 3 have the same coefficients as in R 1 .For completeness, we also give the following expressions:

D. Bare Results
Tables XVI and XVII show the fit results for the ratios of bare three-point function as described in Section III.These quantities are obviously correlated, not only they have been computed on the same gauge ensembles, but they are normalised by the same quantity.Furthermore O 2 and O 3 only differ by their colour structure (and similarly for O 4 and O 5 ), hence one expects them to have similar statistical fluctuations.We find that the correlations depend very mildly on the quark masses, so we only give the correlation matrices for the lightest unitary ensembles.The numerical values can be found in Tables XVIII and XIX.We also   The SM matrix element is computed from For the BSM matrix elements (i > 1), we can either use the ratios R i or the bag parameters B i In Eqs.37,38 and 39, we take m K = 495.6MeV, f K = 156.2MeV.For the value of B K , we take the results obtained in this work, but we checked that if we use the most recent value [5], the results are    The two methods give very consistent results.We also observe that for O 2,3 the ratios R i give more precise results, whereas for O 4,5 , the results obtained from the bag parameters B i have smaller error bars.With this choice, we obtain the matrix element with a precision of 5% or better.
The results for the quantities G ij given in III have been obtained by a linear extrapolation in m 2 π .Combining these results with the numerical value of B 1 , we can can reconstruct the BSM bag parameters (see 11).We observe that effect of the chiral logarithm for B 1 is negligible within our uncertainties.This has been confirmed by our recent computation in which physical quark masses are included [5].Here we find B (γµ,γµ) 1 = 0.523 (11) at µ = 3 GeV, in complete agreement with our new value B (γµ,γµ) 1 = 0.517 (2).Therefore the difference between the direct fit of the BSM bag parameters and the bag parameters reconstructed from   the quantities G ij is a direct indicator of the the chiral logarithms potentially present in the BSM operators.
Using Eq. ( 11) we find that For three of the BSM bag parameters, we implement an alternative strategy, called Method C .We define other combinations of bag parameter (also free of leading chiral logarithm) After extrapolation to the physical point, we extract the corresponding B 3,4,5 by inverting the previous system of equations.
A couple of remarks are in order • The difference between the various methods is smaller than our errors (actually smaller than the statistical error), showing that our chiral extrapolations are well under control within our precision.0.791 (11) stat (45) syst = 0.791(47) combined .However the corresponding matrix element is better determined from the ratio R 3 .
• As mentioned in Section IV, we have also computed G23 (which is denoted by G 23 in [19]), the results are shown in Fig. 7, We observe that G23 exhibit large a 2 lattice artefacts, see the discussion in section IV.Then B 3 can be computed from and we find B 3 = 0.767(82) (43) Not surprisingly, the error quoted here is much larger than the obtained from G 23 .Indeed, by changing the basis, the error varies by a factor two.

G. Correlations
To provide the correlations between measurements we compute the correlation matrix from our bootstrapped data.We represent this data visually by a matrix plot for our various measurement techniques, orange illustrates positive correlation and blue indicates anti-correlation, the darker the colour the stronger the correlation.Black squares are by definition 1.The analysis is done with 500 bootstrap samples.
In Fig. 8 we compare the correlations between the ratios R i for our two SMOM schemes.We observe that R 2 is strongly anti-correlated with all of the others due to the difference in sign with the others, this has operator signature SS − P P and most of the other ratios are strongly correlated with one-another.
We note that the correlations for the (γ µ , γ µ ) and ( / q, / q) schemes are very similar, this is in fact a feature for the other measurements so we will only show the (γ µ , γ µ ) scheme evaluations for the Bs and Gs.
Similarly to what we found for the correlation between bare ratios, we observe that the colour partners there is an relative sign in R 2 compared so the other ratios).However the correlation between operators of different chirality is enhanced compared to the bare rations (∼ 90%).As illustrated in Table XXIV, the correlation matrix does not depend on the renormalisation scheme.Although not shown here, the matching to MS also has almost no visible effects on the correlations.
For the bag parameters B i , we observe a similar pattern, however the correlations between B 4 and B 5 drops to 90%.The correlations between operators of different chirality are significantly lower than the ones for the ratios R i , namely around 50 − 60%.As expected, the quantities G ij do not exhibit significant correlation.

FIG. 1 :
FIG.1: Example of leading order box diagrams that contributes to K 0 − K0 mixing in the SM.

2 P 8 FIG. 7 :
FIG. 7: Alternative definition for one for the combinations, G23 = B 2 /B 3 where B 2 and B 3 are computed in the BMU basis.The discretisation effects are enhanced in the ratio, illustrating the fact the size of the cutoff effects depends on the choice of basis, see the discussion in the text.

) scheme at µ = 3
GeV. Results are obtained from the ratios R i and from the bag parameters B i .The systematic errors combine the chiral and the discretisation errors, the percentage error is obtained by adding statistical and systematic errors in quadrature.compatible within error and that the error remains the same.For the quark masses, we take advantage of the precise values quoted in[5], m d = 3.162(51) MeV, m s = 87.35(89)MeV for the SMOM-(γ µ , γ µ ) scheme, m d = 3.011(50) MeV, m s = 83.19(87)MeV in the SMOM-( / q, / q) scheme.and m s = 81.64(117)MeV and m d = 2.997(49) MeV in MS.Our results are reported in Tables XX,XXI andXXII.

FIG. 8 :
FIG. 8: Ratios R i for our two intermediate SMOM schemes renormalised at µ = 3 GeV, orange indicates positive correlation and blue anti-correlation, darker colours show a stronger correlation.These are a visualisation of the data from Table XXIV.
TABLE XXV: Correlation matrices for the bag parameters B i and the combinations G ij at µ = 3 GeV.
FIG.9: Correlation matrices in the (γ µ , γ µ ) scheme renormalised at µ = 3 GeV for the bag parameters B i and the combinations G ij .This is a visualisation of the data in TableXXV.

TABLE II :
Summary of our lattice ensembles.The heaviest mass of the coarse ensemble is not used in the

TABLE III :
Central values and error budget for our final results renormalised at µ = 3 GeV.Note that for our non-perturbatively renormalised results in the SMOM-(γ µ , γ µ ) and ( / q, / q) scheme, the error Total * does not include any perturbative uncertainty (PT).We also show the error budget for our MS results where only SMOM-schemes have been considered.The central value is obtained using SMOM-(γ µ , γ µ ) as intermediate scheme.For illustration, in the second part of the table, we give the error budget if we only use the RI-MOM scheme.See text for details.

TABLE IV :
Same as table III for our µ = 2 GeV results.

TABLE V :
Final results for R i , B i and G ij renormalised at µ = 3 GeV.The first error is statistical, the second one is an estimate of the systematic error (chiral and discretisation errors combined in quadrature).

TABLE VI :
[41]arison of the bag parameters B i at 3 GeV in the SUSY basis in the MS scheme of[41].

TABLE IX :
χ 2 /d.o.f of the global fits using a chiral fit (χPT) or a linear fit in m 2 P .

TABLE XIII :
Z/Z 2 V factors for the (27, 1) operators at 3 GeV on the fine lattice a = a 32 .

TABLE XVI :
Fit results for the ratio of bare matrix elements on the coarse ensembles.The corresponding correlation matrix can be found in the text.

TABLE XVII :
Same as XVII for the fine ensembles.

TABLE XVIII :
Correlation matrix for the coarse lattice with am s = 0.04 and am ud = 0.005.
observe that the covariance matrices are very similar for the two lattice spacings.We find almost 100% correlation between the colour partners (O 2 , O 3 ) and (O 4 , O 5 ).The remaining correlation coefficients are of order 60 − 70%.E.Matrix elements from Methods A and B

TABLE XXII :
Same as the previous tables but the results have been converted to MS at µ = 3 GeV.The third error is the estimate of the error due the perturbative matching and is kept separate from the other systematic errors.For the percentage error, all the errors have been added in quadrature.

TABLE XXIII :
Collection of results for the bag parameters using different methods.Results are given in (γ µ , γ µ ) scheme at 3 GeV and the errors combine statistical and systematic.

TABLE XXIV :
Correlation matrices for the ratios R i in our SMOM schemes at µ = 3 GeV.