Calculation of $K \to \pi\pi$ decay amplitudes with an improved Wilson fermion action in a nonzero momentum frame in lattice QCD

We present our result for the $K\to\pi\pi$ decay amplitudes for both the $\Delta I=1/2$ and $3/2$ processes with the improved Wilson fermion action. In order to realize the physical kinematics, where the pions in the final state have finite momenta, we consider the decay process $K({\bf p}) \to \pi({\bf p}) + \pi({\bf 0})$ in the nonzero momentum frame with momentum ${\bf p}=(0,0,2\pi/L)$ on the lattice. Our calculations are carried out with $N_f=2+1$ gauge configurations generated with the Iwasaki gauge action and nonperturbatively $O(a)$-improved Wilson fermion action at $a=0.091\,{\rm fm}$ ($1/a=2.176\,{\rm GeV}$), $m_\pi=260\,{\rm MeV}$, and $m_K=570\,{\rm MeV}$ on a $48^3\times 64$ ($La=4.4\,{\rm fm}$) lattice. For these parameters the energy of the $K$ meson is set at that of two-pion in the final state. We obtain ${\rm Re}A_2 = 2.431(19) \times10^{ -8}\,{\rm GeV}$, ${\rm Re}A_0 = 51(28) \times10^{ -8}\,{\rm GeV}$, and $\epsilon'/\epsilon = 1.9(5.7) \times10^{-3}$ for a matching scale $q^* =1/a$ where the errors are statistical. The dependence on the matching scale $q^*$ of these values is weak. The systematic error arising from the renormalization factors is expected to be around $1.3\%$ for ${\rm Re}A_2$ and $11 \%$ for ${\rm Re}A_0$. Prospects toward calculations with the physical quark mass are discussed.


INTRODUCTION
Calculation of the K → ππ decay amplitudes is very important to quantitatively understand the ∆I = 1/2 rule of the K meson decay and to theoretically predict the direct CP violation parameter (ǫ ′ /ǫ) from the Standard Model. A direct lattice QCD calculation of the decay amplitudes for the ∆I = 3/2 process has been attempted for a long time. The RBC-UKQCD Collaboration presented the results at the physical quark mass in Ref. [1], and those in the continuum limit in Ref. [2], which shows good agreement with experiment.
A direct calculation of the decay amplitudes for the ∆I = 1/2 process has been unsuccessful for a long time, due to large statistical fluctuations coming from the disconnected diagrams. A first direct calculation was reported by the RBC-UKQCD Collaboration in Ref. [3] at the pion mass m π = 422 MeV. They also presented a result at a smaller quark mass (m π = 330 MeV) at Lattice 2011 [4]. They used the domain wall fermion action in their calculations which preserves the chiral symmetry on the lattice.
Expanding on the earlier works by Bernard et al. [5] and by Donini et al. [6], we showed in Ref. [7] that mixings with four-fermion operators with wrong chirality are absent even for the Wilson fermion action for the parity odd process in both the ∆I = 3/2 and 1/2 channels due to CP S symmetry. Therefore, after subtraction of an effect from the lower dimensional operator, a direct calculation of the decay amplitudes with the Wilson fermion action is possible without complications from operators with wrong chirality, as for the case with chirally symmetric lattice actions. A potential advantage with the Wilson fermion action over chirally symmetric lattice actions such as the domain wall action is that the computational cost is generally smaller. Hence, with the same amount of computational resources, a statistical improvement may be expected with the lattice calculation of the decay amplitudes, albeit this point has to be verified by actual calculations. We reported a first result with the Wilson fermion action at the pion mass m π = 280 MeV in Ref. [7].
In both the RBC-UKQCD and our calculations of the ∆I = 1/2 process referred above, the kinematics was a K meson at rest decaying to two zero momentum pions at an unphysical quark mass satisfying m K ∼ 2m π . The RBC-UKQCD Collaboration reported a first direct calculation in the physical kinematics, where the pions in the final state have finite momenta, at the physical quark mass in Ref. [8]. They utilized the G-parity boundary condition in the spacial directions, to extract the decay amplitude from the K → ππ time correlation function for the ground state of the two-pion in the physical kinematics.
In the present work, we consider the K meson decay in the physical kinematics by the non-zero momentum frame imposing the periodic boundary condition in the spacial directions. We calculate the decay amplitudes for the decay process K(p) → π(p) + π(0) with momentum p = (0, 0, 2π/L) on the lattice. We set the light quark mass so that the energy of the initial K meson equals to that of the final two-pion state. An advantage of using nonzero momentum in the periodic boundary condition over the G-parity boundary conditions is that gauge configurations, generated previously with the normal boundary condition for the other study and published on some databases, can be used.
Our calculations have been carried out on the HA-PACS and the COMA at University of Tsukuba, the K computer at the RIKEN Advanced Institute for Computational Science, and the Oakforest-PACS at the Joint Center for Advanced High Performance Computing. This paper is organized as follows. In Sec. II we describe the method of calculation used in the present work. The simulation parameters are also given. We present our results in Sec. III. Conclusions of the present work are given in Sec. IV.

