decay amplitude in the continuum limit.
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 two-pion 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.Gc
RBC and UKQCD Collaborations
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 long-standing 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 two-pion final state (which by Bose symmetry is restricted to 0 or 2). In this paper we present our latest results for .
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 valence-quark 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:
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 finite-volume 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 Lellouch-Lüscher factor for finite-volume corrections in the context of chiral perturbation theory. This calculation demonstrates how to disentangle the finite-volume 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 Lellouch-Lü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 domain-wall 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 sea-quark 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 finite-volume 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 low-mode 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 low-lying 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 time-slice 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.
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 two-pion 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
and the two-pion effective mass, , is found by inverting
The two-point correlation functions and are defined explicitly in Eq. (22) below and the differences in the numerator and denominator on the left-hand 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 on-shell (within statistical errors) decay of a kaon in the center-of-mass frame. We now discuss the evaluation of the matrix elements.
Iii Evaluation of the bare matrix elements
decay amplitudes are defined by
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:
where is the Fermi constant, and are Cabibbo-Kobayashi-Maskawa (CKM) matrix elements, the are all the possible dimension-6 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 two-pion final state has total isospin 2. The nonperturbative contribution to the decay amplitude is contained in the matrix elements:
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
The subscripts L and R denote the left- and right-handed spin structures respectively:
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
The relative factor between the two expressions is due to the different Clebsch-Gordan 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 center-of-mass frame with periodic boundary conditions, the ground state for the two-pion system has each pion at rest. The evaluation of matrix elements at physical kinematics therefore corresponds to the contribution from an excited two-pion 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) d-quark in all three spatial directions, and periodic boundary conditions for the u- and s-quarks Kim:2003xt (). We then exploit the Wigner-Eckart theorem to relate matrix elements to those for the unphysical transition . The relation is
The indices and label the two-pion 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 Wigner-Eckart theorem to reformulate the calculation to that of a matrix element with a final state.
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 Clebsch-Gordan coefficients and the Wigner-Eckart 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
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 three-point correlation functions. The three-point functions are
where is one of the three operators in Eq. (14) and and are interpolating operators for the kaon and two-pion state respectively. For and we take Coulomb gauge-fixed wall-source operators defined as follows:
where in (18) we have used the cosine momentum sources for the -quark:
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 three-point functions described in Eq. (16), the zero total momentum of the decaying kaon and three-momentum conservation imply that the nonzero - momenta cannot occur. For the two-point 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 two-pion 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 four-dimensional representation of the cubic symmetry group, which decomposes into two irreducible representations: a singlet () and a triplet (), out of which only contains an s-wave contribution. Since the lowest energy level of the finite-volume s-wave spectrum of the representation is nearly degenerate with the lowest energy level of the d-wave 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.
We have evaluated for a range of values of the source-sink 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 around-the-world 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
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,
which have the following time dependence:
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 around-the-world effects in .
As a check, we can also construct the time-independent ratio of the correlation functions:
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 “around-the-world” 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 .
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 non-perturbatively into schemes which can be simulated, specifically the “regularization-independent symmetric momentum” (RI-SMOM) schemes Sturm:2009kb () as described in detail in Blum:2012uk () and briefly summarized below. The matching between the RI-SMOM and schemes is necessarily performed in perturbation theory and is currently known at one-loop order. (Below we also present the matrix elements in two RI-SMOM 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 five-point amputated Green’s functions of the three operators in Eq. (14) as a three-component vector , and impose a renormalization condition of the form
where is a vector of projectors and the corresponding tree-level matrix. Denoting the tree-level contribution by the superscript and including explicitly the spinor and color labels, the matrix is given by
Here greek letters label spinor components, the uppercase roman letters represent color indices and denote the operators and projectors. For illustration, the tree-level value of the Green’s function of is
For the renormalization we only consider the parity-even component of the four-quark operators.
The choice of projectors is not unique and we implement two different sets known as the and -projectors, given explicitly by
The corresponding matrices read
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 RI-SMOM schemes to the schemes.
The renormalized operators are related to the bare ones by a matrix relation of the form
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 gauge-fixed configurations. We next solve Eq. (26) which we rewrite in the form
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:
where is the three-point 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 left-hand 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 four-quark 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 four-quark 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 domain-wall 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 left-handed - pair in can propagate in the fifth dimension from the left-hand to the right-hand wall, exploiting numerous but exponentially damped modes which even in perturbation theory link the left- and right-hand walls. This will change the operator into one transforming as the representation, but requires the propagation of two quarks from the left-hand to the right-hand wall. This incurs a penalty of since one power of the residual mass results from the fifth-dimensional mixing of the left- and right-handed components of a single quark.
The second mechanism is nonperturbative and more subtle. For this case the propagation results from the left-right tunneling that can be caused by an eigenvector of the five-dimensional transfer matrix with a near-unit eigenvalue. Such eigenvectors permit left-right 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 near-unit 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 chiral-symmetry 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 RI-SMOM 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 RI-MOM scheme Martinelli:1994ty (); Donini:1999sf () in which .) We have already checked that with domain-wall 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 sea-quark masses. However, the light-quark mass is used in all of the valence-quark propagators in the five-point Green’s functions, including those for both light and strange quarks. We do not extrapolate either the sea- or valence-quark masses to zero and, strictly speaking, do not work in the chiral limit. In practice the light-quark 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 domain-wall fermions) where a chiral extrapolation was performed we find agreement at the per-mille 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:
for the ensembles and
for the ensembles. With momentum sources Gockeler:1998ye (), only a few configurations are needed to obtain an excellent statistical precision. The number of Landau gauge-fixed 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 RI-SMOM 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) RI-SMOM 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 RI-SMOM schemes to the one. We repeat that this matching is perturbative and at present is only known to one-loop 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 RI-SMOM 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.