BK-parameter from Nf = 2 twisted mass lattice QCD

We present an unquenched $N_f=2$ lattice computation of the $B_{K}$ parameter which controls $K^0-\bar K^0$ oscillations. A partially quenched setup is employed with two maximally twisted dynamical (sea) light Wilson quarks, and valence quarks of both the maximally twisted and the Osterwalder--Seiler variety. Suitable combinations of these two kinds of valence quarks lead to a lattice definition of the $B_{K}$ parameter which is both multiplicatively renormalizable and O($a$) improved. Employing the non-perturbative RI-MOM scheme, in the continuum limit and at the physical value of the pion mass we get $B^{\rm RGI}_K=0.729\pm 0.030$, a number well in line with the existing quenched and unquenched determinations.


Introduction
Strong interaction effects in neutral K 0 −K 0 meson oscillations are parametrized by the so-called "bag parameter" B K [1]. Several lattice computations of B K based on different lattice fermion discretizations have been performed in the years, mostly in the quenched approximation. Attention has now shifted to unquenched estimates. See refs. [2,3] for an updated review on the subject. A compilation of recent data is provided in Fig. 1 where renormalization group invariant (RGI) values of B K (denoted as B RGI  [4][5][6][7][8][9][10][11][12], respectively. The symbol ⋆ denotes determinations where the continuum extrapolation was carried out.
The standard way of computing B K with Wilson fermions yields limited accuracy due to partial control of two sources of systematic effects.
• First of all, the ∆S = 2 four-fermion operator relevant for the B K calculation, namely mixes under renormalization with four other operators of the same dimension, but belonging to different chiral representations [13] (the so-called "wrong chirality mixing" phenomenon). The reason behind this feature is the lack of chiral symmetry of the regularized lattice action. In other regularizations (e.g. staggered or Ginsparg-Wilson fermions) partial or full chiral symmetry ensures that the operator (1.1) is multiplicatively renormalizable (i.e. does not mix with other operators). An important detail here is that mixing actually only affects the parity-even component to have maximally twisted sea quarks in combination with Osterwalder-Seiler (OS) valence quarks. We recall that any OS fermion can be thought of as a component of a maximally twisted doublet with twist angle +π/2 or −π/2, depending on the value of its Wilson r-parameter.
Then the relevant ∆S = 2 four-fermion operator is defined in terms of four distinct (maximally twisted) valence flavours, three with twist angle, say, +π/2, and the fourth with twist angle −π/2. The discrete symmetries of the valence quark action combined with the structure of the resulting four-fermion operator ensure that the latter is multiplicatively renormalizable, while its K 0 −K 0 matrix element remains automatically O(a) improved [22]. The same is then true for B K .
In the setup of ref. [22] unitarity violations occur as a consequence of the different regularization of sea and valence quarks. However, if renormalized masses of sea and valence fermions are matched, one can prove [22] that unitarity violations are mere O(a 2 ) effects. Moreover, it just happens that the valence quark content of the K 0 -meson consists of a strange/down pair with, say, the same twist angle, while thē K 0 -meson has necessarily a quark-anti-quark pair with opposite twist angles. Thus the two pseudoscalar mesons have masses that at non-zero lattice spacing differ by (numerically important) O(a 2 ) effects, which mainly come from the lattice artifacts (non-vanishing in the chiral limit) of the K 0 -meson mass. TheK 0 -meson mass on the other hand exhibits only O(a 2 (µ ℓ + µ s )) discretization errors (here µ ℓ/s denotes the mass of the light/strange quark). This effect has been recently studied in the quenched approximation adding the clover term in ref. [7], where it is numerically demonstrated that the scaling behaviour of both B K and the decay constants of the K 0 andK 0 mesons is only weakly affected (see also below). The situation about the O(a 2 ) mass-splitting between K 0 andK 0 is similar to the one already met in the unitary Mtm-LQCD framework, where a large O(a 2 ) artefact is only observed in the mass of the neutral pseudoscalar meson as it is not associated to a conserved lattice axial current [23,24].
In this paper we present an unquenched computation of B K based on the strategy of ref. [22] which exploits the lattice data coming from the state-of-the-art tm-LQCD simulations carried out by the ETM Collaboration (ETMC) with N f = 2 dynamical (sea quark) flavours, at three lattice spacings and "light" pseudoscalar meson masses in the range 280 MeV < m PS < 550 MeV. Renormalization is carried out non-perturbatively in the RI-MOM scheme [25]. In the continuum limit and at the physical value of the pion mass we get B RGI K = 0.729 (25) (17) ⇔ 0.729 (30) @ N f = 2 , (1.2) where in the first expression the two errors quoted are the statistical one (0.025) and (our estimate of) the systematic uncertainty (0.017, resulting from the quadratic sum of 0.004 from renormalization, 0.009 from control/removal of cutoff effects and 0.014 from extrapolation/interpolation to the physical quark mass point). In the second expression the total error is obtained by a sum in quadrature of the statistical and systematic ones. The error budget is further discussed in sect. 3.2. The specification @ N f = 2 is to remind that the whole computation of B RGI K has been carried out in QCD with two light (u and d) dynamical quark flavours. The quality of this result is extremely satisfactory and fully comparable with other quenched and unquenched numbers (see Fig. 1 for a compilation of recent determinations).
The outline of the paper is as follows. In sect. 2 we briefly recall the theoretical framework which allows to obtain an O(a) improved and multiplicatively renormalizable lattice expression of the K 0 −K 0 matrix element of the four-fermion ∆S = 2 operator (1.1). In sect. 3 we give the details of our numerical analysis. Conclusions can be found in sect. 4. We defer to a couple of Appendices a few more technical considerations. Preliminary reports about the present work have already appeared in refs. [3,26,27] 2 Twisted and Osterwalder-Seiler valence quarks at maximal twist The lattice setup of our study is that of ref. [22]. The regularization of sea quarks is different from that of valence quarks. The former is a standard tm-LQCD regularization with a doublet of degenerate, maximally twisted Wilson fermions (representing up and down flavours). In the so-called "physical" basis, the Mtm-LQCD fermionic action takes the form (2.1) where the Wilson's r parameter has been set to unity, ψ(x) is the quark flavour doublet, ∇ µ and ∇ * µ are nearest-neighbour forward and backward lattice covariant derivatives, µ sea is the (twisted) sea quark mass and M cr the critical mass. The lattice pure gauge action is the tree-level improved action of ref. [28], routinely adopted by the ETM Collaboration [29] in studies with two light sea quark flavours. Valence quarks are regularized as Osterwalder-Seiler (OS) fermions [30]. The action of each OS valence flavour q f reads where aM cr is the same number as in eq. (2.1). The critical mass M cr is nonperturbatively tuned to its optimal value [31] as described in ref. [29]. This ensures O(a)-improvement of physical observables and control of the leading chirally enhanced cutoff effects 1 .

Lattice operators
For the computation of interest, it is convenient to introduce four OS valence quark flavours q 1 , q 2 , q 3 , q 4 . Eventually q 1 and q 3 will be identified with the strange quark by setting µ 1 = µ 3 ≡ µ h , and q 2 and q 4 with the down quark by setting The key point is that the four-fermion operator is multiplicatively renormalizable once the Wilson parameters appearing in the action of the OS fermions are taken to satisfy the relations The factor 2 in the r.h.s. of eq. (2.3) guarantees that the correlator (2.11), where Q V V +AA is inserted between the interpolating K 0 =q 2 q 1 andK 0 =q 3 q 4 fields, is normalized as in continuum QCD (as one can check by direct application of the Wick theorem). The proof of the multiplicative renormalizability of the operator Q V V +AA is based on the discrete symmetries of the OS action and is provided in full detail in ref. [22]. A simpler proof is readily obtained in the so-called "twisted basis". In fact, upon performing the chiral field rotation from the "physical" (q f ) to the "twisted" (q tm f ) basis, the valence quark action in the massless limit (all µ f 's set to zero) takes the standard (untwisted) Wilson form, while the parity-even operator Q V V +AA gets rotated onto the parity-odd Q V A+AV . The latter, being expressed in terms of quark fields with standard Wilson action, is known [14,15] to be multiplicatively renormalizable. Moreover, the argument implies that the renormalized operator is given by the formula where, for consistency with literature on Wilson fermions, the renormalization constant is named after the form (V A + AV ) the operator takes in the "twisted" basis. In actual simulations we fixed the common value in eq. (2.4) by setting r 1 = 1. In order to compute B K , we also need to consider the axial currents which, as a consequence of the relations between the Wilson parameters r 1 , . . . , r 4 specified in eq. (2.4), when passing from the "physical" to the "twisted" basis take has no practical consequences in the evaluation of correlation functions of operators where no ghost fields appear.
the form of axial and vector operators, respectively. Consequently they renormalize according to the relations

Lattice correlation functions
Next we proceed to construct the correlation functions, involving the operators that are necessary to evaluate B K . The lattice is taken to be of size L 3 · T , with periodic boundary conditions on all fields and in all directions, except for quark fields in the time direction which satisfy antiperiodic boundary conditions. At a reference time slice y 0 , we define a "K-meson wall" with pseudoscalar quantum numbers andq 2 and q 1 quark fields, namely A second "K-meson wall", W 43 (y 0 + T /2), withq 4 and q 3 quark fields, is placed at the time slice y 0 + T /2. The four-fermion operator, Q V V +AA , is inserted at position x = ( x, x 0 ), with x 0 in the range y 0 ≪ x 0 ≪ y 0 + T /2. In this way lattice estimators of the bare B K -parameter can be calculated at several values of x 0 from the ratio involving the correlation functions 2 The signal-to-noise ratio is greatly enhanced by summing over the spatial position of the quark and antiquark fields in each "K-meson wall" as well as summing over that of the four-fermion operator. The x 0 -behaviour of the B K -estimator excludes significant contaminations from excited states (see also sect. 3.1). On the other hand in order to reduce the cost of computing the quark propagators we chose to employ a stochastic method. This consists in computing quark propagators by inverting the lattice Dirac operator for the relevant OS valence quark flavours on random sources of the form with free indices α (spin), x 0 (time), i (colour), x (space), while γ and z 0 are fixed as indicated. The Z 2 -valued η vectors carry only colour and three-space indices (colour and three-space "dilution") and are normalized according to Given the pattern of Wilson r-parameters specified in eq. (2.4) and the invariance under γ 5 -hermitian conjugation combined with r f → −r f of the OS lattice Dirac operator, it is enough to compute for each gauge configuration the propagator of the fields q 1 andq 2 on the random source at y 0 , and the propagator of the fields q 3 and q 4 on the random source at y 0 + T /2. The conditions (2.15) guarantee that unbiased estimators of the correlators C 3 , C 2 and C ′ 2 are obtained 3 . We found that, for an ensemble of a few hundred gauge configurations, employing one set of Z 2 -valued random sources (namely those in eq. (2.14)) per gauge configuration is sufficient to achieve a good statistical precision for our B K -estimators.

Lattice estimator of (renormalized) B K
As we said, the masses of (up/down) sea quarks and the masses [µ 2 ] R , [µ 4 ] R of (down) valence quarks must be tuned to the same (physical) value. At the same time the valence mass parameters [µ 1 ] R and [µ 3 ] R should be adjusted to the strange quark mass in order for the external states to be identified with K-particles 4 . In this way the renormalized ratioR( valence flavours in the partially quenched set up specified above, will differ from its continuum limit by only O(a 2 ) discretization errors. One peculiar feature of our approach which we have to keep in mind in the analysis of simulation data is the fact that kaon and anti-kaon masses are not degenerate as the two mesons involve differently regularized valence quarks. In our setting the kaon is made up of the valence quarks q 2 → q d and q 1 → q s with r 1 = r 2 , while the anti-kaon consists of the valence quarks q 4 → q d and q 3 → q s with r 3 = −r 4 . The K 0 −K 0 mass difference M 12 − M 34 → M K − MK is an O(a 2 ) effect which has an origin similar to that of the well known mass splitting one encounters in the pion sector [24]. Naturally the relative effect is much less important in the kaon case as the kaon mass is substantially larger than the pion mass 5 .
Taking this mass difference into account, we write for the three-point correlation function the large time expansion (y 0 ≪ x 0 ≪ y 0 + T /2) In this expression |P ij denotes a zero three-momentum pseudoscalar state with quark content q i andq j (i, j = 1, · · · , 4) and mass M ij . Similarly we have Furthermore, defining the pseudoscalar decay constants F ij through the equations one finds that they are not equal, but again they differ by O(a 2 ) effects. The renormalized lattice B K parameter which we will have to extrapolate to the continuum limit is finally obtained by averaging over the x 0 -plateau (which lies between y 0 and y 0 + T /2) of the estimatorR(x 0 ), viz.

Scaling checks
In this section we present some tests aimed at checking the size of the discretization errors affecting the quantities that enter our computation of B K (i.e. those that appear in eq. (2.20)). In order to keep the issue of lattice artifacts separated from those related to the extrapolation to "physical quark mass point", the present scaling tests have been performed at fixed reference values of the renormalized "light" and "heavy" (would be strange) quark masses in the MS-scheme at the scale of 2 GeV. For convenience we takeμ * ℓ ∼ 40 MeV andμ * h ∼ 90 MeV.    ) and ∆ F (panel (d)) as functions of (af 0 ) 2 at the reference quark massesμ * ℓ ∼ 40 MeV andμ * h ∼ 90 MeV (in the MS scheme at 2 GeV). We display the best linear fit in (af 0 ) 2 and the corresponding continuum limit result.
In Fig. 2, panels (a) and (b), we show the scaling behaviour of (M 34 /f 0 ) 2 and F 34 /f 0 , i.e. the mass squared and the decay constant of the "kaon" made up of valence quarksq 4 and q 3 , normalized by the chiral limit pion decay constant f 0 ≃ 121 MeV determined in ref. [35]. Owing to the choice r 3 = −r 4 the quantities M 34 and f 34 are our lattice estimators of the "kaon" mass and decay constant that should exhibit best scaling properties. The plots in panels (a) and (b) fully confirm this expectation, showing a nice a 2 -scaling with cutoff effects for the lattice "kaon" mass and decay constant estimates that are at the level of few percents (from the largest lattice spacing to the continuum limit).
In Fig. 2, panels (c) and (d), we show instead the scaling behaviour of the difference of the "kaon" mass squared and decay constant lattice estimators entering our B K computation, namely The quantities (2.21) are by construction mere O(a 2 ) lattice artifacts. We see from the panels (c) and (d) that both ∆ M and ∆ F extrapolate nicely to zero in the continuum limit, as expected. The large values of ∆ M in Fig. 2(c) signal the presence of a substantial cutoff effect on M 12 , with M 12 on the coarsest lattice differing by about 30% from its continuum limit value (∼ 570 MeV). On the other hand, the lattice artifact difference between F 12 at the largest lattice spacing and its continuum limit value is only about 5%. This is readily inferred by comparing Fig. 2(b) and Fig. 2(d).
Having taken r 1 = r 2 , the large and positive lattice artifact observed in M 12 was expected since for OS valence quark doublets, (unlike the case of twisted-mass quark pairs) no isospin component of the lattice isotriplet axial current is conserved in the limit of vanishing quark mass. Relying on arguments similar to those given in ref. [24], one can however conjecture that the large discretization error in M 12 is due to dynamically large matrix elements entering the Symanzik description of the lattice OS pseudoscalar meson mass, rather than to large coefficients in front of some terms of the Symanzik local effective action for OS valence fermions 6 . This conjecture suggests that the large cutoff effect detected in the lattice OS pseudoscalar meson mass is peculiar to this quantity (or to others directly related to it) but not a general feature of physical observables built in terms of OS valence quarks. The small cutoff effects observed in F 12 is well in line with such an expectation.
In Fig. 3 we show the scaling behaviour of B RGI K,lat (μ * ℓ ,μ * h ) evaluated according to eq. (2.20) using two different methods (referred to as M1 and M2) for the evaluation of the renormalization constants Z RGI V A+AV and Z A . The two methods (discussed in sect. 3.2.1 and Appendix A) differ only in the way one deals with O(a 2 ) artifacts and other unwanted systematic effects pertaining to the RI-MOM scheme computation of renormalization constants. The somewhat different slope in a 2 of the data-points in Fig. 3 should be ascribed to these effects. Nevertheless in both cases the a 2 -scaling behaviour of B RGI K,lat (μ * ℓ ,μ * h ) is very good and the extrapolated continuum limit values agree very well. Cutoff artifacts in B RGI K,lat (μ * ℓ ,μ * h ) are about 10% with the method M1 and a bit larger if the method M2 is adopted.
Two remarks are in order here. First of all, our choice of the lattice normalization factor (M 12 F 12 M 34 F 34 ) −1 in the r.h.s. of eq. (2.20) turns out to be very 5e-03 4e-03 3e-03 2e-03 1e-03 0e+00  K,lat as function of (af 0 ) 2 , at the same reference quark masses,μ * ℓ and µ * h , as in Fig. 2, is shown together with the best fit linear in a 2 and the corresponding continuum limit result. The two sets of data, labelled as M1 and M2, come from different procedures for evaluating the renormalization constants Z AV +V A and Z A , expected to agree in the continuum limit.
beneficial as it partially compensates the discretization errors coming from the matrix element P 34 |Q V V +AA (0)|P 21 . Secondly, the cutoff effects in the RI-MOM scheme computation of Z A and Z V A+AV are significantly reduced when 1-loop perturbative lattice artifacts (i.e. O(a 2 g 2 0 ) effects) are subtracted. The latter have been computed in refs. [32,33].
In conclusion, in spite of the occurrence of a substantial O(a 2 ) artifact in the mass of the lattice state |P 21 and the ensuing fact that the P 34 |Q V V +AA (0)|P 21 matrix element is evaluated with an O(a 2 ) four-momentum transfer, the scaling of B RGI K,lat (μ * ℓ ,μ * h ) with a 2 is found to be well under control, thereby allowing for a reliable continuum limit extrapolation. The same is thus expected, and found to be true (see sect. 3.2.2), also upon approaching the physical quark mass point. pion mass aM tm ℓℓ 7 , the low-energy constants af 0 and aB 0 and the chiral limit value of a/r 0 , are derived from ETMC data [29,34,35]. In this paper the symbolQ means that the quantity Q is renormalized in the MS scheme at the 2 GeV scale. We note that, since in our analysis the continuum limit valueB 0 = Z P B 0 is employed, in order to evaluate the RGI quantity B 0 µ ℓ =B 0μℓ one must knowμ ℓ = Z −1 P µ ℓ at the β-values of interest. The necessary renormalization constants of quark bilinears, such as Z P , Z A and Z V , have been computed non-perturbatively in the RI-MOM scheme [25] as reported in ref. [32] and converted to the MS wherever necessary. The calculation of the four-fermion operator renormalization constant is new and discussed below in subsect. 3.2.1.

Extracting bare B K estimates from lattice data
In the present computation we deal with hℓ pseudoscalar mesons ("kaons") made of two valence quarks, q h and q ℓ , with bare masses (µ h , µ ℓ ). The values of the mass parameters (aµ h , aµ ℓ ) are chosen so as to have ℓℓ pseudoscalar mesons (charged "pions" with r u = −r d ) with mass in the range 280-550 MeV and hℓ pseudoscalar mesons with mass in the range 450-700 MeV, depending also (see sect. 2) on whether the lattice mass M 34 ↔ MK , or rather M 12 ↔ M K , is considered. The value of the light quark mass parameter, aµ ℓ , is common for sea and valence light (the would-be u/d) quark flavours, while the heavier quark (the would-be s) is quenched. With this choice of mass parameters we will eventually be able, as discussed in sect. 3.2, to make contact with the physical K 0 andK 0 mesons by interpolating our results in µ h → µ s and extrapolating them in the µ ℓ → µ u/d and continuum limit.
Simulation details are gathered in Tables 1, 2, and 3, where we report, at each gauge coupling and for all the available (µ h , µ ℓ )-combinations, the values (in lattice units) of the measured pseudoscalar meson masses and decay constants as well as the bare B K -parameter. The latter (up to a trivial factor 8/3, see eq. (2.20)) is denoted by R and it is obtained by averaging over a plateau in τ ≡ x 0 − y 0 of the lattice estimator R(x 0 ) given in eq. (2.10). In the table captions we specify the plateaux intervals [(x 0 − y 0 ) min /a, (x 0 − y 0 ) max /a]. Indeed, the quantity R(x 0 ) must be evaluated at large time separations. In our case this means requiring that x 0 (the time coordinate where the four-fermion operator is located) is sufficiently far away from y 0 + T /2 and y 0 (with y 0 , y 0 + T /2 the time coordinates of the two "K-meson walls", see sect. 2.2). In this way the K 0 -and theK 0 -state dominate the three-point correlator C 3 (x 0 ). Actually, due to the finite T extension of the lattice, it is also necessary that the contribution to the quantum-mechanical representation of C 3 (x 0 ) from states wrapping around-the-world, such as that from |ππ or the a 2 -suppressed one from |π 0 , be negligible compared to the contribution from the vacuum state.
For this purpose it is important that T is large enough with respect to the inverse lattice pion mass(es).
As one can appreciate from Fig. 4, panels (a), (b) and (c), the quality of our signal for R(x 0 ) is quite good. Wider plateaux, as well as stronger suppression of unwanted vacuum-sector states wrapping around-the-world, are obtained when the time extension of the lattice is increased by a factor 4/3 going from a 24 3 × 48 to a 32 3 × 64 lattice at β = 3.9 (i.e. up to T ∼ 5.8 fm). This is illustrated in Fig. 4 panel (d) where the values of R(x 0 ) and its plateau for the larger lattice is displayed. As shown in Fig. 5, no systematic effect in the plateau value of R is visible at β = 3.90 upon comparing results from the lattices 24 3 × 48 and 32 3 × 64. Similarly by comparing data for C 2 (x 0 ) (C ′ 2 (x 0 )) from the two lattices, we could also establish (see Table 2) the absence of significant finite volume effects in the quantities M 12 and F 12 (M 34 and F 34 ). In particular we checked that the marginally significant finite size effect that one can notice in aM 34 at β = 3.9 and aµ ℓ = 0.0040 has a negligible impact on our determination of B RGI K .

Computing Z AV +V A
Renormalization of B K requires in particular the calculation of the four-fermion operator renormalization constant, Z AV +V A , in the mass independent scheme. This task was carried out using the RI-MOM approach. We followed the general strategy presented in ref. [32] employing the two procedures referred there as M1 (or extrapolation method) and M2 (or p 2 -window method). In this way one obtains for each β two estimates of the renormalization constant at the lattice reference scale, which we will denote as Z RI ′ ,Mj AV +V A (a −2 ), j = 1, 2. Using the known NLO anomalous dimension of the four fermion operator in the RI-MOM scheme, the RGI quantities Z RGI,Mj AV +V A are evaluated. We note that the cutoff effects on the lattice estimators of Z AV +V A are significantly reduced by subtracting their 1-loop perturbative expres- , computed in [33]. Further information on the computation of the four-fermion operator renormalization constant, including a numerical check on the absence of wrong chirality mixing, is given in Appendix A.
In the following we will quote results obtained by employing the estimates of Z A and Z RGI AV +V A from the procedure M2, as well as the very precise determination of Z V obtained in ref. [32]. The latter is determined by comparing (matrix elements of) the local current to the exactly conserved one. We immediately notice that upon using the procedure M1 for the evaluation of Z A and Z RGI AV +V A very similar results for B RGI K are obtained, with differences typically as small as 0.002 (that we include in the systematic error) 8 . The values of the renormalization constants entering the present computation of B RGI K can be found in Table 4.

Continuum and chiral extrapolations
Continuum and chiral (µ ℓ → µ u/d ) extrapolations are carried out simultaneously exploiting the SU(2)-χPT based fit ansatz of ref. [36] which, for our computation 8 The consistency of the continuum limit results for the unphysical quantity B RGI K (μ * ℓ ,μ * h ) evaluated with renormalization constants from the procedure M1 or M2 was noted in sect. 2.4. with µ sea = µ ℓ , reads 9 K,lat (μ ℓ ,μ s ) denotes the RGI lattice estimates of the "bag parameter" at renormalized ℓ-quark massμ ℓ and strange quark massμ s . Note that all dimensionful quantities are expressed in units of f 0 . The dimensionless fit parameters B RGI K,χ (μ s ), b(μ s ) and D B (μ s ) are functions ofμ s . The first of them represents the RGI "bag parameter" atμ s in the SU(2)-χPT limit (μ ℓ → 0).
Several comments are in order here. First, as shortly mentioned above, the renormalized (in M S at 2 GeV) quantitiesB 0 andμ u/d and the chiral limit pion decay constant f 0 as well as estimates of a at β = 3.8, 3.9, 4.05, 3) are obtained by an analysis of the data for the mass and the decay constant of the ℓℓ (charged) pseudoscalar meson following the lines of ref. [35] and using the results on the quark mass renormalization constant Z µ = Z −1 P from [32]. Here we just recall that in this analysis based on SU(2)-χPT at NLO we set the physical scale using the experimental value of f π and take into account finite-size effects and O(a 2 ) artifacts. As we rely on NLO SU(2)-χPT formulae, the ℓℓ pseudoscalar meson data are taken only from aμ ℓ region (see Tables 1, 2 and 3) where the impact of possible NNLO corrections has been checked (via the methods of ref. [35]) to be negligible.
The renormalized strange quark massμ s is extracted by analysing our data on the hℓ pseudoscalar meson mass M 34 at β = 4.05, 3.9 and 3.8 with the help of the SU(2) χPT based ansatz [36]  , although not shown in eq. (3.4), were taken into account in our actual fits. These effects, which are checked at β = 3.9, aµ ℓ = 0.0040 to be ≤ 1% and only marginally significant within errors (see Table 2), have been estimated by using resummed χ-PT formulae [37] and subtracted out from the lattice data. The coefficients C M (μ h ) and D M (μ h ) grow almost linearly withμ h in the aforementioned range, in line with the chiral behaviour expected at finite lattice spacing, namely M 2 34 ∼ (μ ℓ +μ h ) as (μ ℓ +μ h ) → 0 [31]. Since at this stage the value ofμ u/d is already known, once the parameters C M (μ h ), c(μ h ) and D M (μ h ) have been determined by the fit of theμ ℓ -dependence based on eq. More details on this analysis will be given in a forthcoming publication [38]. Also our lattice estimators for B RGI K,lat (μ ℓ ,μ h ) turn out to have a weak dependence onμ h . It is thus straightforward to obtain by linear interpolation inμ h their values at the pointμ h =μ s . By exploiting the available determinations of a(β) in physical units from the analysis of the pion sector data, the interpolation of B RGI K,lat (μ ℓ ,μ h ) to the s-quark mass point can be performed separately at each value of β and aμ ℓ . An example of such interpolations is shown in Fig. 6. In this way we obtain the quantity B RGI K,lat (μ ℓ ,μ s ) for the three lattice spacing and the aμ ℓ -values corresponding to the bare parameters in Tables 1, 2 and 3. This set of numbers can be fit to the formula (3.1) where the quantities B RGI K,χ (μ s ), b(μ s ) and D B (μ s ) are the free parameters to be determined. As it clearly appears from Fig. 7 the quality of the resulting fit is very good.

Final result and error budget
The value of B RGI K at the physical point is finally determined by evaluating eq. (3.1) in the continuum limit at the best-fit values of B RGI K,χ (μ s ) and b(μ s ) and setting µ ℓ = µ u/d . In this way we arrive at the result B RGI K = 0.729 (25) . (3.6) The error quoted in parentheses in eq. (3.6) is of statistical origin. It comes from our fitting procedure and reflects the statistical errors on the raw data for the bare matrix elements and the associated renormalization constants. Its estimate has been performed via a standard bootstrap analysis, thereby taking properly into account all possible cross-correlations. As one checks in Tables 1-3, at fixed quark masses and lattice spacing the typical statistical error on the bare lattice estimator of B K is around 1%. Inspection of Table 4 shows that relative error coming from the B Krenormalization factor Z RGI V A+AV /(Z A Z V ) is about 2.0%. The extrapolation to the β = 3.80 continuum limit and the physical pion mass, plus interpolation to the kaon point, finally leads to the error quoted in eq. (3.6).
In order to try to have an idea of the size of the systematic error other analyses have been performed besides the one discussed above.
i) First of all, the whole analysis was repeated on the same set of raw data but choosing a/r 0 extrapolated to the chiral limit [35] (with r 0 the Sommer scale), instead of af 0 , as the scaling variable used to build dimensionless quantities suited for continuum extrapolation. This change turned out to produce an increase of the value of B RGI K of a few permille. ii) Secondly we replaced 2B 0μℓ by the squared ℓℓ-meson mass (M tm ℓℓ ) 2 in ex-  trapolating to the physical pion mass point and we determined the kaon point by interpolating the quantity (M tm hh ) 2 to the physical value of 2M 2 K − M 2 π . As before B RGI K increases by a few permille. Moreover we have tried a fit function which also includes higher order analytical terms added to the SU(2)−χPT formula. We observe a shift equal to 0.007 in the central value of B RGI K . iii) We also tried to perform the extrapolation to the physical pion mass point using a first or second order polynomial in µ ℓ in alternative to the SU (2)χPT-based ansatz (3.1). We obtained good quality fits where the central value of B RGI K gets shifted by +0.021 or +0.001, respectively. Completely analogous results are obtained by employing a first or second order polynomial in (M tm ℓℓ ) 2 . iv) Finally all these alternative analyses were repeated by using the M1 method for the RI-MOM renormalisation in place of M2. In all cases this produced tiny changes of B RGI K (never larger than ±0.004). For instance, for the main analysis discussed in some detail in the present section, a decrease of B RGI K by 0.002 was observed. Further information on the analyses for the estimate of the systematic error is given in Appendix B.
In this way we arrive at the systematic error budget for B RGI K we anticipated in the Introduction, which is detailed below. We get • 0.004 from the RI-MOM renormalization, evaluated by combining the uncertainty due to the spread coming from methods M1 and M2 and the estimated error due the truncation to NLO of the perturbative series in the conversion from the RI-MOM scheme to RGI.
• 0.009 from the uncertainty in controlling/removing cutoff effects. In this figure we combine the neglected O(a 4 ) lattice artifacts which we estimate to be ∼ 1% 10 and the effect related to the choice of the scaling variable used in the analysis (af 0 or a/r 0 ).
• 0.014 from the systematic uncertainty in the extrapolation to the pion mass point, where we add in quadrature 0.004 from varying the choice of the extrapolation variable (µ ℓ or (M tm ℓℓ ) 2 ), 0.011 from the difference between an SU (2)χPT-based fit and a simple polynomial ansatz (we take the average of the spread resulting from the polynomial fits discussed above) and a shift equal to 0.007 due to a possible inclusion of higher order analytical terms in the SU(2)χPT fit formula. .
Summing in quadrature all these small systematic uncertainties leads to a total systematic error estimate of ±0.017. This number was added in quadrature to the statistical error in eq. (3.6), yielding the final result quoted in eq. (1.2).