II. METHOD
A.

∆S = 1 weak operators
In this section the renormalization of the weak operators for the Wilson fermion action is reviewed. The details have been shown in Ref. [7]. The effective Hamiltonian of the K → ππ decay can be written as [9] , and z i (µ) and y i (µ) (i = 1, 2, · · · 10) are the coefficient functions at renormalization scale µ. Here we consider the case µ ≤ m c , where three light quarks, up, down and strange, are the active quarks in the theory. The ten operators Q i (µ) (i = 1, 2, · · · 10) denote the ∆S = 1 four-fermion operators renormalized at µ. The K → ππ decay amplitude is obtained by A I = K|H|ππ; I for the two-pion state with the isospin I.
The K → ππ matrix elements K|Q i |ππ; I calculated on the lattice are converted to those in the continuum by the renormalization factor for the operators Q i . In the Wilson fermion action, the chiral symmetry is broken explicitly. Hence mixings among different chiral multiplets are in general allowed, and new operators arise through radiative corrections. However, such a problem is absent for the parity odd part of the operators of Q i due to the CP S symmetry.
The mixing to lower dimensional operators is allowed. From the CP S symmetry and the equation of motion of the quark, there is only one operator with the dimension less than 6, which is This operator also appears in the continuum, but does not yield non-vanishing contributions to the physical decay amplitudes, since it is a total derivative operator. However, this is not valid for the Wilson fermion action due to chiral symmetry breaking, and the operator (2) does give non-zero unphysical contributions to the amplitudes. This contribution should be subtracted non-perturbatively, because the mixing coefficient includes a power divergence due to the lattice cutoff growing as 1/a 2 . We can subtract it as by imposing the following condition [10,11], for each operator Q i . The matrix of the renormalization factor of the subtracted operators Q i has the same structure as in the case in which chiral symmetry is preserved.

