Precision computation of the kaon bag parameter

Indirect CP violation in K \rightarrow {\pi}{\pi} decays plays a central role in constraining the flavor structure of the Standard Model (SM) and in the search for new physics. For many years the leading uncertainty in the SM prediction of this phenomenon was the one associated with the nonperturbative strong interaction dynamics in this process. Here we present a fully controlled lattice QCD calculation of these effects, which are described by the neutral kaon mixing parameter B_K . We use a two step HEX smeared clover-improved Wilson action, with four lattice spacings from a\approx0.054 fm to a\approx0.093 fm and pion masses at and even below the physical value. Nonperturbative renormalization is performed in the RI-MOM scheme, where we find that operator mixing induced by chiral symmetry breaking is very small. Using fully nonperturbative continuum running, we obtain our main result B_K^{RI}(3.5GeV)=0.531(6)_{stat}(2)_{sys}. A perturbative 2-loop conversion yields B_K^{MSbar-NDR}(2GeV)=0.564(6)_{stat}(3)_{sys}(6)_{PT}, which is in good agreement with current results from fits to experimental data.


Introduction
Neutral kaon mixing is responsible for indirect CPviolation in K→ππ decays. This violation is quantified by the parameter , which is related to quark flavor mixing parameters and the ratio of hadronic matrix elements where O ∆S=2 = [sγ µ (1 − γ 5 )d][sγ µ (1 − γ 5 )d] (cf. [1] for details). The computation of (1) has some advantages over direct computations of the O ∆S=2 matrix element, such as the partial cancellation of statistical and systematic uncertainties.
Note that a precise determination of B K together with experimental measurements of yields important constraints on the unitary triangle parameters (ρ,η).

Lattice Details
We use the N f =2+1 and N f =3, 2 HEX [2] smeared, tree-level clover improved Wilson ensembles generated for determining light quark masses [3,4]. Out of the five available lattice spacings, we use the four finest covering a range of 0.054 fm < ∼ a < ∼ 0.093 fm (the low momentum cutoff at the largest lattice spacing a≈0.116 fm does not allow for a reliable extraction of the mixing coefficients).
For the N f =2+1 ensembles, the pion masses straddle the physical value. Our lattices have sizes as large as L∼6 fm and finite volume corrections to M π [5] are below 0.5% [4]. We also computed the finite volume corrections to B K using the results from [6]. We found that these effects are even below the 0.3% level and thus fully under control. The N f =3 configurations are used to compute the required renormalization constants nonperturbatively using the RI-MOM method [7,8].

Nonperturbative renormalization
Due to explicit chiral symmetry breaking of the Wilson action, the parity even part of operator A, S, P, T denote vector, axial-vector, scalar, pseudoscalar and tensor ∆S = 1 bilinears respectively. We denote their bare matrix elements K 0 |O i |K 0 by Q i . The renormalization pattern is then given by [8] Q ren Hence, the renormalization matrixZ ij is decomposed into Z ij , which, analogously to the continuum renormalization,  [4]), divided by the same running computed at NLO. We observe agreement of the runnings between µ = 1.8 GeV and 3.5 GeV within a statistical error which grows to 2% at the lower end of the range.
only mixes O 2/3 and O 4/5 respectively, and a correction ∆ jk , quantifying the mixing of different O k due to chiral symmetry breaking. Since O 1 does not mix in the continuum, the only relevant terms in the above expression are Z 11 and ∆ 1k , for k = 2, . . . , 5. Because of (1), the multiplicative renormalization factor for B K is given by We use the nonperturbative running method [4,9,10] to circumvent the window-problem of RI-MOM [7,11] and allow matching to other schemes such as MS and RGI with small perturbative uncertainties. This means that, after an essentially flat extrapolation of Z B K to vanishing quark mass for each lattice spacing, we compute the ratio on the three finest lattices at different µ between 1.8 GeV and 3.5 GeV. This ratio is then extrapolated to the continuum limit yielding the nonperturbative running factor from µ to 3.5 GeV, R RI B K (µ, 3.5 GeV). Fig. 1 shows that our continuum extrapolated results for the running agree with NLO perturbation theory [12,13] in the µ-range considered. Thus we set Z RI can be safely extracted for all four lattice spacings. As an additional improvement, we subtract a contact term from the propagators as described in [14][15][16]. The mixing coefficients ∆ 1k are obtained as described in [8], where the additional subtraction [11] ∆ sub 1k (a, m 1 , m 2 ) = is applied. Here, ∆ 1k (a, m i ) is the mixing coefficient obtained at quark mass m i . This procedure removes O(p −2 ) contributions coming from virtual pion exchanges.
The dominant corrections are then O((ap) 2 ) discretization errors and an O(p −4 ) term attributed to double pion exchanges. We use different fit windows as well as fit functions, which include either an (ap) 2 term or an additional p −4 term to estimate systematic effects coming from this ambiguity. However, all these effects turn out to be very small, as both fit functions give compatible results. We also remove the small remaining quark mass dependence in the same fit. Fig. 2 shows the mixing term ∆ sub 14 prior to and after removing the discretization effects at a≈0.077 fm lattice spacing. The same data at a≈0.054 fm are also shown for comparison. We observe that all mixing coefficients are small and tend to zero as a is decreased.