Conclusions
In this paper we have computed from the N f = 2 tm-LQCD simulation data produced by the ETM Collaboration the value of the strong interaction parameter,B K , which controls the K 0 −K 0 oscillations. To this end we have exploited the partially quenched set up proposed in ref. [22] which, at the price of introducing O(a 2 ) unitarity violations, ensures O(a) improvement and absence of wrong chirality mixing effects [13]. Renormalization is carried out non-perturbatively in the RI-MOM scheme [25,32].
Using data at three lattice spacings and a number of pseudoscalar masses in the interval 280 MeV < m PS < 550 MeV, we get in the continuum limit of N f = 2 QCD and at the physical value of the pion and kaon mass the value of the RGI "bag parameter" reported in eq. (1.2), namely B RGI K = 0.729 ± 0.030. Our determination is quite accurate, with a ∼ 4% total error (with statistical and total systematic errors added in quadrature), and agrees rather well with other existing values. In particular, from the comparison of our result for B K with the results from N f = 2 + 1 dynamical simulations (see Fig. 1), it may be inferred that the quenching of the strange quark leads to an error which is rather small i.e. smaller, at present, than other systematic uncertainties affecting current B K estimates. If taken at face value, our result confirms the tension between the lattice determination (which do not seem to depend significantly on the number of dynamical quark flavours) and the preferred value of recent phenomenological, Standard Model based analyses (see refs. [3,[39][40][41][42]). However some theoretical concerns remain about the way of estimating the long distance contributions in the relation between the strong interaction parameterB K and the experimentally well measured quantity ǫ K that parameterizes the (indirect) CP violation in the K 0 −K 0 system.
In order to get a substantially more accurate determination of B RGI K we would need to increase the statistics and work at even finer lattice spacings with the u/dquark mass closer to the physical point. A further possibly important improvement is the use of data from more realistic simulations where dynamical strange and charm sea quarks are included. Actually the ETM Collaboration is already moving ahead in these directions [43] and work aiming at more precise lattice determinations of the kaon "bag parameter" is in progress.

Appendix A -The RI-MOM computation of Z V A+AV
Here we discuss the non-perturbative computation of the renormalization constant (RC) of the operator Q V V +AA (see eq. (2.3)), which, for the reasons explained in sect. 2.1, we call Z V A+AV . It is first evaluated in the RI'-MOM scheme 11 , at some appropriate scale µ 0 lying in the domain of applicability of (RG-improved) perturbation theory, following the strategy detailed in ref. [32]. Then Z RI ′ V A+AV is 11 The prime in the label RI ′ is to remind the specific defining condition adopted for the quark field RC Zq, which, as convenient in lattice studies, is taken here to coincide with the quark propagator form factor Σ1. For the latter we adopt the definition of eq. (33) of ref. [32].
converted into the (scheme-and scale-independent) quantity Z RGI V A+AV according to the formula which is here evaluated with N f = 2 and to NLO accuracy [44,45]. RGI quantities are normalized following the conventions of ref. [1]. In this appendix we quote results for Z RGI V A+AV but, since the conversion from RI'-MOM to RGI is a straigtforward step, most of our discussion will actually concern the evaluation of Z RI ′ V A+AV (µ 2 0 ) at the three lattice spacings of interest for the present work. The error on B RGI K associated to the NLO perturbative conversion from RI'-MOM to RGI is estimated to be about 0.002 and has been taken into account in the error budget discussion of sect. 3.2.3. This estimate is obtained assuming that the neglected NNLO and higher order relative corrections to the conversion factor are as large as the square of the NLO relative correction, which turns out to be ≃ 5%.
For the purpose of determining Z RI ′ V A+AV , besides the q f -quark propagator 12 we compute the relevant Green function is the momentum space correlator The where lower-case Greek (Latin) symbols denote uncontracted spin (colour) indices. In terms of the corresponding amputated Green function 13 the RI'-MOM RC is obtained by imposing in the chiral limit (µ val = µ sea = 0) the renormalization condition where Z q denotes the (flavour independent) quark field RC and the V V + AA-vertex D lat 11 is given by where P V V +AA αβα ′ β ′ is a suitably normalized spin-projector reflecting the V V +AA Dirac structure of the vertex. Details on this projector and the notation can be found in [15].
Sketch of procedure for extracting Z RGI

V A+AV
Here we summarise the rather standard procedure we follow to extract Z RGI V A+AV , deferring the discussion of further details to subsequent subsections. First of all, for each β and choice of the scalep 2 , the condition (A.6) is enforced at all the values of µ sea and µ val given in Table 2 of ref. [32]. By doing so one obtains lattice approximants at non-zero quark mass(es) of the desired RC, Z RI ′ V A+AV (p 2 ; a 2p2 ; µ val , µ sea ) which are then extrapolated to µ val = µ sea = 0 14 .
Improved and controlled estimates of Z RI ′ V A+AV (p 2 ; a 2p2 ; 0, 0) are obtained by removing the perturbatively leading cutoff effects. Using NLO continuum QCD evolution [44,45], the first argument is brought to a reference scale value, µ 2 0 . The reliability of this step rests on the usual assumption that the scalesp 2 and µ 2 0 are large enough to make NLO-PT accurate (within statistical errors). The remaining a 2p2 dependence in Z RI ′ V A+AV (µ 2 0 ; a 2p2 ; 0, 0) is to be regarded as a mere lattice artifact which will be taken care of by employing either the M1 or M2 method (introduced in ref. [32]), as briefly described below.
In order to reduce the statistical error, the lattice approximants of the Q V V +AA RC have been averaged over equivalent pattern of Wilson parameters (r 1 , r 2 , r 3 , r 4 )  Table 4: Results for Z RGI V A+AV (M1,M2) from the analysis discussed in this appendix, together with the (best) final estimates of Z A (M1,M2) and Z V from ref. [32], for the three β-values of interest here.
of valence quarks (namely (1, 1, 1, −1) and (−1, −1, −1, 1)) as well as over different lattice momenta corresponding to the same value ofp 2 . We have checked that performing these averages before or after taking the chiral limit leads to perfectly consistent (actually almost identical) results. The outcome of the whole analysis is conveniently expressed in terms of Z RGI V A+AV (M1) and Z RGI V A+AV (M2). These are the RGI quantities whose values are quoted in Table 4 together with the values of the other RC's relevant for the B RGI K computation. It is interesting to compare our results with those obtained from a 1-loop perturbative lattice computation. For instance in the case β = 3.9 one gets Z 1−loop,RGI V A+AV = 0.422, if the perturbative gauge coupling is set to g 2 0 / P ( P is the plaquette expectation value), or Z 1−loop,RGI V A+AV = 0.675, if the perturbative gauge coupling is set to the "boosted" value corresponding to the prescription of ref. [46]. These numbers should be compared with the values in Table 4.

Valence chiral limit
For fixed values of β, a 2p2 and aµ sea we fit D 11 (see eq. (A.7)) to the ansatz the form of which we would like now to justify. In general the RI-MOM approach relies on the fact that, at very large values of p 2 (≫ Λ 2 QCD ), the relevant Green functions (and the associated vertices) are polynomial in all the quark mass parameters, because non-perturbative contributions which can violate this property vanish in this limit. However, at finite values of p 2 special care must be taken in studying the quark mass behaviour of Green functions that admit one Goldstone boson (GB) intermediate states in their spectral decomposition. Indeed, such non-perturbative contributions to Green functions, though suppressed by a factor 1/p 2 for each GB pole (see below), are divergent in the chiral limit and must hence be disentangled and removed. As explained in Appendix B of ref. [47], this is precisely the case for the correlator G V V +AA (see eq. (A.4)) and the associated vertex D 11 . By using the LSZ reduction formulae one finds that the GB pole contributions relevant for (the spectral decomposition of) D 11 (p; µ val , µ sea ) are of the form 15 where, for simplicity, all the coefficients that are regular in the GB squared mass and weakly (at most logarithmically) depending on p 2 have been set to unity. Furthermore immaterial (for this discussion) lattice artifact corrections are neglected. As discussed in sect. 2, owing to the choice (2.4) of the valence quark Wilson parameters, the axial currentq 4 γ µ γ 5 q 3 is exactly conserved on the lattice (only broken by soft mass terms), while the conservation ofq 2 γ µ γ 5 q 1 is spoiled by lattice artifacts.
From the form of Q V V +AA it follows in particular that the matrix elements vanish as (M 34 ) 2 ∼ (µ 3 + µ 4 ) in the limit µ 3,4 → 0. This property can be checked by using for these matrix elements the soft pion theorem associated to the conservation of the lattice currentq 4 γ µ γ 5 q 3 in order to reduce the GB-particle P 34 . Hence in the spectral decomposition of D 11 no terms with two GB-poles occur. Moreover, the terms (ii) and (iii) above, which both contain one GB-pole, are strongly suppressed by factors of the kind a 2 /p 2 and (1/p 2 ) 2 , respectively. In line with these arguments a very smooth dependence of D 11 on µ val (or equivalently on M 2 12 ) is observed for all the momenta of interest in our analysis, see Fig. 8 for two typical examples. By fitting the D 11 data to the ansatz (A.9) the valence chiral limit is thus safely determined. A practically identical result is obtained by employing a fit ansatz analogous to (A.9), with µ val substituted by M 2 12 . Combining the valence chiral limit lattice estimator of D 11 with the corresponding approximant of Z q (the valence chiral extrapolation of which poses no problems [32]) leads to reliable estimates of the intermediate quantities Z lat V A+AV (p 2 ; a 2p2 ; 0, aµ sea ).  Figure 9: The quantity Z RI ′ V A+AV (p 2 ; a 2p2 ; 0, aµ sea ), taken at the valence chiral limit, as a function of a 2 µ 2 sea , for a typical lattice momentum choice (see inset) giving a 2p2 ∼ 1.06, for the three β values considered in this paper.

Sea chiral limit
This limit is taken at fixed β and a 2p2 by fitting the Z RI ′ V A+AV (µ 2 0 ; a 2p2 ; 0, aµ sea )-data to a polynomial of first order in a 2 µ 2 sea . Within our statistical errors (typically at the level of 1 to 2%) on the data at fixed aµ sea , the dependence on aµ sea is hardly significant, as one can see e.g. from Fig. 9. We checked that repeating the whole analysis with Z RI ′ V A+AV (µ 2 0 ; a 2p2 ; 0, aµ sea ) fitted instead to a constant in aµ sea leads to compatible results for the desired RC, but with smaller errors and often (but not always, in particular at β = 3.9) acceptable χ 2 -values. Based on these findings we conservatively decided to perform the sea chiral extrapolation via a linear fit in a 2 µ 2 sea . In this way we get the RC-estimators Z RI ′ V A+AV (µ 2 0 ; a 2p2 ; 0, 0). O(g 2 a 2 ) Figure 10: The effect of subtracting from Z RI ′ V A+AV (p 2 ; a 2p2 ; 0, 0) at β = 4.05 (blue dots) the O(a 2 g 2 ) correction (A.10), setting either g 2 = g 2 0 (red squares) or g 2 =g 2 (green diamonds).
The coefficients c (j) q , j = 1, 2, 3 can be found in eq. (34) of ref. [32], while the coefficients c (j) have been recently computed by evaluating the vertex D 11 (p) (eq. (A.7)) in lattice perturbation theory, to 1-loop and including cutoff effects up to the second order in a [33]. Here we just recall that in lattice perturbation theory the vertex D 11 (p) can be expressed in the form .12) and refer to [33] for the values of the (here irrelevant) coefficients b (1,2) D 11 and all details. The form of eq. (A.10) follows directly from eqs. (A.6) and (A.12) above, as well as from eq. (32) of [32], where the perturbative expression of Σ 1 (p) = Z q (p), including the O(a 2 g 2 ) corrections, is given.
In the numerical evaluation of the perturbative correction in eq. (A.10) g 2 was taken as the simple boosted couplingg 2 ≡ g 2 0 / P , where the average plaquette P is computed non-perturbatively. The values of P we employed here are [0.5689, 0.5825, 0.6014] for β = [3.8, 3, 9, 4.05], respectively. The (nice) effect of the (A.10) correction in removing the unwanted a 2p2 dependence is illustrated, in the case β = 4.05, in Fig. 10. In the figure the uncorrected values of Z RI ′ V A+AV (p 2 ; a 2p2 ; 0, 0) are compared with the values of Z RI ′ −impr V A+AV (p 2 ; a 2p2 ) obtained by setting either g 2 = g 2 0 = 6/β or (as we chose in the end) g 2 =g 2 .