B. Simulation parameters
Our calculations are carried out with N f = 2 + 1 gauge configurations generated with the Iwasaki gauge action and non-perturbatively O(a)-improved Wilson fermion action on a 48 3 × 64 at β = 1.9. The hopping parameters are κ ud = 0.137712 for the up and the down quark, and κ s = 0.136400 for the strange quark. We use the same algorithm for a generation of the gauge configurations, and the same value of β and the improved factor C SW = 1.7150 as in Ref. [12].
We measure hadron Green's functions and the K meson decay amplitudes at every 10 trajectories. The total length of the run is 9800 trajectories and the total number of gauge configurations employed for the measurement is 980. We estimate statistical errors by the jackknife method with bins of 20 configurations (200 trajectory). The parameters determined from the spectrum analysis are a = 0.091 fm (1/a = 2.176 GeV) for lattice spacing, La for the pion and the K meson mass.
In the present work, we consider the K → ππ decay in the physical kinematics, where the pions in the final state have finite momenta. In order to realize this, we calculate the decay amplitudes for the decay process K(p) → π(p) + π(0) in the non-zero momentum frame with momentum p = (0, 0, 2π/L) on the lattice. We set the light quark mass so that the energy difference between the initial K meson and the final two-pion state vanishes. On our gauge configurations the energy difference is ∆E ≡ E K − E I ππ = −4.7(1.4) MeV for I = 2 and 9.8(7.8) MeV for I = 0 as shown in the following section. In the present work, we assume that these mismatches of the energy give only small effects to the decay amplitudes. We extract the matrix element K|Q i |ππ; I from the time correlation function for the K → ππ process, where Q i (t) is the subtracted weak operator at the time slice t defined by with the subtracted operator Q i (x, t) at the space-time position (x, t) defined in (3). The operator W K 0 (t) in (6) is the wall source for the K 0 meson with the momentum p = (0, 0, −p) (p = 2π/L) at the time slice t, where the wall source for the quark q = u, d, s with the momentum k = (0, 0, k) is given by We adopt K 0 = −dγ 5 s as the neutral K meson operator, so our correlation function has an extra minus from the usual convention. The operator W I ππ (t 1 , t 2 ) in (6) is the wall source for the two-pion state with the isospin I with the momentum p = (0, 0, p) (p = 2π/L), with k 1 = p = 2π/L and k 2 = 0, wherẽ The operator W π i (k, t) is the wall source for π i meson with the momentum k = (0, 0, k) at the time slice t, with k 1 = k and k 2 = 0, wherẽ The periodic boundary condition is imposed in all directions. The summation over δ, where T = 64 denotes the temporal size of the lattice, is taken in (6) to improve the statistics. The time slice of the K meson is set at t K = 29. The wall source of each pion is separated by four lattice unit according to t 1 = 0 and t 2 = 4 in (6) to improve the statistics following the suggestion by Ref. [8]. The gauge configurations are fixed to the Coulomb gauge at the time slice of the wall source t K + δ, t 1 + δ and t 2 + δ for each δ.
The mixing coefficient of the lower dimensional operator, β i in (3), is obtained from the following ratio of the time correlation function K → 0, in the large t K − t region, where W 0 K 0 (t) is the wall source for the K 0 meson with the zero momentum given from (8) by setting p = 0, and Q i (t) = x Q i (x, t) and Q P (t) = x Q P (x, t), with the operator Q i (x, t) and Q P (x, t) defined in (2) at the space-time position (x, t).
The quark contractions of the K → ππ time correlation function G I i (t, t 1 , t 2 ) in (6) are given by (57), (58) and (59) in Ref. [7]. A subtraction of the the contributions of the vacuum diagram, , is not necessary in the present case of the non-zero momentum frame. Those of the K → 0 time correlation function in (18) are given by (60), (61) and (62) in Ref. [7].
The quark loop at the weak operator, the quark propagator starting from the position of the weak operator and ending at the same position, appears in the quark contractions for the K → ππ and the K → 0 processes. We calculate them by the same method as the previous work in Ref. [7], where the stochastic methods with the hopping parameter expansion technique (HPE) and the truncated solver method (TSM) proposed in Ref. [13] are used. As found in the previous work, we find that the first term of TSM in (77) in Ref. [7] is negligible for all time correlation functions also in the present work, when the stopping condition is set at R < 1.2 × 10 −6 . Thus we adopt N R = 0 and N T = 10 with this stopping condition in the calculation of the quark loop.

D.
Time correlation function for ππ → ππ We calculate two types of time correlation functions for the ππ → ππ process to obtain the normalization factors which are needed to extract the matrix elements K|Q i |ππ; I from the time correlation function G I i (t, t 1 , t 2 ) in (6). These are point-wall and wall-wall time correlation functions, which are defined by The operator (ππ) I (t) is the operator for the two-pion with the isospin I with the momentum with k 1 = −p and k 2 = 0, where The operator π i (k, t) is the operator for π i meson with the momentum k = (0, 0, k), defined by The operator W I ππ (t 1 , t 2 ) (I = 0, 2) is defined by (12) and (13). In the present work, we fix t 1 = 0 and t 2 = 4 in (19) and (20) as in the K → ππ process.
The quark contractions for the ππ → ππ time correlation function are given by (89) and (90) in Ref. [7]. Subtractions of the the contributions of the vacuum diagrams, 0|(ππ) , are not necessary in the present case of the non-zero momentum frame.

