Neutral kaon mixing beyond the Standard Model with nf=2+1 chiral fermions part II: Non Perturbative Renormalisation of the $\Delta F=2$ four-quark operators

We compute the renormalisation factors (Z-matrices) of the $\Delta F=2$ four-quark operators needed for Beyond the Standard Model (BSM) kaon mixing. We work with nf=2+1 flavours of Domain-Wall fermions whose chiral-flavour properties are essential to maintain a continuum-like mixing pattern. We introduce new RI-SMOM renormalisation schemes, which we argue are better behaved compared to the commonly-used corresponding RI-MOM one. We find that, once converted to MS, the Z-factors computed through these RI-SMOM schemes are in good agreement but differ significantly from the ones computed through the RI-MOM scheme. The RI-SMOM Z-factors presented here have been used to compute the BSM neutral kaon mixing matrix elements in the companion paper [1]. We argue that the renormalisation procedure is responsible for the discrepancies observed by different collaborations, we will investigate and elucidate the origin of these differences throughout this work.


I. INTRODUCTION
Numerical simulations of Quantum ChromoDynamics (QCD) allow for first-principle evaluations of hadronic matrix elements, which play a crucial rôle in theoretical calculations as they encapsulate the low-energy physics of a process.Computation of such matrix elements is usually done in two steps: Firstly, the bare quantities of interest are computed at finite lattice spacing a, whose inverse plays the rôle an ultra-violet regulator.Secondly, these quantities have to be renormalised in order to be divergence-free and have a well-defined continuum limit (a 2 → 0).There are two known nonperturbative methods to perform this renormalisation: the Schrödinger Functional (SF) scheme and the other being the Rome-Southampton method [2].We choose to work with the latter, for practical reasons (the interested reader can find a recent study of the same set of operators within the SF in [3] and [4]).In phenomenological applications the renormalised quantities are then matched to a scheme in which the corresponding short distance contributions can be computed, this is commonly performed in the modified minimal subtraction scheme MS, see for example [5,6].
Let us begin by considering the matrix element of an operator O which renormalises multiplicatively, and with O bare (a) being a bare matrix element computed at finite lattice spacing a.
We denote Z RI the corresponding renormalisation factor computed on the same lattice (following the Rome-Southampton method) in a regularisation independent (RI) scheme.The precise definition of the schemes (RI-MOM or a RI-SMOM) will be given in the next section.Within our conventions, at some renormalisation scale µ, the renormalised matrix element is given by which now has a well defined continuum limit O RI (µ, a) . ( Suppose now that this operator occurs in the determination of some physical quantity, say an amplitude.For example in a typical phenomenological application the hadronic matrix element has to be combined with a Wilson coefficient C(µ) computed in continuum perturbation theory (the hadronic matrix element describes the long-distance effetcts and the Wilson coefficient the short-distance ones).Both of these must be computed in a common scheme, MS, to be matched where R is the conversion factor from the RI scheme to MS. Eq. 3 can easily be generalised to the operator mixing case where O and C become vectors, and R and Z become matrices.We remind the reader that although the renormalisation is performed non-perturbatively, the matching to MS from the RI scheme (R MS←RI (µ)) has to be done using continuum perturbation theory as MS is not possible to implement on the lattice.
Accurate matching of lattice operators using the Rome-Southampton technique requires the matching scale µ (given by the magnitude of a momentum µ = p 2 ) to be well-separated from both the scales where non-perturbative effects of QCD such as chiral symmetry-breaking become important and the (inverse) lattice scale where cut-off effects dominate; ideally one would impose The first condition ensures that a perturbative treatment of the matching to MS is justified, while the latter ensures that the lattice artifacts are under control 1 .The MS renormalisation factors should be independent of the intermediate (RI) scheme used; however, in practice there will be some dependence due to systematic uncertainties in the lattice matching step as well as perturbative truncation errors in the continuum matching.
We compute the Z-matrix needed to renormalise the operators required for the determination of neutral kaon mixing beyond the Standard Model (BSM).As is usually done by the RBC-UKQCD collaboration, we implement momentum sources and partially-twisted boundary conditions.The use of momentum sources (introduced by QCDSF in [9]) results in very low statistical noise while the use of partially-twisted boundary conditions allows us to change µ = p 2 smoothly while keeping the orientation of p fixed [10][11][12].In this way we do not discontinuously 'jump' into different hypercubic representations as p 2 varies, resulting in Zs which are smooth functions of p 2 .
In principle, after extrapolation to the continuum and conversion to MS (or any common scheme) at a given scale, the results should be universal -up to truncation error of the perturbative series -and in particular should not depend on the details of the discretisation.The physical results |K / K|O SM K , in MS at 3 GeV in the SUSY basis.The statistical and systematic errors have been combined in quadrature.Although in principle these quantities should agree up to α 2 s errors, the RI-MOM results differ significantly from the (γ µ , γ µ ) and (/ q, / q) ones, which are consistent with each other.
The latter are RI-SMOM schemes whose precise definitions are given in this work.The N f = 2 + 1 results quoted here are obtained with exactly the same framework apart from the intermediate renormalisation scheme, see [1].We argue that the difference comes the renormalisation and we suggest to discard the results obtained with the RI-MOM scheme.Not included are results obtained with N f = 2 + 1 + 1 flavours by the ETM collaboration [13], which are roughly consistent with the N f = 2 RI-MOM results, and by SWME with N f = 2 + 1 [14,15], which are in a good agreement with our RI-SMOM results, see text for discussion.
could still depend on the number of dynamical flavours but, by experience, we do not expect this dependence to be important for the weak matrix elements discussed in this work.In the past few years, these matrix elements have been computed by three different collaborations and some discrepancy has been observed for two of the four relevant four-quark operators.The first results with dynamical quarks was reported by our collaboration in [16], it was done with N f = 2 + 1 flavours of dynamical quarks at a single value of the lattice spacing.Shortly after our work was published, the ETM collaboration published their study with N f = 2 flavours and several lattice spacings [17], they found compatible results (within 2σ for O 5 ).Then the SWME collaboration [14] reported on their computation, using N f = 2 + 1 flavours of improved staggered and again several lattice spacings.They find an important disagreement for two of these matrix elements.The ETM collaboration has then repeated their computation with N f = 2 + 1 + 1 flavours [13] and found roughly the same results as in their previous study (again only within ∼ 2σ for O 5 and the new result is now in perfect agreement with our old result).
In [1,18], we added another lattice spacing and investigated the origin of the discrepancy.In particular for the non-perturbative renormalisation procedure, in addition to the traditional RI-MOM scheme, we have implemented new intermediate renormalisation schemes, called (γ µ , γ µ ) and (/ q, / q) which satisfy the RI-SMOM condition, and therefore exhibit non-exceptional kinemat-ics.As summarised in Table I, we find that the results depend significantly on the intermediate renormalisation scheme: • If we use the traditional RI-MOM scheme with exceptional kinematics, we reproduce our old result and are compatible with ETMc, who used the same RI-MOM scheme .
• With the RI-SMOM schemes, our results for O 4 and O 5 are significantly different from our old RI-MOM results, but are consistent with each other.
• Our new RI-SMOM results are also in good agreement with SWME, who perform the renormalisation at one-loop in perturbation theory.This has been confirmed by the update of SWME [15].Therefore, one of our main conclusions in [1,18] is that the renormalisation procedure is the source of the discrepancy and we suggest to discard the results obtained with exceptional kinematics due to the systematic uncertainty in the pion pole subtraction.
In Table I, we choose to compare the results for R i [19] as they give directly the deviation of new physics with respect to the SM contribution.Since we could not find these quantities in [13][14][15], we do not show the results from ETMc (N f = 2 + 1 + 1) and SWME.However such a comparison for the bag parameters can be found in [1].
The main purpose of this work is the definition of RI-SMOM schemes for the BSM operators, generalising what has been done for the Standard Model B K and for K → ππ matrix elements [8,[20][21][22][23][24][25][26][27].These RI-SMOM schemes use non-exceptional kinematics with a symmetric point and have much better infrared behaviour, resulting in the suppression of pion pole contribution and wrong-chirality operator mixing [28,29].We argue in this work that at this point, results obtained using the RI-MOM scheme should be approached with skepticism or, if possible even discarded, at least for these quantities (the renormalisation of BSM kaon mixing operators).In addition, we define two new NPR schemes which both have different perturbative truncation systematics, upon comparing the two we can cleanly estimate the systematic from the renormalisation procedure The paper is organised as follows: in the next section we explain our procedure to obtain the Z factors.In Section III we give the explicit definitions of the projectors, which complete the definition of the schemes.The numerical results can be found in Section IV.In Section V we discuss the pole subtraction and the advantages of using the RI-SMOM schemes.Section VI contains our conclusions.Further details can be found in the appendices, where we give the relevant Z-factors for the bag parameters, the non-perturbative scale evolution of our renormalisation matrices, its comparison with perturbation theory, and finally the Fierz relations for the operators considered here.

II. METHODOLOGY
The Non-Perturbative-Renormalization (NPR) procedure works as follows: we compute numerically the Landau-gauge-fixed Green's functions of the operators of interest between incoming and outgoing quarks in a given kinematic configuration.After amputation of the external legs, projection onto the Dirac-colour structure and extrapolation to the chiral limit, we require that the renormalised Green's functions are equal to their tree-level values.Since we renormalise a set of four-quark operators which can mix, this renormalisation condition defines a matrix of renormalisation factors.We will discuss importance of the choice of kinematics; in particular the renormalisation condition is imposed for a certain momentum transfer p which defines the renormalisation scale µ = p 2 .For comparison we will also implement the original RI-MOM scheme [2], for which results at a single lattice spacing were presented in [16], but we chose to discard them for our final result in [1] as we will argue herein they appear to suffer from large systematic errors.
In the Standard Model only one operator contributes to neutral kaon mixing (a and b are colour indices) Beyond the Standard Model, under reasonable assumptions, four other four-quark operators are required (seven if parity is not conserved).Different choices of basis are possible but since we are concerned here with renormalisation, we find it convenient to only consider color-unmixed operators, i.e. those with the same colour structure as Q 1 .In Appendix VII D we give the relation between the colour-mixed and colour-unmixed operators.In order to simplify the equations, we do not explicitly write the colour indices, the contraction over spin and colour indices is simply indicated by the parentheses.We define the BSM operators (see for example [5]): where In practice we only consider the parity-even part of these operators, We will refer to Eq. 7 as the NPR basis (the relation between the SUSY and the NPR basis can be found in Appendix VII D).The factor 1/4 in Q 5 of Eq. ( 6) ensures that our definition matches the ususal lattice convention: These four-quark operators mix under renormalisation and -in a massless scheme -the mixing pattern is given by the chiral properties of these operators.They belong to three different repre- close to the continuum one.However the effects of spontaneous chiral symmetry breaking will be present at some level and could introduce some forbidden mixing (mixing between operators which belong to different representations of SU L (3) × SU R (3)).These unwanted infrared contaminations decrease as the renormalisation scale is increased beyond the typical interaction scale of QCD (Λ QCD ).
Such unphysical mixings are strongly suppressed in SMOM schemes (compared to the RI-MOM scheme), where the choice of kinematic prevents the contribution of exceptional momentum configurations [28].In practice, we will take the degree to which the expected continuum mixing pattern is satisfied as a quantitative indicator of the degree to which the NPR condition Eq (4) is satisfied.
The choice of kinematic for the RI-SMOM schemes is illustrated in Fig. (1).There are two different momenta p 1 and p 2 such that the momentum transfer is p2 = (p 2 − p 1 ) 2 .In this way a single renormalisation scale µ = p 2 is maintained and momentum flows through the vertex, which suppresses unwanted non-perturbative behaviour compared to the original RI-MOM scheme 2 .In practice we need two (momentum source) propagators, we associate a momentum to a given flavour, here p 1 for the d-quark and p 2 for the s-quark.The momenta used are of the form (the Euclidiantime component is the last coordinate) so that p = p 2 − p 1 = 2π L (m, m, 0, 0) .Since we use twisted boundary conditions in the valence sectors the momenta are not restricted to the Fourier modes.Our conventions are such that where θ is the twist angle of the boundary condition and n is an integer Fourier mode.
Our choice of convention is the following: with respect to the position of the vertex x, 1.An incoming s quark with momentum p 2 is denoted by 2. An outgoing d quark with momentum −p 1 is denoted by For each operator Q i of Eq.6 we compute the following Green's function (where we define xi = where the Greek letters denote combined spin-color indices.The color-Dirac structure of the four- (i) , (there is no summation over i in Eq. 13).For the numerical implementation, we have only considered four-quark operators that are color unmixed (the color partners can be obtained by Fierz transformation, see Appendix VII D).For example, if i, j, k, l are Dirac indices and a, b, c, d are color indices, then for the operator Q 2 , we have The vertex functions are then amputated where we have introduced the inverse of the "full momentum propagators" We still have to project these amputated vertex functions in order to obtain the renormalisation matrix.This is described in the next section.

B. Projection
Following Eq. 1, we introduce the renormalisation matrix Z which relates the renormalised four-quark operators to the bare ones (we drop the superscript RI for the Z factors) Denoting by Π bare i the bare amputated Green's function of the four quark operator Q i , the matrix Z ij is defined by imposing the renormalisation condition3 : where Z q is the quark wave function renormalisation.In the previous equation, P i projects onto the tree-level spin-colour structure of Q i : where the superscript (0) denotes the tree-level value.The fact that there is a non-vanishing momentum transfer in the vertex gives us more freedom for the choice of projectors.In this work, we introduce two different sets of projectors: P (γ µ ) and P (/ q) , they are defined explicitly below.We also need a prescription for the quark wave function Z q .This is done in two steps: first we cancel the factors of Z q in (18) using the vertex function of the local vector current.The value of Z V is then determined from some Ward identity in [30].We implement two projectors P (γµ) V and P V to obtain Z V /Z q .The choices of projectors for the four-quark operators and for the vector current define the non-perturbative scheme.Denoting by A and B the choices of projectors, ie (γ µ ) or (/ q), for both the four-quark operators and the vector current, the NPR condition for the scheme (A, B) reads The matrix Z (A,B) converts the bare four-quark operators onto the renormalised four-quark operators in the RI-SMOM scheme (A, B).
In [1] the primary quantities we presented were the ratios of particular BSM matrix elements over the SM one4 So we now consider the Z factors needed for these ratios.Introducing some notation for the projected vertex functions From Eq. ( 20), neglecting the mixing of the (27, 1) with the other operators, one finds that the quantity is independent of B, which is the choice of the projector for the denominator of Eq.( 20).Therefore, although in principle we have defined four RI-SMOM schemes (γ µ , γ µ ), (γ µ , / q), (/ q, γ µ ), (/ q, / q), in this work we mainly consider the "diagonal" schemes, for which A = B, namely (γ µ , γ µ ) and (/ q, / q).

III. NON-EXCEPTIONAL SCHEMES A. Choice of projectors
For the quark wave function renormalisation, we make use of two different definitions of Z q .
The factors Z q /Z V are determined by imposing the condition The two projectors we use are P (γ µ ) V and P (/ q) V , they are defined explicitly by: where Π V is the amputated Green's functions of the vector and axial-vector current.
The basis of the four-quark operators is given in Eq. ( 6), our convention is such that that all the operators are "colour-unmixed".The definition of the γ µ -projectors is straighforward: they are defined with the same spin-colour structure as their respective operators.Explicitly, for the SM operator we have For the / q schemes, following [31], we replace the γ µ matrices by / q/ q 2 , for example Similarly for the (8, 8) doublet we have For the / q projectors, in the case of P 2 , we apply the same recipe as the previous operator.For P 3 , we take advantage of the Fierz arrangements to "trade" the S and P Dirac matrices for the vector and axial ones.Explicitly we define Where the latter is now "colour-mixed" (this set of projector has already been introduced in [25,26] in the context of K → ππ decays).
Finally for the (6, 6) operators we define and where P R,L = 1 2 (1 ± γ 5 ).Imposing Eq. 18 with the projectors given above defines the various schemes (A, B) where A and B are either γ µ or / q

B. Tree-level values
For SM operator the tree-level vertex function reads: and equivalently for the other Dirac structures.The projectors act on the vertex functions by simply tracing over the Dirac and colour indices, explicitly the tree-level version of Eq. 18 is The corresponding tree-level matrices (N = 3 is the number of colours) are and IV. NUMERICAL RESULTS

A. Non-perturbative Z factors
The renormalisation is performed on the same ensembles as in [1], the parameters are summarised in Table II.We implement numerically Eq. 20 and obtain the Λ matrices (as defined in are shown for the SMOM − (γ µ , γ µ ) scheme on the 24 3 lattice.
Eq 22) at finite quark mass for the the list of momenta listed in Tables III and IV.The parameters for these ensembles are summarised in Table II.We perform a chiral extrapolation, invert the result and then interpolate to the desired scale of 3 GeV.Strictly speaking, there is mismatch from m sea s = m sea ud , however the quark mass dependence is dominated by the valence sector, the sea contribution plays very little rôle here.Furthermore, for the RI − SMOM schemes the light quark mass-dependence is very mild, practically invisible at our renormalisation scale even within our high statistical resolution, and so we consider any associated systematic to be negligible.
Due to the use of partially twisted boundary conditions, we can simulate momenta arbitrarily close to the targeted point, hence only a very small, well controlled, interpolation (performed with a quadratic Ansatz) is required.We illustrate these points in Fig. 2. The numerical results for the Z factors at 3 GeV are given in tables V,VI, VII, VIII,IX and X.
In principle we only need momenta close to the scale we wish to present our final results at (here µ = 3 GeV), however it is useful to compute the Z factors for a larger range, say between 2 and 3 GeV.We can then compare the non-perturbative scale evolution to its perturbative approximation and estimate the effects of truncating the perturbative series for the various schemes.Furthermore, since the running has a continuum limit, we also obtain a nice handle on the discretisation effects.

B. Conversion to MS
It is commonplace to convert the renormalised matrix elements computed on the lattice to the MS scheme.In that way, the Wilson coefficients can be combined with the matrix elements to produce phenomenological predictions.The conversion from the RI − MOM or RI − SMOM to MS is done in continuum perturbation theory.The matching coefficients are known at the one-loop level for RI − MOM from [32] and [6].The situation is different for the RI − SMOM schemes: the relevant conversion factors of the (27, 1) operator have been computed in [31].The conversion matrix for the (8,8) operators can be extracted from [33] where the conversion was computed for the ∆S = 1 K → ππ four-quark operators.For the (6,6) operators, the coefficients were unknown and have been computed for this work.The full expression can be found in Appendix VII B.
To obtain α s at µ = 3 GeV in the three-flavour theory, we start from α s (M Z ) = 0.1185 (6), we use the four-loop running given in [34,35] to compute the scale evolution down to the appropriate charm scale, while changing the number of flavours when crossing a threshold, and then run back up to 3 GeV in the three-flavour theory.
The values of the one-loop conversion matrices and the Z factors in MS (ie the Z factors which convert our bare matrix elements to MS) are given in tables V,VI, VII, VIII,IX and X.For completeness, we also give the conversion factor for the original RI − MOM scheme (the equivalent of the second columns of the above-listed tables) The conversion to MS is then given by Z MS = R MS←(scheme) × Z (scheme) where (scheme) can be RI − MOM, (γ µ , γ µ ), (/ q, γ µ ),(γ µ , / q) or (/ q, / q).
We observe that in general, the "diagonal" schemes (γ µ , γ µ ) and (/ q, / q) have a better perturbative convergence than the off-diagonal ones.At 3 GeV, the conversion matrices are rather close to the identity (which probably explains why our results agree so well with SWME).For our two favorite schemes, we find that after conversion to MS, the numbers agree rather well.The convergence of the perturbative series and the effects of the lattice artefacts could also be estimated by looking at the step-scaling matrices, which we do in the next section (see also Appendix VII C).

C. Non-perturbative scale evolution and comparison with perturbation theory
The scale evolution matrix, σ(µ 1 , µ 2 ) is a rich source of information, in particular it helps us to estimate the systematic errors affecting the renormalisation procedure.We define where Z is the 5 × 5 matrix defined in Eq. 20. (Although in practice we take the chiral limit of the right hand side of Eq. 37, once again in order to simplify the notation, we discard any reference to the quark masses.) The scale evolution matrix has a universal continuum limit and may be directly compared to continuum perturbation theory.The continuum extrapolation is performed assuming a linear behaviour in a 2 .For this step the use of twisted boundary conditions is essential, since it allows us to vary µ continuously holding the momentum orientation (and O(a 2 ) coefficients) fixed.
The continuum extrapolation of σ ii (2 GeV, µ), where 2 GeV ≤ µ ≤ 3 GeV, is shown in Figs.7-9, compared with continuum perturbation theory.We find in general good agreement with the perturbative series, indicating that the a 2 extrapolation is valid and discretisation effects are under control.An example of off-diagonal matrix elements can be found in Fig 10.
By comparing the non-perturbative running to its perturbative approximation, we can estimate the quality of the perturbative series for the various schemes.This is important in view of the perturbative macthing of the NPR factors to MS.In order to compare the scale evolution matrix to the perturbative estimates, it is useful to construct the quantity σ(µ 1 , µ 2 )σ −1 PT (µ 1 , µ 2 ), which is equal to 1 5×5 up to higher-order terms not included in the perturbative expansions, residual discretisation effects, and non-perturbative contributions.These quantities are shown in Figs.14.When running from 3 to 2 GeV, we find that these effects are typically of order a few percents, and in many instances much less.

V. RI-MOM RENORMALISATION SCHEME
In addition to the RI-SMOM renormalisation schemes used to obtain our main results [1], we also implemented RI-MOM renormalisation conditions for the intermediate scheme.The RI-MOM scheme differs in the kinematic configuration of the vertex functions, which depend on a single momentum vector (obtained by setting p 1 = p 2 in Eq. ( 13)).Vertex functions in this "exceptional" configuration can have large contributions from infrared poles which go as inverse powers of the quark mass (m 2 π ) and momenta; as our renormalisation matrices are defined in the chiral limit m → 0 (here and throughout this section m = m bare + m res ) we have an unphysical divergence due to this scheme, which must be subtracted.These pole contributions are suppressed by powers of p 2 but in practice turn out to be large for momenta accessible in our Rome Southampton window.
As the m → 0 limit is approached the raw RI-MOM data clearly suffers from pole contamination, the effect of these pion poles is clearly visible in our data, in particular in the Λ i3 and Λ i4 elements (Fig. 3); in contrast the RI-SMOM data have only a weak mass dependence and tend to Z −1 in the m → 0 limit (Fig 4).
At large µ = p 2 the matrix of vertex functions Λ will become block diagonal in the chiral limit if the effects of spontaneous chiral symmetry breaking are suppressed.In the RI-MOM scheme chiral symmetry breaking effects can be extremely enhanced in the m → 0 limit; as a result the chiral structure is strongly broken.This can be seen for example in Fig. 5

(right).
We focus first on the chiral extrapolation and work at fixed momentum.In order to extract Z ij from the RI-MOM data, we fit the mass dependence of the vertex functions Λ ij .In principle we expect the vertex function to exhibit poles which go like 1/m and 1/m 2 (see for example [17]), and so will be described by the general form First, we observe that not all the matrix elements require a pole subtraction.In that case, we just perform a linear fit in the quark mass (ie B = C = 0) with the three (lightest) unitary 5 .quark masses: am bare = 0.005, 0.01, 0.02 on the 24 3 and 0.004, 0.006, 0.008 on the 32 ij .This gives results equivalent to fitting Λ to the form A + B/m.We observed that amΛ is consistent with linear am behavior to justify neglecting the (am) 2 term, and also found the data after subtracting the pole contribution is linear.After subtracting the pole we find good restoration of the chiral block structure for the 32 3 ensemble (Table XV), The chiral restoration is not as good on the 24 3 ensemble, the residual matrix elements are of the order of a few %.However we observed that they do affect the physical matrix elements, and that different fit procedure give the same residual, see below and Table XIV.
Since this infrared contamination completely dominates some of the raw data in the RI-MOM scheme, we investigated the effect of this pole subtraction, in particular we want to have a reasonable estimate of the systematic error associated with the procedure.On the 24 3 , we used another ensemble m val light = m sea light = 0.03 and have implemented different fit forms.We fit each of the Λ ij with the forms A + Dm (fit-0), A + B/m + Dm (fit-1), and A + C/m 2 + Dm (fit-2).We find that in cases with significant singular behavior, the fit-1 has χ 2 < 1 and fit-2 has χ 2 1.For j = 1, 2, 5 there is no evidence of 1/m behavior and the results are compatible with fit-0.The fits are shown in Fig. 3 for the chirally allowed elements Λ 23 , Λ 33 , Λ 44 , and Λ 54 .From this we conclude that any 1/m 2 dependence is to a large degree suppressed in the range of m val for which we have data, and we determine Z −1 ij assuming the form of fit-1.As a check on the procedure we compare the fit-1 results on the 24 3 ensemble to the linear fit procedure.For the linear fit method we threw out the heaviest (am = 0.03) mass point because we found a degradation in the χ 2 (though central values remain consistent), and which we attribute to neglecting the quadratic term.The results from the two subtraction methods are shown in Fig. 3 and Table XIV.There is a slight tension in the extrapolated results which highlights that some uncontrolled systematic due to specifics of the subtraction procedure may remain.
We also implemented Bayesian fits using the lsqfit package6 to include additional terms from Eq. ( 39) without requiring the number of data points to exceed the number of fit parameters.
Table XIV compares results of frequentist and Bayesian fitting on the 24 3 ensemble, for both chirally allowed and forbidden elements.The Bayesian fit of the full form (39) is consistent with the results of the other methods but with larger uncertainties.For the chirally-forbidden elements, the single pole fits (fit-1 and Lin.fit) find values which differ significantly from zero, whereas the Bayesian method finds best-fit values very close to zero, but with errors comparable to the size of the central values in the single pole case.As another consistency check of the method, we should also find an approximate recovery of the block diagonal structure expected from chiral symmetry after removing the singular parts of the data.Although to a decent approximation the terms that are chirally-forbidden are suppressed after the pole subtraction, we find that the values are statistically non-zero and the magnitude of chirally-forbidden elements tend to be larger for the pole-subtracted (Λ i,3/4 ) compared to elements that do not require pole subtraction (Λ i,1/2/4 ).Fig. 5 shows the mass and µ dependence of chirallyforbidden RI-MOM vertex functions for a case without discernible singular structure (Λ FIG. 4. Same as Fig. 3, from left to right: Λ 23 , Λ 33 (first row) and Λ 44 , Λ 54 (second row), for the nonexceptional (γ µ , γ µ ) scheme.Here we fix the momentum µ close to 3 GeV..In that case we observed a very mild, linear, quark mass dependence.In contrast to the RI-MOM case, no pole subtraction is required (we show the vertex function without applying any pole strubaction procedure).
and where the pole behavior is clearly visible (Λ 24 , right).These results should be contrasted with the RI-SMOM results shown in Fig. 6, where in all cases the chirally-forbidden elements extrapolate very nearly to zero.
On the 32 3 ensemble we also compare results of including the single pole or both poles using a Bayesian fit, and results from the linear fit method, shown in Table XV.The results again agree with the linear fit results but have larger associated uncertainties.Note here the chirally-forbidden elements obtained from the linear fit method are much smaller than in the 24 3 case and are in fact zero within errors.We also tried including the 1/m terms in 'global' fits by constraining the 1/m coefficient coefficient in Λ i3 to be the negative of the coefficient in Λ i4 , which we observed to be the case.Althought this strategy seems to improve somewhat the fit quality, the numerical resuls  were essentially unchanged.
We choose two options when we compute our renormalisation matrices: firstly we invert the whole matrix of fit parameters for all Λ ij and secondly we invert only the block diagonal elements of the matrix, zeroing by hand the chirally-forbidden elements.We will label these as the Not Block-Diagonal (NBD) and the Block Diagonal (BD).
Here we list the results for Z-matrices obtained in RI-MOM from the linear fit method.On the And the 32 3 at µ = 3.01 GeV, It is evident that the BD and NBD Z-matrices are not too dissimilar, we take the difference in results of the operators renormalised using either of these as a systematic for our final RI-MOM results.Our results in the RI-MOM scheme after chiral extrapolation and interpolation to µ = 3 GeV read In conclusion, the infrared contamination in some of the RI-MOM vertex functions makes it difficult to extract the Z-factors precisely in the m → 0 limit, where these contributions diverge.
These effects also strongly break the chiral structure one expects to recover for µ Λ QCD , though this structure is restored (albeit imperfectly) after subtraction of the pole contributions.For these reasons, we find that the RI-MOM scheme (with exceptional kinematics) suffer from systematic errors which are difficult estimate.Applying different strategies to subtract the poles, we find that final results vary by 5% in the worse case.
In contrast the SMOM procedure strongly suppresses these infrared effects -evidence of chiral symmetry breaking disappears in the am → 0 limit at sufficiently large µ (Fig. 6), and the chirally-allowed Z-factors have very mild linear mass dependence.We also note that the SMOM to MS matching factors are much closer to unity, suggesting a better behaved perturbative series and a reduced perturbative matching uncertainty.Therefore we strongly advocate using SMOM renormalisation conditions, which are theoretically much cleaner.
We have argued that the discrepancies from results of [13,16,17] come the renormalisation procedure.Because these discrepancies appear in those matrix elements affected by these issues, we suggest avoiding the RI-MOM renormalisation conditions, at least for this set of operators. 7ven we assess a rather conservative 5% systematic error from the renormalisation procedure in RI-MOM, our results are still not compatible with the RI-SMOM ones.It remains the possibility of a conspiracy between these infrared artefacts and omitted term in the perturbative matching.
Even if the latter should be of order α 2 s , the anomalous dimensions of those operators are rather large.Since a computation at the next order is technically very challenging, this systematic error is difficult to control without using multiple schemes.

VI. CONCLUSIONS
In this work we have defined and investigated new RI-SMOM intermediate schemes for the renormalisation of ∆F = 2 four-quark operators needed for neutral kaon mixing beyond the standard model studies.These schemes can easily be generalised to other processes.We have implemented these different schemes and shown that they lead to consistent results after continuum extrapolation and conversion to MS.These results are, however, inconsistent with those obtained using the intermediate RI-MOM scheme.
Although the theoretical advantages of the RI-SMOM schemes -as compared to RI-MOMhave been known for a long time, we have provided further numerical evidences in the context BSM kaon mixing: • No pole subtraction is required.
• The chirally-forbidden matrix elements are largely suppressed.
• The Z and conversion matrices are closer to the identity matrix; the scale-evolution between 2 and 3 GeV is relatively close to the perturbative prediction (known at next-to-leading order).
On the other hand, in the RI-MOM scheme the effects of chiral symmetry breaking can be large even at large momentum, and a procedure must be used to remove infrared contributions that dominate some vertex functions in the chiral limit.We investigated the effect of different subtraction procedures in our RI-MOM data and found some dependence on the procedure, which may be at least partly responsible for the discrepancies in O 4 and O 5 .These effects are particularly important in the (S + P ) and (S − P ) channels.We have shown that the RI-SMOM procedure is superior because the unwanted infrared behaviour is nearly completely suppressed (and has better pertubative behaviour).
Our study indicates these discrepancies in O 4 and O 5 could be due to a conspiracy of systematic errors in the RI-MOM scheme, the dominant ones being the infrared contamination and the truncation error of the perturbative series in the matching to MS (as these operators have rather large anomalous dimension).
In other to have a better control on the physical point extrapolation, we are currently investigating the effects of including physical pion-mass ensembles and a finer lattice spacing.Our preliminary analysis [38,39] shows that our results are stable and we hope to decrease the un-certainties on the BSM matrix elements by at least a factor of two.We are also investigating a strategy to run through the charm threshold with n f = 2 + 1 + 1 flavours [40,41].

B. Matching factors between the RI-SMOM schemes and MS
The conversion between the RI-SMOM schemes and MS (of [6]) is given at one-loop order.We define (we chose a negative sign for historical reasons) In the following expressions, the constant , where ψ is the PolyGamma function, N is the number of colors and ξ the usual gauge parameter (the non-perturbative Zfactors have been computed in the Landau gauge, ξ = 0).Note that the coefficients for the (27,1) and the (8, 8) operators were already known or could be derived from [31,33].
For the (8,8) doublet: ∆r and for the (6, 6) doublet: ∆r ) operator for the various schemes; left: (γ µ , γ µ ), right: (/ q, / q).We show the non-perturbative running computed on the coarse lattice, on the fine lattice and extrapolated to the continuum.We also compare with the perturbative prediction at leading-order (LO) and next-to-leadingorder (NLO).
non-perturbative scale evolution is qualitatively well-described by the Next-to-Leading perturbative prediction.In the worse cases we observe a deviation of around 5 % at 2 GeV.In a future we will include a finer lattice spacing to have a better handle on the discretisation effects.


(γ µ , γ µ ) 3 .The chirally-allowed elements which suffer from those pole contaminations are Λ 23 , Λ 33 , Λ 44 and Λ 55 and the chirally-forbidden are Λ 24 , Λ 34 , Λ 44 and Λ 53 .In this case, our main results are obtained on the same data with a single pole fit Ansatz C = D = 0. 5 Strictly speaking the setup is unitary in the light quark sector m val light = m sea light , but partially quenched for the strange as m val s = m val light = m sea s , however we have checked this effect is negligible within our systematic errors RI-MOM vertex functions.Here we multiply the data by am and fit amΛ ij = (am)Z −1 ij + B ij + O((am) 2 ) to a straight line to determine Z −1

FIG. 3 .
FIG.3.Chirally-allowed RI-MOM vertex functions with singular behaviour from the 24 3 ensemble.The result of fitting the raw data (circles) to fit-1 (dotted line) and a fit to the lightest three points with the form a + b/m (solid line), along with the result of subtracting the single pole contribution from each of the fits (same line type as respective fits through data).Quantities shown from left to right are Λ 23 , Λ 33 (first row) and Λ 44 , Λ 54 (second row) at fixed momentum close to 3 GeV.
FIG. 5. Left: Example of an amputated and projected Green function in the exceptional RI-MOM scheme at finite quark mass (on 32 3 ensemble) for different momenta.This specific quantity should vanish if chiral symmetry is exact.Right: Example of a RI-MOM vertex function with strong singular behavior.This specific quantity should also vanish if chiral symmetry is exact but is affected by large infrared contaminations.

FIG. 10 .
FIG.10.Same as the previous plot for the scale evolution of the non-diagonal(8,8) mixing matrix element σ 32 and σ 33 .

TABLE I .
Example of results for the ratio of the BSM matrix elements over the SM one R i = − p 2 ) 2 .This configuration prevents the existence of a channel with zero-momentum transfer.

TABLE II .
Summary of the lattice ensemble used in this work.Since the renormalisation is performed with momentum sources, only a few configurations are needed (between ten and twenty for each ensemble).

TABLE III .
List of momenta for the 243lattices.Here we fix the Fourier mode to n = 3 and only change the twist angle θ, see Eq.10.

TABLE IV .
List of momenta for the 32 3 lattices.
FIG. 2. Example of amputated and projected vertex functions at the simulated momenta and quark masses(left) and interpolation of a Z matrix element to the 3 GeV-scale after chiral extrapolation (right).Results

TABLE V .
Z/Z 2 V factors for the (27, 1) operators at 3 GeV for a = a 24 .
TABLE XIII.Same for the running matrix of the (6, 6) operators.

TABLE XIV .
Comparison of fit results on 24 3 using "linear fit method", frequentist fit with 1/m term (1/m 2 term set to zero), Bayesian fit with only 1/m term (* result uses only lightest three masses), and Bayesian fit with both 1/m and 1/m 2 terms.The lower set of values corresponds to chirally-forbidden elements.
The others are new, they have been computed for this work.First we have the matching factors for the (γ µ , γ µ ) scheme, (6,6) 2 ] Same as the previous plot for the scale evolution of the diagonal(6,6)mixing matrix element σ 44 and σ 55 .