Absence of wrong chirality mixings
According to the analysis of ref. [22], at maximal twist wrong chirality mixings can affect the renormalization of Q V V +AA only at order a 2 (or higher). The relation between the renormalized lattice and bare operators has then the form where (using the operator basis of ref. [15]) we have In eq. (A.13) the mixing coefficients ∆ 1j are O(a 2 ) quantities and . . . stand for higher order lattice corrections. In Fig. 11 we plot ∆ 1j , j = 2, 3, 4, 5 as function of a 2p2 for β = 3.8, β = 3.9 and β = 4.05. It is clearly seen that in all cases the mixing coefficients are zero within errors and that this is systematically more so as β increases.

Final estimates from M1 and M2 methods
The first step is to convert Z RI ′ −impr V A+AV (p 2 ; a 2p2 ) to Z RI ′ −impr V A+AV (µ 2 0 ; a 2p2 ) by using the known formula for the NLO running [44,45]. This step is necessary in order to disentangle the O(a 2p2 ) cutoff effects from the genuine continuum p 2 dependence, but the actual value of µ 0 has no impact on the RGI result for the RC. As customary, we take µ 0 = 1/a(β) for each value of β, namely 1/a(3.8) = 2.0 GeV, 1/a(3.9) = 2.3 GeV and 1/a(4.05) = 2.9 GeV.
The method M1 consists in fitting Z RI ′ −impr V A+AV (µ 2 0 ; a 2p2 ) to the ansatz in the large momentum region, 1.0 ≤ a 2p2 ≤ 2.2. As expected, the slope λ V A+AV exhibits a very mild gauge coupling dependence which (unlike the case of quark bilinear RC's) is not significant within our current statistical errors. We thus treat λ V A+AV as a β-independent quantity. According to the ansatz (A.15) a simultaneous linear extrapolation to a 2p2 = 0 at the three β values of interest was performed. The extrapolated values, Z RI ′ V A+AV (µ 2 0 ), are finally used to evaluate via eq. (A.1) to NLO accuracy the quantities Z RGI V A+AV (M1) quoted in Table 4.   Figure 11: The behaviour of the mixing coefficients ∆ 1j , j = 2, 3, 4, 5 as function of a 2p2 for β = 3.8, β = 3.9 and β = 4.05, respectively.
The numbers for Z RI ′ −impr V A+AV (µ 2 0 ; a 2p2 ) at the three β's and the lines of simultaneous best linear fit in a 2p2 are shown in Fig. 12. We recall that in this figure, as in all other figures presented in this Appendix, only data points corresponding to the momenta p satisfying the constraint (A.8) appear. As we said above, data that refer to momenta p violating the condition (A.8) are never used in our RC analysis.
When the method M2 is used, we separately average at each β the values of Z RI ′ −impr V A+AV (µ 2 0 ; a 2p2 ) over a narrow interval of momenta where data are linear in a 2p2 . In physical units we take the same momentum interval for all β's, i.e.p 2 ∈ [8.0, 9.5] GeV 2 .

Appendix B -Chiral and continuum extrapolations
Here we give some details and discuss a few technical aspects of the auxiliary analyses that were performed, as discussed in sect. 3.2.3, to estimate the systematic error on our B RGI K computation. In order to evaluate the systematic errors affecting the extrapolation in the u/d-quark mass, which was performed together with that of the continuum limit, we also tried an analysis where the role of the renormalized quark masses,μ ℓ andμ h , is played by (M tm ℓℓ ) 2 and (M tm hh ) 2 i.e. the masses of the pseudoscalar mesons made out of two mass-degenerate quarks (regularized with opposite Wilson parameters, as implied by the label "tm"). The fit ansatz (3.1) was then replaced by the following one  M tm hh , namely 5.10f 0 , 5.50f 0 , 5.90f 0 . Then three continuum and chiral fits for B RGI K,lat at these three reference values of M tm hh were performed, based on the ansatz (B.1). The fits turned out to be of good quality and yielded continuum limit results for B RGI K (M π , M hh ) at M π = 135 MeV and M hh = (5.10, 5.50, 5.90)f 0 , respectively. For illustration we report in Fig. 13 the results for the case of M tm hh = 5.50f 0 . The extrapolated value, B RGI K (M π , 5.50f 0 ), is represented by an open circle in the figure. As one can see from Fig. 14, the dependence of these extrapolated results on M 2 hh is very mild and hardly significant compared to the 4%-level of our statistical errors. In this situation the interpolation of B RGI K (M π , M hh ) in (M tm hh ) 2 to the physical kaon mass (which we identified as the point where (M tm hh ) 2 = 2M 2 K − M 2 π ) poses no particular problem (see Fig. 14), and yields B RGI K = 0.732 (24). Notice that owing to the weak dependence of B RGI K (M π , M hh ) on (M tm hh ) 2 , the approximation induced by the use of the LO SU(3)-χPT formula to relate (M tm hh ) 2 to the physical value of the 2M 2 K − M 2 π mass difference turns out to have a completely negligible impact on the final result.
With the purpose of estimating the systematic uncertainty involved into the extrapolation of our B K data to the pion (or u/d-quark) physical point, we also studied the effect of performing the chiral fit assuming a first or second order polynomial dependence on µ ℓ (or (M tm ℓℓ ) 2 ). The analysis of theμ h andμ ℓ quark mass while always takingμ s = 92 (5) MeV. In the case where only a linear dependence onμ ℓ is assumed, we set C 2 (μ s ) ≡ 0. Besides B RGI K,χ (μ s ), the other fit parameters in the formula (B.2) are C L (μ s ), C 1 (μ s ) and (only for the case of a second order polynomial) C 2 (μ s ).
In the case of an assumed quadratic dependence of B RGI K,lat data onμ ℓ , in the continuum limit and at the physical u/d-quark mass point one finds B RGI K = 0.730(51). The central value of this determination is perfectly in line with all our previous evaluations, only the associated statistical error is somewhat larger because here we have been dealing with a four (instead of a three) parameter fit. The good quality of the fit is well illustrated by Fig. 15. A good fit yielding B RGI K = 0.750 (26) is also obtained if just a simple linear dependence of the data onμ ℓ is assumed. Not unexpectedly (given the distribution of data points) the central value of B RGI K in this case turns out to be slightly larger than in chiral fits that are based on SU(2) χPT or an ansatz allowing for some curvature in theμ ℓ behaviour. This effect, which is hardly larger than one (statistical) standard deviation, is taken into account in our systematic error budget (see sect. 3 The label M2 has the same meaning as in Fig. 13. For comparison the continuum limit curve, corresponding to the ansatz (3.1) based on SU(2) χPT is also shown.
For completeness we also show in Fig. 16, the chiral and continuum extrapolation of our B RGI K data when the method M1 instead of M2 (see ref. [32] and Appendix A) is employed in the RI-MOM computation of Z A and Z RGI V A+AV . By an analysis otherwise identical to that discussed in sect. 3.2.2 we obtain in the continuum limit and at the physical pion mass point B RGI K = 0.727 (30). As expected from the RC values reported in Table 4, in the present case the cutoff effects at finite lattice spacing are somewhat smaller but the statistical errors a bit larger compared to the case of Fig. 7.
As stated in sect. 3.2.3, the effect of replacing with a/r 0 the scaling variable af 0 to build dimensionless quantities in intermediate steps of our analyses and measure the distance from the continuum limit, is as small as 0.004-0.005. The curves illustrating our interpolation to the s-quark mass point and extrapolation to the u/d-quark point and continuum limit in the case where a/r 0 is used as scaling variable look almost identical (within graphical resolution) to those we presented above and are hence not shown.
Other analyses were also performed corresponding to different variants of our reference analysis, discussed in sect. 3 physical point were found to be very close to the value quoted in eq. (3.6), with a spread that is well covered by our final estimate, ±0.017, of the systematic error.