A. Energy of two-pion state
The time correlation functions for the ππ → ππ process G I P W (t, t 1 , t 2 ) in (19) and G I W W (t, t 1 , t 2 ) in (20) behave in the large time region as with t π = (t 1 + t 2 )/2, where E I ππ is the energy of the two-pion state with the isospin I, S π = m 2 π + p 2 − m π , A I is a constant whose form is irrelevant, and with ∆ = (t 2 − t 1 )/2. The terms with the constant C I and the D I in (27) and (28) arise from the process in which either of the two pions propagates in the opposite time direction to the other pion (i.e., around-the-world effect for the two-pion operator). The effective energy of the point-wall time correlation function G I P W (t) is plotted in Fig. 1, where the effective energy E ef f at t is obtained by f (t + 4, t π , S π ) · G I P W (t + 1, t 1 , t 2 ) − f (t + 1, t π , S π ) · G I P W (t + 4, t 1 , t 2 ) f (t + 3, t π , S π ) · G I P W (t, t 1 , t 2 ) − f (t, t π , S π ) · G I P W (t + 3, t 1 , with determined value of S π = m 2 π + p 2 − m π from the correlation function of the pion. In the present work we set the source operator for the two-pion at (t 1 , t 2 ) = (0, 4) as mentioned in the previous subsection.
We find plateaus in the time region t ≥ 13 for both I = 0 and 2. The value in the non-interacting two-pion case obtained by adding the effective energy of the pion with the momentum p = (0, 0, 2π/L) and that for the zero momentum, which corresponds to E free = m 2 π + p 2 + m π , is plotted for comparison. We find that the two-pion energy for I = 2 is larger than E free signifying repulsive interaction of the two-pion system, whereas that for I = 0 is smaller showing attractive interaction.
In the extraction of the matrix elements K|Q i |ππ; I from the time correlation function G I i (t, t 1 , t 2 ) in (6), the values of E I ππ and A I ππ are needed. Since the statistical error of the point-wall correlation function G I P W (t, t 1 , t 2 ) is smaller than that for the wall-wall function G I W W (t, t 1 , t 2 ), we first extract the energy E I ππ from G I P W (t, t 1 , t 2 ) with the determined value of S π = m 2 π + p 2 − m π from the correlation function of the pion. We then extract the amplitude A I ππ from G I W W (t, t 1 , t 2 ) by fitting to (28) in the lattice unit, where we adopt the fitting range t = [14,26] for I = 2 and t = [14,20] for I = 0. The energy difference between the initial K meson and the final two-pion state ∆E I=0 = 0.0045(36) ( 9.8(7.8) MeV ) .
In the present work, we assume that these violations of energy conservation yield only small effects to the results for the K → ππ decay amplitudes.

B. K → ππ matrix elements
In order to extract the K → ππ matrix element, we consider an effective matrix element M I i (t), which behaves as M I i (t) = M I i ≡ K| Q i (x, 0) |ππ; I in the time region (t 1 , t 2 ) ≪ t ≪ t K . It can be constructed from the time correlation function G I i (t, t 1 , t 2 ) in (6) by where t π = (t 1 + t 2 )/2. Here, the energy of the K meson E K , and the energy of the two-pion state E I ππ are fixed at the values obtained from the time correlation function of the K meson and the ππ → ππ. The factor (−1) comes from our convention of the K 0 operator in (8).
The constant A K = 0|W K |K 2 / K|K is estimated from the wall-wall propagator of the K meson, with the value A K = 2.3544(48) × 10 10 in the lattice unit. The constant A I ππ is defined by (30) and its value is given by (32) or (33).
The constant F I in (36) is the Lellouch-Lüscher factor in the non-zero momentum frame with the momentum |p| = p given in Refs. [14,15], where V is the lattice volume V = L 3 , δ I (p c ) is the two-pion scattering phase shift for the isospin I at the scattering momentum p 2 c = ((E I ππ ) 2 − p 2 )/4 − m 2 π , and the function φ(q c ) at q c = p c · L/(2π) is defined in Refs. [14,15]. In the non-interacting two-pion case, the factor takes the form (F I ) 2 ≡ (F | free ) 2 = (m π + E π ) 2 E K V , where E π = m 2 π + p 2 and E K = m 2 K + p 2 (p = 2π/L). For the interacting case, we need the value of the first derivative of the phase shift δ I (p c ) with the momentum p c to calculate the factor. However, the phase at only one momentum is calculated in the present work and the derivative cannot be evaluated. In the present work, we estimate the factor with an approximation, neglecting the cubic term, leaving a precise estimation of the factor to study in the future.
The calculated values of the Lellouch-Lüscher factors are tabulated in Table I. We find that the main contribution to the factor comes from the term q c · (∂φ(q c )/∂q c ) while the term with the derivative of the phase shift is small (2% for the I = 2 and 12% for the I = 0 case). We consider from this that our approximation for the phase shift in (38) does not make a large systematic error. Our results for the effective matrix elements for the ∆I = 3/2 process are plotted in Fig. 2 and Fig. 3. We do not find plateaus. Rather they show a large and systematic time dependence in the time region (t 1 = 0, t 2 = 4) ≪ t ≪ (t K = 29). This time dependence can be explained by the effect of the two pions propagating in the opposite time directions (i.e., around-the-world effect for the two-pion operator) as follows.
Including the process in which either of the two pions propagates in the opposite direction to the other pion, the time dependence of the correlation function is given by with t π = (t 1 + t 2 )/2, where the label of the isospin I and i for the operator Q i are omitted. States K(k) and π(k) are the K meson and the pion with the momentum k = (0, 0, k) (k = p, 0, −p, (p = 2π/L)), E K = m 2 K + p 2 and E π = m 2 π + p 2 . The second and third terms in (39), F 2 (t) and F 3 (t), are the contributions from the around-the-world effect for the two-pion operator. The constant A ′ corresponds to the K → ππ matrix element, and B ′ and C ′ to the Kπ → π. The around-the-world effects for the K meson, which are suppressed by exp(−E K T ), are neglected in these equations. From (40), (41) and (42), the time dependence of the effective matrix element M(t) in (36) can be written as with S π = E ππ − E π + m π , where the constant A is the K → ππ matrix element, M I i = K| Q i (x, 0) |ππ; I . The other constants B and C correspond to unphysical matrix elements for the Kπ → π process, including exponential factors for t π and T .
We try to extract the K → ππ matrix elements by fitting the effective matrix elements M I i (t) to the fitting function (43) with the determined value of S π from the correlation function of the pion, regarding A, B, and C as unknown parameters. Results of the fitting are plotted in Fig. 2 and Fig. 3 by curved lines, and the extracted values of the matrix elements A = M I i are shown by straight lines with one sigma band. We adopt the fitting range t = [13,18]. The chi-squares of the fitting take very small values, χ 2 /n.o.d = 9.6 × 10 −2 , 3.0 × 10 −4 , 4.6 × 10 −4 for the M I=2 1,7,8 . Therefore the time dependence of the effective matrix element is explained well by the around-the-world effect for the two-pion operator and a value m π T = 7.55 in the present work is not large enough to neglect this effect. We also find that the K → ππ matrix element can be safely extracted by the fitting with (43). The results of the M I i are also tabulated in Table II, where the relations among the matrix elements, M I=2 10 due to the Fierz identity, are used. We should consider this effect also for the ∆I = 1/2 process as for the ∆I = 3/2 process. In Fig. 4 the effective matrix elements M I=0 i (t) for the LL-type four-fermion operators Q 1,2,3,4,9,10 which have the spin structureq 1 γ µ (1 − γ 5 )q 2q3 γ µ (1 − γ 5 )q 4 , are plotted. We do not observe clear time dependence within the statistical error. In Fig. 5 the matrix elements M I=0 i (t) for the LR-type four-fermion operators Q 5,6,7,8 which have the spin structurē q 1 γ µ (1 − γ 5 )q 2q3 γ µ (1 + γ 5 )q 4 , are plotted. We see large time dependence for these cases.
Ideally one should extract the matrix elements by the three parameter fit with (43) regarding A, B and C as unknown constants as was done for the ∆I = 3/2 process, to control the around-the-world effect for the two-pion operator. But the statistics in the present work is not sufficient and the time range for the fitting cannot be taken large to obtain reliable results of the matrix elements. In the present work, we therefore extract the matrix elements by a constant fit for the LL-type four-fermion operators Q 1,2,3,4,9,10 assuming that the around-the-world effects for these operators are small, given that the results in Fig. 4 do not show clear time dependence. For the LR-type four-fermion operators Q 5,6,7,8 , for which we find a large time dependence, we use a three parameter fit with (43). We adopt the fitting range t = [13,18]. The fits are carried out well and the chi-squares of the fitting are small for all cases. Results of the fitting are shown in Fig. 4 and Fig. 5, and tabulated in Table III. C.
K → ππ decay amplitudes The decay amplitudes can be written in terms of the bare matrix element M I i obtained on the lattice by with the effective Hamiltonian H in (1), where with τ = − (V * ts V td ) / (V * us V ud ). The functions C i (µ) are the coefficient functions at the renormalization scale µ and the functions U ij (q * , µ) are the running factor of the operators Q i from the scale q * to µ, which have been given in Ref. [9]. In the present work we set µ = m c = 1.3 GeV and adopt the Standard Model parameters tabulated in Table IV in estimations of these functions. The factors Z ij (q * a) in (45) are renormalization factors which have been calculated by perturbation theory in one-loop order in Ref. [17]. A nonperturbative determination is not yet available. For the renormalization in the continuum theory, we adopt the modified minimal subtraction scheme (MS) with naive dimensional regularization scheme (NDR). We choose two values q * = 1/a and π/a as the matching scale from the lattice to the continuum theory in order to estimate the systematic error coming from higher orders of perturbation theory. The values of the C i have been evaluated in our previous work [7].
Our final results of the decay amplitudes at m π = 260 MeV and m K = 570 MeV are tabulated in Table V. Here, the direct CP violation parameter ǫ ′ /ǫ is obtained by with ω = ReA 2 /ReA 0 , where the experimental value of the indirect CP violation parameter |ǫ| = 2.228 × 10 −3 is used in the estimation. From Table V we learn that the dependence on q * is negligible for most of the decay amplitudes, but it is very large for ImA 2 . A nonperturbative determination of the renormalization factor is necessary to obtain a reliable result for this value. We find a large enhancement of the ∆I = 1/2 process over that for the ∆I = 3/2 at our quark mass m π = 260 MeV. The RBC-UKQCD Collaboration found a numerical mechanism for the enhancement in Ref. [18]. We confirm this numerical mechanism also in our case. Our result for A 0 , particularly for the imaginary part, still has a large statistical error so that we do not obtain a non-zero result for Re(ǫ ′ /ǫ) over the error. Improving statistics by devising a more efficient operator for the two-pion state is an important work reserved for the future. The contributions of the bare matrix element M I i to the decay amplitude A I (A I (i) in (44)) are tabulated in Table II for the ∆I = 3/2 and in Table III for the ∆I = 1/2 process. We find that the main contribution to ReA 2 comes from the operator Q 1 and Q 2 , and that to ImA 2 from Q 8 . The main contribution to ReA 0 comes from the operator Q 2 and that to ImA 0 from Q 6 .

IV. CONCLUSIONS
We have reported on a calculation of the K → ππ decay amplitudes with Wilson fermion action. In order to realize the physical kinematics, where the pions in the final state have finite momenta, we consider the decay process K(p) → π(p) + π(0) in the non-zero momentum frame with momentum p = (0, 0, 2π/L) on the lattice.
We have been able to show a large enhancement of the ∆I = 1/2 process over that for the ∆I = 3/2 at our quark mass (m π = 260 Mev). However, our result for A 0 , particularly for the imaginary part, still has a large statistical error so that we have not obtained a non-zero result for Re(ǫ ′ /ǫ) over the error. Improving statistics by devising a more efficient operator for the I = 0 two-pion state is necessary for a precise determination of the decay amplitudes. Our calculation is carried out away from the physical quark mass. Calculations near or on the physical quark mass are also necessary. We leave these issue to studies in the future. Figures   FIG. 1: Effective energy of the time correlation function G I P W (t, t 1 , t 2 ) for the ππ → ππ with the isospin I = 0 and I = 2. The source operator for the two-pion is located at (t 1 , t 2 ) = (0, 4). A value in the non-interaction case given by adding the effective energy of the pion with the momentum p = (0, 0, 2π/L) and that for the zero momentum, which corresponds to E free = m 2 π + p 2 + m π , is plotted by blue symbol for a comparison.    (44) ) for q * = 1/a and π/a.   III: Decay amplitude for the ∆I = 1/2 process, following the same convention as in Table II. q * = 1/a q * = π/a i M   [16]). τ = − (V * ts V td ) / (V * us V ud ) and Λ (5) MS is the lambda QCD for N f = 5 theory. The standard representation of the CKM matrix of Ref. [16] is adopted, where the CP violation enters entirely through a complex phase of V td , thus τ .