decay amplitude in the continuum limit.
Abstract
We present new results for the amplitude for a kaon to decay into two pions with isospin : Re; Im. These results were obtained from two ensembles generated at physical quark masses (in the isospin limit) with inverse lattice spacings GeV and GeV. We are therefore able to perform a continuum extrapolation and hence largely to remove the dominant systematic uncertainty from our earlier results Blum:2011ng (); Blum:2012uk (), that due to lattice artifacts. The only previous lattice computation of decays at physical kinematics was performed using an ensemble at a single, rather coarse, value of the lattice spacing [GeV]. We confirm the observation reported in Boyle:2012ys () that there is a significant cancellation between the two dominant contributions to Re which we suggest is an important ingredient in understanding the rule, Re/Re, where the subscript denotes the total isospin of the twopion final state. Our result for implies that the electroweak penguin contribution to is Re(.
pacs:
11.15.Ha, 11.30.Rd, 12.15.Ff, 12.38.GcRBC and UKQCD Collaborations
I Introduction
Nonleptonic decays continue to be an important class of processes in the phenomenology of the standard model of particle physics. Historically it was in these decays that both direct and indirect violation were discovered and the challenges for theoretical physicists include an explanation of the longstanding puzzle of the rule and an ab initio computation of . Developments in the theoretical framework of lattice QCD and in efficient algorithms, together with the availability of the latest computing power, have made meeting these challenges feasible. A significant element of the current joint research program of the RBC and UKQCD collaborations is the evaluation of the amplitudes and , where the subscript represents the isospin of the twopion final state (which by Bose symmetry is restricted to 0 or 2). In this paper we present our latest results for .
In Blum:2011ng (); Blum:2012uk () we reported on the first results from a lattice determination of the amplitude for decays, where is the total isospin of the twopion final state:
(1) 
This was the first quantitative calculation of an amplitude for a realistic hadronic weak decay and hence extended the framework of lattice simulations into the important domain of nonleptonic weak decays. As explained in the Introduction of Blum:2012uk (), in order to obtain the result in Eq. (1) it was necessary to overcome a number of theoretical problems and exploit recent improvements in algorithms and the opportunities provided by increases in computing resources. The systematic errors in (1) are dominated by the fact that the calculation was performed at a single, rather coarse, value of the lattice spacing (fm). We estimated these errors to be .
In this paper we repeat the calculation at two finer values of the lattice spacing and perform the continuum extrapolation.The simulations are carried out at physical pion masses (with unitary sea and valencequark masses) using our two new ensembles with lattice spacings fm and fm. Our new result is presented in Eq. (63) and we reproduce it here for the reader’s convenience:
(2) 
A very interesting feature of our earlier calculation of was the observation that the two dominant contributions to Re show a significant numerical cancellation Boyle:2012ys (). We argued in Boyle:2012ys () that this cancellation is an important element in the explanation of the rule, Re/Re. We confirm this cancellation in the present calculation. Of course, before we can claim that we fully understand the =1/2 rule, we need to compute at physical quark masses and momenta; this calculation is even more challenging than the evaluation of but is under way. For the status of this calculation we refer the reader to ckdaiqian ().
The structure of the remainder of this paper is as follows. In the next section we present the parameters of the two ensembles used in this calculation. The evaluation of the bare matrix elements and the renormalization of the lattice operators are discussed in Secs. III and IV respectively. We consider finitevolume effects in Sec. V and present an overview of the different sources of systematic uncertainty in Sec. VI. We perform the continuum extrapolation in Sec. VII and present our final result in Eq. (63). Section VIII contains our conclusions and a brief discussion of the prospects for the reduction of the errors in as well as for the calculation of . There is one appendix in which we reproduce the calculation from Lin:2002nq () of the LellouchLüscher factor for finitevolume corrections in the context of chiral perturbation theory. This calculation demonstrates how to disentangle the finitevolume corrections which decrease exponentially with increasing lattice volume (a source of systematic error) from those which decrease as a power of the volume (which are corrected by the LellouchLüscher factor). This calculation also clarifies a misunderstanding of these effects in the literature Aubin:2008vh ().
Ii Details of the simulation
The calculations described below have been performed on two new 2+1 flavor ensembles generated with the Iwasaki gauge action and with Möbius domainwall fermions Blum:2014tka () . The parameters of the ensembles are

with ( GeV);

with ( GeV).
These two ensembles use the Möbius variant of domain wall fermions Brower:2012vk () with a Möbius scale factor . For compactness of notation we will refer to these ensembles as and respectively. The lattice spacing and quark masses were set by choosing the masses of the pion, kaon and the baryon to be equal to their physical values. The corresponding seaquark masses are and , with the residual mass for the ensemble and , and for the ensemble. The two ensembles have approximately the same physical volume with spatial extent fm, enabling the continuum extrapolation to be separated from finitevolume effects which we estimate separately. For more details on these ensembles see Blum:2014tka () and we will return briefly to the determination of the lattice spacings in the context of the continuum extrapolation in Sec. VII.
The results presented below were obtained using 76 gauge configurations on the ensemble and 40 on the ensemble. The large statistical uncertainty one expects with a relatively small number of gauge configurations can be significantly reduced if we perform many measurements on each configuration in which the sources and sinks are simply translated in space and time Blum:2014tka (). Performing multiple measurements on the same configuration offers two important opportunities for increased efficiency. First if we can use a lowmode deflation method such as eigCG Stathopoulos:2007zi () we will be able to amortize the setup costs of such an approach over a large number of inversions. Second we can use the all mode averaging technique Blum:2012uh () and perform most of these many inversions at reduced precision and use a relatively few accurate inversions to determine a correction that guarantees systematic double precision but with an additional (usually small) statistical error that reflects the small number of accurate solves. Specifically for the ensemble, the eigCG method was used in single precision with 600 approximate lowlying eigenvectors and a stopping residual of . The approximate (wall source) propagators were computed on all 96 time slices. The accurate solves used to correct the approximation were computed on time slices 0, 76, 72, 68, 64, 60 and 56 with Conjugate Gradient (CG) stopping residual . (This choice of timeslice separations is not related to the calculation presented here but to an accompanying calculation of Blum:2014tka ().) To ensure that no bias results from the choice of inexact solves for which the correction is calculated, this complete pattern of source time slices for the accurate solves was shifted by a different random time displacement on each configuration. A similar procedure was used on the ensemble but with 1500 low modes and a stopping residual of for the approximate solves and accurate solves on time slices 0, 103, 98, 93, 88, 83, 78 and 73. On both ensembles, the accurate CG solves were also computed using eigCG, exploiting the approximate eigenvectors created during the inaccurate applications of eigCG.
(lattice units)  

(lattice units)  
(MeV)  
(MeV) 
Measurements on the and ensembles are separated by 20 and 40 molecular dynamics (MD) units respectively. In order to study the effects of autocorrelations we bin the data. We find that the effects are small, typically leading to a variation of the statistical errors of less than 10%. The results presented below were obtained after binning the 76 configurations of the ensemble into 19 bins of 4 configurations and the 40 configurations of the ensemble into 8 bins of 5 configurations. The 40 configurations from the ensemble are precisely those used in the global analysis reported in Blum:2014tka (). The 76 configurations from the ensemble include 73 of the 80 used in Blum:2014tka (). We have however, repeated the relevant analysis of Blum:2014tka (), including the determination of the lattice spacings, using precisely the 76 configurations for which we have computed . This makes it possible to compute standard jackknife errors for our physical results which necessarily depend upon the value of the lattice spacing.
The pion () and kaon masses () as well as the energies of the twopion state () obtained on the two ensembles are shown in Table 1. The fitting ranges used for pion and kaon masses as well as two pion energies were from 10 to 86 on the ensemble and from 10 to 118 on the ensemble. These choices were motivated by the plateaus in the effective mass plots shown in Figs. 1  2. The effective mass of the kaon, , is defined numerically by the ratio
(3) 
and the twopion effective mass, , is found by inverting
(4) 
The twopoint correlation functions and are defined explicitly in Eq. (22) below and the differences in the numerator and denominator on the lefthand side of Eq. (4) are introduced to eliminate the constant in Eq. (23).
The pion and kaon masses correspond closely to their physical values. We will explain below that the pions are given a momentum in each of the three spatial directions and from the table we see that with this choice and the matrix elements correspond to the onshell (within statistical errors) decay of a kaon in the centerofmass frame. We now discuss the evaluation of the matrix elements.
Iii Evaluation of the bare matrix elements
decay amplitudes are defined by
(5) 
where is the component of the weak Hamiltonian which changes the strangeness by one unit. The weak Hamiltonian can be separated into short and long distance contributions by using the operator product expansion:
(6) 
where is the Fermi constant, and are CabibboKobayashiMaskawa (CKM) matrix elements, the are all the possible dimension6 operators which contribute to the decay and are the corresponding Wilson coefficients which contain information about the short distance physics. The take the form where is the ratio of CKM matrix coefficients .
In this paper we only consider decays where the twopion final state has total isospin 2. The nonperturbative contribution to the decay amplitude is contained in the matrix elements:
(7) 
There are only three operators which contribute to , which we label according to their chiral transformation properties. We have one (27,1) operator and two electroweak penguin operators labeled (8,8) and , where the subscript mx denotes a color mixed operator. Explicitly, the operators are given by
(8)  
(9)  
(10) 
The subscripts L and R denote the left and righthanded spin structures respectively:
(11) 
The Lorentz indices are understood to be contracted between the two parentheses in each of the operators in Eqs. (8)  (10) and are color indices which are summed from 1 to 3.
Below we will confirm the feature found in our earlier work Blum:2011ng (); Blum:2012uk () that the dominant contribution to Re() comes from the (27,1) operator, while the dominant contribution to Im() in the scheme at 3 GeV comes from the operator. We can now write the expressions for the amplitude, which are
(12) 
The relative factor between the two expressions is due to the different ClebschGordan coefficients.
A major challenge in the calculation of (and even more so in the calculation of ) is to ensure that the pions have physical momenta. In the centerofmass frame with periodic boundary conditions, the ground state for the twopion system has each pion at rest. The evaluation of matrix elements at physical kinematics therefore corresponds to the contribution from an excited twopion state resulting in a considerable loss of precision. We can avoid the necessity of multiexponential fits to extract the excited state contribution by utilizing the technique suggested in Kim:2003xt (); Kim:2005gka () and applied successfully in our original calculation of A2 Blum:2011ng (); Blum:2012uk (): we introduce antiperiodic boundary conditions for the (valence) dquark in all three spatial directions, and periodic boundary conditions for the u and squarks Kim:2003xt (). We then exploit the WignerEckart theorem to relate matrix elements to those for the unphysical transition . The relation is
(13) 
The indices and label the twopion state’s total and third component of isospin respectively. With antiperiodic boundary conditions in three spatial directions, the ground state has total momentum , with each pion having momentum . It can be seen from Table 1 that is very close to on both the and ensembles. (For the smaller physical volume in our original calculation Blum:2011ng (); Blum:2012uk (), we imposed antiperiodic boundary conditions for the quark in two spatial directions in order to achieve .) Note that with both periodic and antiperiodic boundary conditions on the quark, the lowest momentum of the meson is zero; this motivates the use of the WignerEckart theorem to reformulate the calculation to that of a matrix element with a final state.
The operators which appear on the righthand side of Eq. (13), and which correspond to the operators in Eqs. (8)  (10), are
(14) 
To simplify the notation we have dropped the labels and on the operators in Eq. (14); this will be implicit in the following. In this paper we compute the matrix elements of the three operators in Eq. (14).
The factor of in Eq. (13) is a combination of coming from the ClebschGordan coefficients and the WignerEckart theorem, and a further corresponding to the simple choice for the normalization of operators in Eq. (14). The amplitude is given in terms of the matrix elements by
(15) 
Since it is the matrix elements which we compute directly in this paper, we choose the compact notation . The label runs over the three operators in Eq. (14).
iii.1 Evaluation of the correlation functions
The bare matrix elements are obtained from the computation of two and threepoint correlation functions. The threepoint functions are
(16) 
where is one of the three operators in Eq. (14) and and are interpolating operators for the kaon and twopion state respectively. For and we take Coulomb gaugefixed wallsource operators defined as follows:
(17)  
(18) 
where in (18) we have used the cosine momentum sources for the quark:
(19) 
represents the quark field and the components of momenta satisfy . Just as for the quark source in Eq. (17), the quark sources in shown in Eq. (18) are given zero momentum by summing them over the full spatial volume, evaluated in the Coulomb gauge. As explained in Ref. Blum:2012uk () the cosine source described above creates quarks with both signs for each component of the three momentum , for , and . This will then produce pairs of pions with total momentum in each direction of in addition to the desired value of . For the threepoint functions described in Eq. (16), the zero total momentum of the decaying kaon and threemomentum conservation imply that the nonzero  momenta cannot occur. For the twopoint function defined in Eq. (22) below we use a  sink which is different from the source and which explicitly projects onto  states with zero total momentum, as described in Ref. Blum:2012uk (). A further subtlety, not described in that reference, relates to the possible angular momentum of the twopion state. For our two identical bosons which carry equal but opposite momenta, there are actually four possible states given our boundary conditions. Specifically, the which carries may have four possible values for the other momentum components: and . These four states form a fourdimensional representation of the cubic symmetry group, which decomposes into two irreducible representations: a singlet () and a triplet (), out of which only contains an swave contribution. Since the lowest energy level of the finitevolume swave spectrum of the representation is nearly degenerate with the lowest energy level of the dwave spectrum of the representation, it is important that we use the cubically symmetrical source specified in Eq. (19) which couples only to the state of interest.
The spinor and color labels are contracted within each set of square parentheses in Eq. (18). A schematic diagram of the correlation function is shown in Fig. 3.
We have evaluated for a range of values of the sourcesink separations . For the () ensemble we performed the calculations for values of between 24 and 39 (26 and 36). These separations were chosen to be large enough for the plateau region to give a reliable fit and small enough for the aroundtheworld effects to be small. The fitting ranges were chosen to be from 10 to for both ensembles. These choices are motivated by the locations of plateau regions in Fig. 4.
For sufficiently large time separations and , the expected time dependence of is
(20) 
where
(21) 
We have introduced the label “bare” as a reminder that are matrix elements of the bare operators in the lattice regularization which we are using. The renormalization of the operators is discussed in the following section. For illustration, in Fig. 4 we plot computed on each of the two ensembles for . The observed plateaus are a manifestation of the fact that the volumes have been tuned so that [cf. Eq. (20)].
We obtain the matrix elements by fitting Eq. (20), using the values of , , and obtained from fitting (under the jackknife) the correlation functions,
(22) 
which have the following time dependence:
(23)  
(24) 
The “” limit should be understood as taking a sufficiently large time separation so that excited state contributions are negligible. Introducing the constant in Eq. (23) allows one to account for possible aroundtheworld effects in .
As a check, we can also construct the timeindependent ratio of the correlation functions:
(25) 
This ratio is plotted for in Fig. 5. As anticipated, all three operators exhibit a constant behavior in the region where the contribution from excited states is negligible. Equation (25) is expected to hold in the region , where is the total time extent of the lattice. In this region “aroundtheworld” effects arising from different time orderings of the operators can be neglected.
The values of the bare matrix elements are shown in Table 2. The entries have been obtained by performing weighted averages (under the jackknife) over the values obtained for each choice of .
ensemble  

ensemble 
Iv Renormalization of the Operators
Having determined the matrix elements of the bare operators in the lattice regularization we now have to combine them with the remaining factors in Eq. (6) to obtain . The Wilson coefficients [] and composite operators [] appearing in Eq. (6) are separately renormalization scheme and scale () dependent. To obtain the physical amplitudes they must be combined in the same scheme and at the same scale. The are calculated in perturbation theory for which it is convenient to use the NDR scheme (called in the following). NDR stands for “naive dimensional regularization” prescription for the matrix, which preserves the anticommutation relations with other gamma matrices Buchalla:1995vs (). The matrix elements calculated in Sec. III, on the other hand, were obtained using bare operators with the lattice spacing as the ultraviolet regulator with the lattice discretization of QCD. The operators can be renormalized nonperturbatively, but only into schemes for which the renormalization condition can be imposed on lattice Green’s functions. The scheme, which is based on dimensional regularization cannot be simulated in a lattice computation. Our procedure is to start by renormalizing the operators nonperturbatively into schemes which can be simulated, specifically the “regularizationindependent symmetric momentum” (RISMOM) schemes Sturm:2009kb () as described in detail in Blum:2012uk () and briefly summarized below. The matching between the RISMOM and schemes is necessarily performed in perturbation theory and is currently known at oneloop order. (Below we also present the matrix elements in two RISMOM schemes so that if the perturbative coefficients are calculated to higher order in the future, these matrix elements can be used to reduce the systematic uncertainty in due to the truncation of the perturbation series.)
We now briefly summarize the renormalization procedure. We write the fivepoint amputated Green’s functions of the three operators in Eq. (14) as a threecomponent vector , and impose a renormalization condition of the form
(26) 
where is a vector of projectors and the corresponding treelevel matrix. Denoting the treelevel contribution by the superscript and including explicitly the spinor and color labels, the matrix is given by
(27) 
Here greek letters label spinor components, the uppercase roman letters represent color indices and denote the operators and projectors. For illustration, the treelevel value of the Green’s function of is
(28) 
For the renormalization we only consider the parityeven component of the fourquark operators.
The choice of projectors is not unique and we implement two different sets known as the and projectors, given explicitly by
(29) 
and
(30) 
The corresponding matrices read
(31) 
and
(32) 
where is the number of colors.
The final result for the amplitude is, of course, independent of the choice of intermediate scheme defined by , but comparing the results obtained with different projection operators gives us an estimate of the systematic uncertainty due to the truncation of perturbation theory in relating the RISMOM schemes to the schemes.
The renormalized operators are related to the bare ones by a matrix relation of the form
(33) 
In order to extract the renormalization constants we follow the standard procedure Martinelli:1994ty (); Donini:1999sf () and compute numerically the amputated Green’s functions of the bare operators in Eq. (14) with particular choices of external momenta (as discussed below) on Landau gaugefixed configurations. We next solve Eq. (26) which we rewrite in the form
(34) 
where is the quark field renormalization constant and is the renormalization scale, which we ultimately choose to be 3 GeV.
The choice of is also not unique, and we use the following two cases:
(35) 
where is the threepoint amputated Green’s function of the local vector current and is the renormalization constant of the local vector current. In practice, we multiply each side of Eq. (34) by the square of the corresponding side of Eq. (35). This eliminates and after this multiplication the lefthand side of Eq. (34) contains the ratio of renormalization factors . is then calculated by imposing the Ward identity , where is the local vector current and is the state of a pseudoscalar meson at rest with mass ; this is explained in detail in Blum:2014tka ().
The choice of projection operator for the fourquark operator and defines a renormalization scheme, which we will label with for the choice of and . In particular, we consider the (,) and (,) schemes, having found in earlier studies that the perturbative conversion to the scheme is more precise in these schemes. This is based on the observation that the nonperturbative running is generally closer to the perturbative one for these schemes for the fourquark operators in Eq. (14) Blum:2012uk (); Arthur:2011cn (). As explained below, we follow our previous practice and choose the (,) scheme for our central value and the (,) scheme to estimate the error due to the perturbative conversion to the scheme.
Chiral symmetry suppresses mixing of operators in different irreducible representations of the chiral symmetry group, so that if the symmetry is exact, is a block diagonal matrix with a block corresponding to the renormalization of the operator and block corresponding to the mixing of and operators. In a massless renormalization scheme with a chiral discretization such as the domainwall action, we expect a mixing pattern very similar to this, but with a small mixing between the blocks.
The mixing of the operator with either of or due to explicit chiral symmetry breaking induced by finite is proportional to (which is in this work). Such mixing can result from two mechanisms Aoki:2005ga (); Christ:2005xh (). First, both quarks in a lefthanded  pair in can propagate in the fifth dimension from the lefthand to the righthand wall, exploiting numerous but exponentially damped modes which even in perturbation theory link the left and righthand walls. This will change the operator into one transforming as the representation, but requires the propagation of two quarks from the lefthand to the righthand wall. This incurs a penalty of since one power of the residual mass results from the fifthdimensional mixing of the left and righthanded components of a single quark.
The second mechanism is nonperturbative and more subtle. For this case the propagation results from the leftright tunneling that can be caused by an eigenvector of the fivedimensional transfer matrix with a nearunit eigenvalue. Such eigenvectors permit leftright mixing but are rare and therefore give a small contribution to . Under some circumstances such modes can simultaneously allow a number of quark flavors to flip chirality. However, to change a (27,1) representation into an (8,8) one, both a quark and an antiquark must flip chirality which requires two distinct transfer matrix eigenvectors and is therefore also doubly suppressed by a factor . Such doubled suppression will not occur for the mixing between the operator and, for example, an operator in the representation. Here a single transfer matrix eigenvector with nearunit eigenvalue can result in a mixing between and by allowing both a  and a quark (localized near this eigenvector) to flip chirality. This kind of mixing has been studied for example in Boyle:2011kn () and it was found to be largely suppressed by our choice of kinematics, as explained below.
In order to suppress physical infrared chiralsymmetry breaking effects we choose to impose the renormalization conditions with the kinematics indicated in Fig. 6 with . We compute the Green’s functions for several momenta and interpolate to GeV using a quadratic Ansatz. Using partially twisted boundary conditions, we have a good resolution around the targeted momentum. The momenta in such RISMOM schemes are chosen so that there are no “exceptional” channels, i.e. no channels in which the square of the momenta is small Sturm:2009kb (). (This is in contrast with the original RIMOM scheme Martinelli:1994ty (); Donini:1999sf () in which .) We have already checked that with domainwall fermions and this choice of kinematics the chirally forbidden matrix elements are numerically negligible Blum:2012uk (). In the present computation, we use the and ensembles which have physical light and strange seaquark masses. However, the lightquark mass is used in all of the valencequark propagators in the fivepoint Green’s functions, including those for both light and strange quarks. We do not extrapolate either the sea or valencequark masses to zero and, strictly speaking, do not work in the chiral limit. In practice the lightquark masses are sufficiently small that their effects are negligible as is the nonzero mass of the strange sea quark. Comparing our results with those of our previous work (with Shamir domainwall fermions) where a chiral extrapolation was performed we find agreement at the permille level or better.
We find that all the chirally forbidden renormalization factors are smaller than , so we set the corresponding matrix elements of to zero and finally obtain the renormalization matrices:
(36)  
(37) 
for the ensembles and
(38)  
(39) 
for the ensembles. With momentum sources Gockeler:1998ye (), only a few configurations are needed to obtain an excellent statistical precision. The number of Landau gaugefixed configurations used to obtain these results varies between 5 and 15. The statistical errors were estimated with 200 bootstrap samples. The matrices in Eqs. (36) – (39) are the ones used in Eq. (33) to obtain the operators renormalized in the RISMOM schemes at the scale from the corresponding lattice bare operators.
The procedure described above enables us to calculate the matrix elements of the operators in Eq. (14) in the (continuum) RISMOM schemes with a very small systematic uncertainty due to the renormalization. The Wilson coefficients however, are computed in the scheme and so we have to match the RISMOM schemes to the one. We repeat that this matching is perturbative and at present is only known to oneloop order Lehner:2011fz (); this limitation amplifies the uncertainty due to the renormalization. This uncertainty could be reduced by extending the perturbative calculations to higher orders. Future lattice calculations could also help here by using step scaling to run the renormalization constants obtained in the RISMOM schemes nonperturbatively to larger momentum scales. The perturbative matching to the scheme can then be performed at these larger scales where the coupling constant is smaller, leading to smaller uncertainties. We now estimate the current uncertainty due to the matching.
ensembles  ensembles  

Re()  GeV  GeV 
Im()  GeV  GeV 
Re()  GeV  GeV 
Im()  GeV 