Matrix Elements
To obtain the matrix elements relevant for the computation of B K , we use color-random U(1) wall sources at t = 0 and t = T /2 [17] and vary the time slice τ of the operator insertion between 1 and T − 1. The relevant operator insertions are O 1...5 . The matrix elements Q i are determined by performing constant fits of the time-symmetrized plateaus in τ as shown in Fig. 3. We use three different fitting ranges in order to estimate systematic effects due to excited states. Combining the results from these fits and the ∆ 1k determined before, we can decompose B K into the contributions from the individual Q i . The chiral symmetry breaking contributions of operators O 2,..., 5 are O(α s ) suppressed relative to that of the leading operator O 1 . However, the corresponding matrix elements Q 2,...,5 are chirally enhanced and grow relative to Q 1 as the SU(3) chiral limit is approached. Thus it is important to control the subtraction of these chiral symmetry breaking contaminations. The good chiral properties of our fermion action mitigate this problem since the contribution of Q 1 largely dominates. As shown in Fig. 4, it is 98.1(1.2)% for a≈0.077 fm, M π ∼120 MeV and m s very close to its physical value.

Extraction of Physical B K
We perform a combined chiral and continuum fit to extract the renormalized B K at the physical mass point and in the continuum limit. Since we simulate at or below the physical light quark masses, we can perform a safe interpolation in the quark masses instead of relying on extrapolation formulas. For this combined fit we choose the following functional form  (4) and extrapolated to the chiral limit at a≈0.077 fm lattice spacing (left). The next panel (middle) shows the same data with the O((ap) 2 ) discretization error removed. The same procedure was applied to ∆ sub 14 at a≈0.054 fm (right). The dashed vertical bars indicate the corresponding fitting regions and the horizontal line corresponds to the extracted ∆ 14 along with its statistical 1σ error band. Note the long plateaus in which the data agree with the fit. The extracted mixing coefficient tends to zero as a is reduced.
where f (x, y) with x=M 2 π and y=2M 2 K − M 2 π describes the quark mass dependence. The generic form of f is f (x, y) = 1 + a 10 x + a 20 x 2 + a 01 y We use in total five different fit forms: three Taylor fits  with different powers in x and y, i.e. with a χ =0 and a 10 , a 01 left free as well as either a 20 or a 11 free or set to zero. Additionally, we apply two SU(2) χPT fits (cf. [18,19]) with a χ =1 and a 10 = a 20 = 0 as well as f 0 , µ, a 01 free and a 11 either left free or kept fixed to zero. All fits have good fit quality and show full agreement with the expected SU(2) chiral behavior for the case a χ =1. Since the ratio (1) is tree-level O(a) improved, d(a) is chosen proportional to either α s a or a 2 , as discussed in [4]. We do not include terms whose coefficients are compatible with zero in our final fits. In addition, we use the expressions given in [6] to correct the data for the remaining small finite volume effects. We further apply two different pion mass cuts of 380 MeV and 340 MeV (cf. [4]). Together with two different fit ranges for pion and kaon mass extractions, different fit ranges and functions for obtaining the mixing terms and three different scales for extracting the renormalized matrix elements, we end up with 5760 different analyses. All of our fit results are very precise and compatible with each other. A sample (Taylor) fit with α s a scale dependence in Fig. 5 shows an essentially flat chiral behaviour and extremely small discretization effects. This is also evident from the continuum extrapolation shown in Fig. 6. Using the method from [4,20] for a controlled determination of all systematics as well as the statistical error, our full nonperturbative main result reads where the individual contributions to the systematic error originate from the subtraction of the mixing terms (0.0021), excited state uncertainties (0.0007), extrapolating to the continuum and interpolating to physical quark masses (both 0.0006), as well as the extraction of the renormalization constant (0.0002). Fig. 7 shows the statistical and systematic error distributions of our values for B RI K (3.5 GeV). Both distributions are fairly symmetric and clearly show that our final result is dominated by statisti-  For the reader's convenience, we convert our main result of (7) into the MS-NDR scheme and into the RGI valuê B K . We do so by using the NLO anomalous dimensions of [12,13] and the beta function at the highest available loop order [21]. It is notoriously difficult to reliably assess the truncation error of a perturbative series, particularly in the 68% probability sense of our systematic error treatment. As the NLO contributions to the conversion factors are 2% and NNLO contributions are typically much smaller [22], we add a rather conservative 1% truncation error, which is larger than a variety of perturbative estimates that we have tried. Because this truncation error does not fall into our fully controlled systematic error framework, we list it separately and do not combine it with  (9) with the value forB K obtained by CKMfitter (ICHEP 10 update to [23], vertical line). The dark and light bands correspond to CKMfitter's 1σ and 2σ confidence intervals, respectively. The results from different N f =2 ( [18], 1st) and N f =2+1 ( [19,24,25], 2nd to 4th) lattice computations are also shown.
In Fig. 8 we compare our result to Standard Model expectations and other recent lattice results. Our result is in good agreement with indirect B K determinations from global Standard Model fits of flavor mixing data obtained by CKMfitter (ICHEP 10 updates to [23]). It is consistent with expectations obtained by UTfit [26] by either including all decays (B K,all =0.94(17)) or neglecting the semileptonic channels (B K,no−sl =0.88 (13)) in the fits. Therefore, we find no evidence for new fundamental contributions to indirect CP-violation in K→ππ decays. This is in-line with the findings of [27]. Moreover, we hope that the high precision of our result will encourage our colleagues responsible for the determination of the other contributions to epsilon to work on reducing their uncertainties.