Numerical Study of $ K^0$--$\bar{K}^0$ Mixing and $ B_K $

TL;DR

This study computes the kaon B-parameter ($ B_K $) using staggered fermions in full QCD, comparing different methods and sources, and finds consistent results with no effect from non-degenerate quark pairs.

Contribution

It introduces a comparison of two computational methods and source techniques for calculating $ B_K $ in lattice QCD, demonstrating their agreement and robustness.

Findings

Renormalized $ B_K $ results agree across methods

Improved wall source effectively isolates pseudo-Goldstone bosons

No effect of non-degenerate quark pairs on $ B_K $

Abstract

We have computed $ B_K $ with staggered fermions, using two different methods: both in the one spin trace form and two spin trace form. Renormalized results in both forms are in good agreement. The numerical simulations were performed on a $ 16^3 \times 40 $ lattice in full QCD with $ \beta = 5.7 $. We also tried an improved wall source method in order to select only the pseudo-Goldstone bosons and compare the numerical results obtained with those from the conventional wall source method. We have studied $ B_K $ with a series of non-degenerate quark anti-quark pairs and saw no effect on $ B_K $, although dramatic effects in the chiral limit were seen on the individual terms making up $ B_K $.

Tables

Table 1. Table 1: Here, we present the numerical results for the vacuum saturation amplitude ( ℳ K V subscript superscript ℳ 𝑉 𝐾 {\cal M}^{V}_{K} ) with the quark mass pair (0.01, 0.01), calculated both in the even-odd wall source method and in the cubic wall source method. All the values in the table have been obtained through the covariance fitting to a constant on the bootstrap ensembles.

	Cubic Wall Source		Even-Odd Wall Source
Fitting Range	${ℳ}_K^V$	${\chi }^2/\left(dof\right)$	${ℳ}_K^V$	${\chi }^2/\left(dof\right)$
$14\le t\le 21$	12.17(66)	0.40(63)	49.2(30)	0.40(68)
$13\le t\le 22$	12.10(58)	0.37(52)	49.3(25)	0.33(53)
$12\le t\le 23$	12.15(60)	0.37(48)	48.7(24)	0.61(65)
$11\le t\le 24$	12.10(63)	0.45(48)	49.0(25)	0.68(56)
$10\le t\le 25$	11.76(59)	1.06(51)	47.0(27)	1.24(76)
$9\le t\le 26$	11.70(59)	1.27(66)	46.7(27)	1.34(70)
$8\le t\le 27$	11.61(58)	1.31(62)	45.2(24)	1.40(68)
$7\le t\le 28$	11.20(58)	1.65(70)	44.2(24)	1.52(70)
$6\le t\le 28$	10.74(60)	1.72(65)	42.9(24)	1.47(67)

Table 2. Table 2: Here, we present the numerical results for the ℳ K subscript ℳ 𝐾 {\cal M}_{K} with the quark mass pair (0.01, 0.01), calculated both in the even-odd wall source method and in the cubic wall source method. All the values in the table has obtained through the covariance matrix fitting to a constant on the bootstrap ensembles.

	Two Spin Trace Form				One Spin Trace Form
Fitting	Cubic Wall Source		Even-odd Wall Source		Even-odd Wall Source
Range	${ℳ}_K$	${\chi }^2/(d.o.f.)$	${ℳ}_K$	${\chi }^2/(d.o.f.)$	${ℳ}_K$	${\chi }^2/(d.o.f.)$
$14\le t\le 21$	7.20(48)	1.11(99)	29.8(15)	1.05(102)	33.4(21)	1.45(106)
$13\le t\le 22$	7.23(29)	0.88(87)	29.9(15)	0.87(83)	33.6(17)	1.18(82)
$12\le t\le 23$	7.26(28)	0.75(75)	28.7(14)	1.47(80)	34.0(17)	1.19(68)
$11\le t\le 24$	7.19(28)	0.73(67)	28.6(13)	1.24(70)	34.1(18)	1.03(62)
$10\le t\le 25$	7.16(29)	0.66(60)	28.1(14)	1.29(59)	33.7(16)	0.96(58)
$9\le t\le 26$	7.11(30)	0.62(51)	28.1(13)	1.14(54)	33.7(16)	1.04(49)
$8\le t\le 27$	7.11(30)	0.60(52)	27.7(12)	1.27(61)	33.6(16)	0.95(47)
$7\le t\le 28$	7.04(27)	0.57(54)	27.1(11)	1.47(52)	33.9(16)	0.95(43)
$6\le t\le 29$	6.97(26)	0.61(50)	27.3(12)	1.52(54)	33.8(15)	0.90(45)

Table 3. Table 3: Here, we present the optimal fitting range with respect to quark mass pairs. The optimal fitting range implies that we can get the same average with smaller error bar for χ 2 ≤ 1.0 superscript 𝜒 2 1.0 \chi^{2}\leq 1.0 .

	Fitting Range
Quark Mass	Two Spin Trace Form		One Spin Trace Form
Pair	Cubic Wall Source	Even-odd Wall Source	Even-odd Wall Source
(0.004, 0.01)	$9\le t\le 26$	$13\le t\le 22$	$13\le t\le 22$
(0.004, 0.02)	$11\le t\le 24$	$13\le t\le 22$	$13\le t\le 22$
(0.004, 0.05)	$14\le t\le 21$	$15\le t\le 20$	$15\le t\le 20$
(0.01, 0.01)	$11\le t\le 24$	$13\le t\le 22$	$13\le t\le 22$
(0.01, 0.03)	$11\le t\le 24$	$13\le t\le 22$	$13\le t\le 22$
(0.02, 0.02)	$13\le t\le 22$	$13\le t\le 22$	$13\le t\le 22$
(0.03, 0.03)	$13\le t\le 22$	$13\le t\le 22$	$13\le t\le 22$

Table 4. Table 4: K 𝐾 K meson mass and unrenormalized B K subscript 𝐵 𝐾 B_{K} versus quark mass pair: the B K subscript 𝐵 𝐾 B_{K} data are analyzed by the single-elimination jack-knife method over the optimal fitting range. The K 𝐾 K meson mass, m K subscript 𝑚 𝐾 m_{K} was obtained by analyzing the results for ⟨ 0 ∣ S ¯ ​ ( γ 0 ​ γ 5 ⊗ ξ 5 ) ​ D ∣ K ⟩ quantum-operator-product 0 ¯ 𝑆 tensor-product subscript 𝛾 0 subscript 𝛾 5 subscript 𝜉 5 𝐷 𝐾 \langle 0\mid\bar{S}(\gamma_{0}\gamma_{5}\otimes\xi_{5})D\mid K\rangle .

		unrenormalized $B_K$
Quark Mass	$m_K$	Two Spin Trace Form		One Spin Trace Form
Pair		Cubic Source	Even-odd Source	Even-odd Source
(0.004, 0.01)	0.219(2)	0.557(32)	0.547(54)	0.656(73)
(0.004, 0.02)	0.277(2)	0.641(32)	0.630(39)	0.698(58)
(0.004, 0.05)	0.406(2)	0.731(34)	0.687(44)	0.759(88)
(0.01, 0.01)	0.253(2)	0.600(27)	0.579(33)	0.688(46)
(0.01, 0.03)	0.348(2)	0.710(26)	0.689(30)	0.748(37)
(0.02, 0.02)	0.347(1)	0.709(25)	0.689(28)	0.746(33)
(0.03, 0.03)	0.421(1)	0.768(23)	0.753(26)	0.781(28)

Table 5. Table 5: Unrenormalized B K subscript 𝐵 𝐾 B_{K} for each quark mass pair. The B K subscript 𝐵 𝐾 B_{K} data is analyzed by the covariance fitting over the optimal fitting range on the bootstrap ensembles.

	Unrenormalized $B_K$
Quark Mass	Two Spin Trace Form		One Spin Trace Form
Pair	Cubic Source	Even-odd Source	Even-odd Source
(0.004, 0.01)	0.546(23)	0.548(42)	0.675(63)
(0.004, 0.02)	0.627(28)	0.663(37)	0.700(47)
(0.004, 0.05)	0.732(30)	0.715(40)	0.771(78)
(0.01, 0.01)	0.595(22)	0.607(34)	0.682(45)
(0.01, 0.03)	0.717(23)	0.708(31)	0.771(37)
(0.02, 0.02)	0.727(24)	0.707(29)	0.773(36)
(0.03, 0.03)	0.791(21)	0.758(26)	0.804(25)

Table 6. Table 6: We compare our numerical results for naive B K subscript 𝐵 𝐾 B_{K} (cubic wall source, two spin trace form) with those of other groups (Kilcup and Ukawa et al. ). β = 5.7 𝛽 5.7 \beta=5.7 . m sea ​ a = 0.01 subscript 𝑚 sea 𝑎 0.01 m_{\rm sea}a=0.01 . JN implies that the errors are obtained through single-elimination jack-knife method. BS means that the data is analyzed by covariance fitting on the bootstrap ensembles.

	Naive $B_K$
$m_{\mathrm{valence}}a$	this work (JN)	this work (BS)	Kilcup	Ukawa et al.
0.01	0.600(27)	0.595(22)	0.658(18)	0.69(2)
0.02	0.709(25)	0.727(24)	0.771(11)	0.75(1)
0.03	0.768(23)	0.791(21)	0.818(09)	0.79(1)

Table 7. Table 7: We compare our numerical results of B A ​ 1 subscript 𝐵 𝐴 1 B_{A1} and B A ​ 2 subscript 𝐵 𝐴 2 B_{A2} (cubic wall source, two spin trace form) with those of Kilcup. β = 5.7 𝛽 5.7 \beta=5.7 . m sea ​ a = 0.01 subscript 𝑚 sea 𝑎 0.01 m_{\rm sea}a=0.01 . The errors are estimated through the standard jack-knife procedure.

	$B_{A1}$		$B_{A2}$
$m_{\mathrm{valence}}a$	this work	Kilcup	this work	Kilcup
0.01	1.210(64)	1.225(60)	1.053(46)	1.155(24)
0.02	0.565(23)	0.575(15)	0.860(30)	0.932(13)
0.03	0.407(14)	0.418(7)	0.812(24)	0.866(10)

Table 8. Table 9: The full QCD calculation of naive (unrenormalized) B K subscript 𝐵 𝐾 B_{K} (cubic wall source, two spin trace form) is compared with that of quenched QCD (Sharpe et al. , Kilcup and the Ukawa et al. group). JN means that the errors are estimated through the standard jack-knife procedure. BS means that the data is analyzed by covariance fitting on the bootstrap ensembles.

	Naive $B_K$
$m_{\mathrm{valence}}a$	full QCD ($N_f=2$)		quenched QCD
	this work (JN)	this work (BS)	Sharpe et al.	Kilcup	Ukawa et al.
0.01	0.600(27)	0.595(22)	0.68(2)	0.697(29)	0.69(2)
0.02	0.709(25)	0.727(24)	0.73(1)	0.749(16)	0.74(1)
0.03	0.768(23)	0.791(21)	0.78(1)	0.777(13)	0.78(1)

Table 9. Table 10: Covariance fitting of unrenormalized B K subscript 𝐵 𝐾 B_{K} with 3 parameters: Here 2TR and 1TR represent the two spin trace form and the one spin trace from respectively. Cubic and Even-Odd imply the cubic wall source and the even-odd wall source respectively and B K ​ ( m q ) subscript 𝐵 𝐾 subscript 𝑚 𝑞 B_{K}(m_{q}) means interpolation to the physical quark mass.

Trace Form (Source)	${\alpha }_1$	${\alpha }_2$	${\alpha }_3$	${\chi }^2$/d.o.f.	$B_K\left(m_q\right)$
2TR (Cubic)	0.336(65)	-11.2(33)	-24.1(93)	1.05	6.600(21)
2TR (Even-Odd)	0.346(89)	-10.1(39)	-21.9(112)	1.87	6.375(29)
1TR (Even-Odd)	0.477(119)	-9.35(482)	-21.9(137)	0.47	7.270(30)

Table 10. Table 12: Covariance fitting of unrenormalized B K subscript 𝐵 𝐾 B_{K} with 2 parameters: Here 2TR and 1TR represent the two spin trace form and the one spin trace from respectively. Cubic and Even-Odd imply the cubic wall source and the even-odd wall source respectively and B K ​ ( m K ) subscript 𝐵 𝐾 subscript 𝑚 𝐾 B_{K}(m_{K}) means interpolation to the physical K meson mass.

Trace Form (Source)	${\alpha }_1$	${\alpha }_2$	${\chi }^2$/d.o.f.	$B_K\left(m_K\right)$
2TR (Cubic)	0.293(14)	0.397(66)	0.897	0.638(21)
2TR (Even-Odd)	0.289(21)	0.334(82)	1.53	0.619(28)
1TR (Even-Odd)	0.346(24)	0.191(84)	0.70	0.702(35)

Equations

$$\epsilon =C{\hat{B}}_{K}F(V_{CKM},\frac{m_t^2}{M_W^2},\frac{m_c^2}{M_W^2},\frac{m_u^2}{M_W^2}),$$

$$C=\frac{G_F^2{f}_K^2{m}_K{M}_W^2}{6\sqrt{2}{\pi }^2[m(K_L)-m(K_S)]}$$

$$B_K=\frac{⟨K^0\mid [\bar{s}{\gamma }_{\mu }(1-{\gamma }_5)d][\bar{s}{\gamma }_{\mu }(1-{\gamma }_5)d]\mid K^{\prime 0}⟩}{\frac{8}{3}⟨K^0\mid \bar{s}{\gamma }_{\mu }{\gamma }_{5}d\mid 0⟩⟨0\mid \bar{s}{\gamma }_{\mu }{\gamma }_{5}d\mid K^{\prime 0}⟩}$$

$${\hat{B}}_K=\xi ({\alpha }_s(\mu ))B_K(\mu ).$$

$$\sum _x(D+M{)}_{(z,x)}\cdot G_P(x,y)={\delta }_{z,y}$$

$$\chi (x)=\frac{1}{2}\sum _y(G_P(x,y)+G_A(x,y))\cdot h(y),$$

$$\chi (x+N_t\hat{t})=\frac{1}{2}\sum _y(G_P(x,y)-G_A(x,y))\cdot h(y)$$

$$\sum _x(D+M{)}_{(z,x)}\cdot \chi (x)=h(z)$$

$$\sum _x(D+M{)}_{(z,x)}\cdot \chi (x+N_t\hat{t})=0.$$

$$\bar{\mathbf{a}}\in \left\{(0,0,0),(1,0,0),(0,1,0),\cdots ,(1,1,1)\right\}.$$

$$W_{\bar{\mathbf{a}}}(\vec{n})\equiv \sum _{\vec{m}\in Z^3}{\delta }_{\vec{n},2\vec{m}+\bar{\mathbf{a}}}.$$

$${\mathcal{O}}_{\mathrm{source}}=\sum _{\vec{n},{\vec{n}}^{\prime },\bar{\mathbf{a}}}{\bar{\chi }}^a(\vec{n},t=0)W_{\bar{\mathbf{a}}}(\vec{n})ϵ(\bar{\mathbf{a}}){\delta }_{a,b}{W}_{\bar{\mathbf{a}}}^{†}({\vec{n}}^{\prime }){\chi }^b({\vec{n}}^{\prime },t=0),$$

$${\mathcal{O}}_{\mathrm{sink}}(y)=\bar{\chi }(y_A){\overline{({\gamma }_S\otimes {\xi }_F)}}_{AB}\chi (y_B).$$

$$-⟨{\mathcal{O}}_{\mathrm{sink}}(y){\mathcal{O}}_{\mathrm{source}}⟩$$

$$=\sum _{\vec{n},{\vec{n}}^{\prime },\bar{\mathbf{a}}}{\overline{({\gamma }_S\otimes {\xi }_F)}}_{AB}{G}^{c,a}(y_B,\vec{n})W_{\bar{\mathbf{a}}}(\vec{n})ϵ(\bar{\mathbf{a}}){\delta }_{a,b}{W}_{\bar{\mathbf{a}}}^{†}({\vec{n}}^{\prime })G^{b,c}({\vec{n}}^{\prime },y_A)$$

$$=\sum _{\vec{n},{\vec{n}}^{\prime },\bar{\mathbf{a}}}{\overline{({\gamma }_5{\gamma }_S\otimes {\xi }_5{\xi }_F)}}_{AB}{G}^{c,a}(y_B,\vec{n})W_{\bar{\mathbf{a}}}(\vec{n}){\delta }_{a,b}{W}_{\bar{\mathbf{a}}}^{†}({\vec{n}}^{\prime })[G^{c,b}(y_A,{\vec{n}}^{\prime }){]}^{\ast }$$

$$={\overline{({\gamma }_5{\gamma }_S\otimes {\xi }_5{\xi }_F)}}_{AB}\sum _{c,b,\bar{\mathbf{a}}}{\psi }_{\bar{\mathbf{a}}}^{cb}(y_B)[{\psi }_{\bar{\mathbf{a}}}^{cb}(y_A){]}^{\ast }$$

$${\psi }_{\bar{\mathbf{a}}}^{cb}(y_A)\equiv \sum _{\vec{n}}{G}^{c,b}(y_A,\vec{n})W_{\bar{\mathbf{a}}}(\vec{n}).$$

$$B_K(\mu )\equiv \frac{⟨{\bar{K}}^0\mid \bar{s}{\gamma }_{\mu }(1-{\gamma }_5)d\bar{s}{\gamma }_{\mu }(1-{\gamma }_5)d\mid K^0⟩}{\frac{8}{3}⟨{\bar{K}}^0\mid \bar{s}{\gamma }_{\mu }{\gamma }_{5}d\mid 0⟩⟨0\mid \bar{s}{\gamma }_{\mu }{\gamma }_{5}d\mid K^0⟩.}$$

$$A_{\mu }={\bar{\chi }}_s(\overline{{\gamma }_{\mu 5}\otimes {\xi }_5}){\chi }_d.$$

$${\mathcal{O}}_{2TR}^{Latt}$$

$${\mathcal{O}}_{1TR}^{Latt}$$

$${\mathcal{O}}_{\mathrm{chiral}\mathrm{partner}}^{1TR}$$

$${\mathcal{M}}_K^V\equiv \frac{8}{3}⟨{\bar{K}}^0\mid \bar{s}{\gamma }_{\mu }{\gamma }_{5}d\mid 0⟩⟨0\mid \bar{s}{\gamma }_{\mu }{\gamma }_{5}d\mid K^0⟩.$$

$${\mathcal{M}}_K\equiv ⟨{\bar{K}}^0\mid \bar{s}{\gamma }_{\mu }(1-{\gamma }_5)d\bar{s}{\gamma }_{\mu }(1-{\gamma }_5)d\mid K^0⟩$$

$$g_{MF}^2\equiv \frac{g_0^2(a)}{\mathrm{Re}⟨\frac{1}{3}\mathrm{Tr}{U}_{\square }⟩}$$

$$g_{\overline{MS}}^2({\mu }_{\overline{MS}})=g_{MF}^2[1-{\beta }_0{g}_{MF}^2\left\{2\ln (a{\mu }_{\overline{MS}})+2\ln \left(\frac{{\Lambda }_{\mathrm{Latt}}}{{\Lambda }_{\overline{MS}}}\right)\right\}$$

$$-\frac{1}{3}{g}_{MF}^2+O(g_{MF}^4)]$$

$${\beta }_0=\frac{11}{16{\pi }^2}\left(1-\frac{2N_f}{33}\right).$$

$$g_{MF}^2=1.82$$

$$g_{\overline{MS}}^2\left(\frac{\pi }{a}\right)=1.78.$$

$$g_{\overline{MS}}^2({\mu }_{\overline{MS}})=\frac{g_0^2(a)}{1-t{\beta }_0{g}_0^2(a)}$$

$$\text{mbox}wheret=-2\left[\ln (a{\mu }_{\overline{MS}})+\ln \left(\frac{{\Lambda }_{\mathrm{Latt}}}{{\Lambda }_{\overline{MS}}}\right)\right],$$

$$g_0^2(a)=1.05263$$

$$g_{\overline{MS}}^2\left(\frac{\pi }{a}\right)=1.61.$$

$$B_{V1}$$

Full text

1 Introduction

In the standard model, there are two kinds of CP violation: the indirect CP violation and the direct CP violation. [1, 2, 3, 4]. The indirect CP violation comes from the fact that electro-weak $K^{0}$-$\bar{K}^{0}$ mixing causes the neutral kaon eigenstates not to respect CP symmetry [5]. The direct CP violation comes from the CP violating effective operators (e.g. penguin operators) [3, 4], which can directly contibute to the K decay amplitude. The indirect (direct) CP violation is parametrized by the phenomenological quantity $\varepsilon$, while the direct CP violation is parametrized by $\varepsilon^{\prime}$ [1, 2, 3, 4]. Both $\varepsilon$ and $\varepsilon^{\prime}$ can be measured experimentally.

The low energy effective Hamiltonian of the electro-weak interaction is derived by decoupling heavy particles such as $W$ and $Z$ bosons and the t, b, c quarks in the standard model [3, 6, 7, 8]. The coefficients in this effective Hamiltonian are functions of the Cabibbo-Kobayashi-Maskawa flavor mixing matrix elements $V_{CKM}$. This low energy effective Hamiltonian includes a $\Delta S=2$ term which belongs to the (27,1) representation of the $SU(3)_{L}\otimes SU(3)_{R}$ flavor symmetry group and to the (1,0) representation of the $SU(2)_{L}\otimes SU(2)_{R}$ isospin symmetry subgroup. This $\Delta S=2$ electro-weak effective Hamiltonian causes the $K^{0}$-$\bar{K}^{0}$ mixing in the standard model [3, 6, 7, 8]. Therefore, the $\Delta S=2$ effective Hamiltonian is connected to the $\varepsilon$ parameter. Using this technique, $\varepsilon$ can be expressed in terms of the standard model parameters as follows:

[TABLE]

where

[TABLE]

Here $\hat{B}_{K}$ is defined as a renormalization-group invariant quantity, and the function $F(V_{CKM},m_{t}^{2}/M_{W}^{2},\cdots)$ is given in Ref. [1, 3, 4]. Once $B_{K}$ is determined theoretically, we can narrow the domain of $\mid V_{td}\mid$ and the top quark mass [1], using the experimental determination of $\varepsilon$.

There have been a variety of methods used to calculate $B_{K}$, including chiral perturbation theory, hadronic sum rules, QCD sum rules, $1/N_{C}$ expansion and lattice gauge theory. Lattice gauge theory has the virtue of making the smallest number of assumptions and is exactly equivalent to QCD in the limit of infinite volume and vanishing lattice spacing. For this reason, we have adopted lattice gauge theory to obtain $B_{K}$.

In order to calculate $B_{K}$ in lattice gauge theory, one needs to find a set of operators which can describe on the lattice the same physics as the continuum $\Delta S=2$ four-fermion operator. There have been two methods to implement fermions in lattice gauge theory: the staggered fermion method and Wilson fermion method. For the weak matrix elements involving the pseudo-Goldstone bosons, it is very useful to take advantage of the exact $U_{A}(1)$ symmetry of the staggered fermion action, which is not manifest in Wilson fermion action [9, 10]. We have used the staggered fermion method to calculate $B_{K}$ because of this advantage.

There are two methods to transcribe the continuum $\Delta S=2$ weak matrix elements to the lattice with staggered fermions [9, 11]: the one spin trace formalism and the two spin trace formalism. The four fermion operators can be expressed as products of operators bilinear in the fermion fields. In the one spin trace formalism, each external hadron is contracted with both bilinears of the four fermion operator simultaneously. In the two spin trace formalism, each external hadron is contracted with only one of the bilinears in the four fermion operator. Until now, the two spin trace formalism (2TR) has been used exclusively in calculations of weak matrix elements with staggered fermions, [12, 13, 14], [15, 16, 17, 18] and [19, 20]. Recently, the one spin trace formalism has been developed to a level which permits it to be used for numerical simulations of weak matrix elements such as $B_{K}$ [11]. We have used both formalisms to calculate $B_{K}$ and the results are compared in this paper.

Lattice calculations introduce their own systematic artifacts and errors which can be sizable. One of the principal sources of systematic errors is finite volume. The results of a finite volume comparison were reported in Ref. [16], where it was argued that the finite volume effects on $B_{K}$ are quite small when $V\times T\geq 16^{3}\times 40$ at $\beta=6.0$ ($a^{-1}\sim 2$ GeV), $N_{f}=0$. Another source of systematic errors comes from the quenched approximation, neglecting all of the internal quark loops. Quenched effects were studied in Ref. [15], which concluded that the effect of quenching can not be large. Another systematic error which is called scaling violation comes from the finite lattice spacing ($a\neq 0$). In Ref. [14], the scaling corrections were argued to be of $a^{2}$ order. In other words, $O(a)$ corrections do not exist at all for the staggered fermion operators of $B_{K}$. Related to possible scaling violations is a possible discrepancy between gauge-dependent Landau gauge operators and gauge invariant operators. The Landau gauge operator implies that the quark anti-quark propagators are fixed to Landau gauge and that gauge links between the quark and anti-quark fields are omitted. The question was raised as to whether the Landau gauge operators might cause the large scaling violation which had been noticed originally in Ref. [12]. However, it was reported in Ref. [19] that the results of both Landau gauge operators and the gauge invariant operators were numerically found to be consistent with each other.

The purpose of this paper is to report and analyze the numerical data for $B_{K}$ as well as its individual terms and to interpret the numerical results in terms of various physical models. Part of the preliminary results have already appeared in Ref. [21]. In this paper we will address the following five issues through the interpretation of our numerical results. The first issue is how to select the pseudo Goldstone boson state exclusively. For hadron spectrum measurements, the sink operator picks up a specific hadronic state exclusively. In contrast to hadron spectrum measurements, the operators for weak matrix element measurements do not select a particular hadronic state. We need to impose a symmetry requirement on the wall source such that all the unwanted states are decoupled in the weak matrix element measurement. We have tried an improved wall source method (called cubic wall source [22]) to do this. The second issue is whether the numerical results in the one spin trace form are in agreement with those in the two spin trace form. Theoretically the difference between the two formalisms is supposed to vanish in the limit of $a\rightarrow 0$ [11] for $B_{K}$. We have tried both the one spin trace formalism and the two spin trace formalism to calculate $B_{K}$. The results are compared in this paper. The third issue is the effect of non-degenerate quark antiquark pairs on $B_{K}$ and the individual components making up $B_{K}$. The kaon is composed of $s$ and $d$ valence quarks. Here a non-degenerate quark anti-quark pair implies that the $s$ valence quark mass is different from the $d$ valence quark mass, whereas a degenerate quark anti-quark pair implies that both valence flavors have the same mass. The effect of non-degenerate valence quarks on $B_{K}$ in quenched QCD was mentioned briefly in Ref. [13] where a small difference of only marginal significance was found. In this paper, we investigate in detail the effects of non-degenerate valence quark anti-quark pairs on $B_{K}$ and its individual components. From the theoretical point of view, one effect of non-degenerate valence quark anti-quark pairs can be related to the $\eta^{\prime}$ hairpin diagram in (partially) quenched QCD [23, 24, 25]. The forth issue is whether quenched chiral perturbation theory is compatible with the numerical results of $B_{K}$ [23]. Quenched chiral perturbation theory also predicts the chiral behavior of the individual terms making up $B_{K}$ [23]. It is good to know how reliable these theoretical predictions are numerically. The final issue is whether there is any dynamical fermion effect on $B_{K}$. This question was addressed originally in Ref. [15]. We will re-visit this question and see how important the internal fermion loops are to $B_{K}$. It is important to understand the difference between full QCD and quenched QCD in $B_{K}$ both theoretically and numerically.

This paper is organized as follows. In section 2, we will describe the technical details in doubling the lattice for quark propagators and explain the improved (cubic) wall source method in a self-contained manner. In Section 3, we will specify the lattice operators for $B_{K}$ in brief while leaving the details to adequate references. In Section 4, we describe the parameters for the gauge configurations we generated and the measurement parameters for $B_{K}$. In Section 5, we present the numerical data for ${\cal M}_{K}$ (numerator of $B_{K}$) and the vacuum saturation amplitude (denominator of $B_{K}$). The main emphasis is put on the determination of the plateau region. We also discuss the improved wall source with the wrong flavor channel, which is supposed to vanish in the limit of $a\rightarrow 0$. In Section 6, the numerical results from the data analysis are interpreted. The improved (cubic) wall source method is compared with the conventional method. The data of the one spin trace form and of the two spin trace form are compared. Enhanced chiral logarithms in the individual terms making up $B_{K}$ appears to be seen numerically as well as the additional divergence which arises from the effect of the non-degenerate quark anti-quark pairs. We emphasize the large effect of non-degenerate quark anti-quark pairs on the individual components of $B_{K}$. We compare our data with other groups (Kilcup, Sharpe and Ukawa et al.). We describe the covariance fitting procedure for $B_{K}$ and the negligible effect of non-degenerate quark anti-quark pair on $B_{K}$. Section 7 contains a brief summary and our conclusions.

2 Quark Propagators and Wall Sources

Here, we explain the technical details of the quark propagators and the improved wall source method.

Because the lattice size in the time direction is finite, the pions propagating around lattice in the time direction can in principle contaminate the measurements of $B_{K}$. There have been two proposals for avoiding this contamination. The first proposal is to impose Dirichlet boundary conditions in time on the quark propagator and to place the wall source near the boundary [10, 15, 19]. The second proposal is to double the lattice along the time direction and to use periodic boundary conditions in time for the quark propagator [10, 13]. Dirichlet boundary conditions cause a certain number of time slices near the wall source to be contaminated by reflections off the boundary. The time slices lost due to reflections overlap with those contaminated by $\rho$ mesons. Lattice doubling in the time direction suppresses the backward propagating pions which must travel over the whole lattice size in the time direction before they can contribute to the measurements. In contrast to Dirichlet boundary conditions, lattice doubling does not introduce any reflection from the boundary. Hence, the signal from lattice doubling is much cleaner than that obtained using Dirichlet boundary conditions. Unfortunately, lattice doubling takes twice the computational time since the lattice doubling needs two undoubled quark propagators for a given wall source (forward and backward) instead of the one quark propagator required when using Dirichlet boundary conditions.

We adopt the second solution of doubling the lattice for our numerical simulation of $B_{K}$. Here, we describe this second method in detail. The periodic and anti-periodic quark propagators in time on the undoubled lattice are:

[TABLE]

where $(D+M)$ represents the Euclidean Dirac operator, and $G_{P}(x,y)$ and $G_{A}(x,y)$ represent the Green’s functions with periodic and anti-periodic boundary conditions on the undoubled lattice. The source $h(y)$ is introduced in the numerical simulation as follows: for $0\leq x_{t}<N_{t}$,

[TABLE]

where $\hat{t}$ is the unit vector in time direction and $0\leq y_{t}<N_{t}$ ($N_{t}$ is the undoubled lattice size in the time direction). The $\chi(x)$ field in the above satisfies the following: for $0\leq x_{t}<N_{t}$,

[TABLE]

and for $N_{t}\leq x_{t}<2N_{t}$,

[TABLE]

We use periodic and anti-periodic boundary conditions on the undoubled lattice and take the sum and the difference to obtain a quark propagator which is periodic on the doubled lattice in the time direction. Note that the fermionic field $\chi(x)$ is periodic on the doubled lattice in the time direction. One might ask whether periodic boundary conditions in the time direction make any difference compared to the anti-periodic boundary conditions, which are the physical ones at a finite temperature. The answer is that periodic boundary conditions make no difference as long as the volume is large enough. The physical eigenstates in the confined phase of QCD are hadrons, not quarks nor gluons. In the confined phase, a quark must be confined with other quarks or anti-quarks within a small volume of typical hadronic size ($\approx$ ${\cal O}(1fm)$) before it can acknowledge the existence of the boundary with the volume much larger than the hadronic size.

There have been a number of attempts to enhance the overlap with the lightest particle compared to that of the excited states so that one can see the asymptotic signal ($\exp(-M\mid t\mid)$) at smaller time separations and over the longer plateau [15, 22, 26]. The wall source method has an advantage because it enhances the signal of the hadron propagators with respect to the point source method. For hadron mass spectrum measurements, the sink operator possesses the symmetry of a specific hadronic state. In contrast to hadron spectrum measurements, since the operator in the electro-weak effective Hamiltonian does not select any particular hadronic state by itself, the weak matrix element measurements require the symmetry properties of the wall source to determine the specific hadronic state. There have been two attempts to improve the wall source such that it can exclude contamination from unwanted hadronic states: the even-odd wall source method [17, 26] and the cubic wall source method originally proposed by M. Fukugita, et al. in Ref. [22].

Here, we explain the cubic wall source method in detail. The $B_{K}$ measurements requires that the pseudo-Goldstone Kaon mode should be selected exclusively. Hence, we will restrict our discussion to the cubic wall source operator of pseudo-Goldstone mode. Let us start with definitions and notation. The symbol $\bf\bar{a}$ represents one of the 8 vertices in the unit spatial cube (i.e. ${\bf\bar{a}}\in\{0,1\}^{3}$):

[TABLE]

We define $W_{\bf\bar{a}}(\vec{n})$ as follows:

[TABLE]

The cubic wall source operator for the pseudo-Goldstone pions can then be expressed as follows

[TABLE]

where $\epsilon(x)\equiv(-1)^{x_{1}+x_{2}+x_{3}+x_{4}}$, $a$ and $b$ are color indices and $\chi$ is a staggered fermion field. As an example, let us choose the sink operator to be the bilinear operator with spin structure $S$ and flavor structure $F$:

[TABLE]

$y_{A}=2y+A$, where $y\in Z^{4}$ is the coarse lattice coordinate and $A\in\{0,1\}^{4}$ is the hypercubic coordinate. Then the correlation function is

[TABLE]

where $G(x,y)\equiv\frac{1}{2}[G_{P}(x,y)+G_{A}(x,y)]$ for the domain $0\leq x_{t},y_{t}<2N_{t}$ and

[TABLE]

Here $\psi_{\bf\bar{a}}^{cb}(y_{A})$ is actually what we calculate on the computer using a conjugate gradient method with the wall source set to $W_{\bf\bar{a}}$. For the $f_{K}$ and the vacuum saturation part of $B_{K}$, Eq. (15) is used in our numerical simulation. The idea of the above example (bilinear sink operators) can be extended as a whole to the four fermion operator measurements (for example, $B_{K}$) without loss of generality.

3 Operators Computed

Here, we present a set of lattice operators which describe the same physics as the continuum $\Delta S=2$ operator and which have been used for our numerical calculation of $B_{K}$. In the continuum, $B_{K}$ is defined as

[TABLE]

The numerator of $B_{K}$ is a matrix element of the $\Delta S=2$ four-fermion operator with the neutral $K$ meson (kaon) states. The denominator of $B_{K}$ represents the vacuum saturation approximation of the numerator, which inserts the vacuum state between the two bilinears of the $\Delta S=2$ four-fermion operators. For the denominator, one needs to prescribe a lattice bilinear operator which corresponds to the continuum axial current:

[TABLE]

Here the notations are the same as those adopted by Ref. [9, 10, 11]. The axial current is chosen such that the flavor structure is identical to that of the pseudo Goldstone kaon of the exact $U_{A}(1)$ symmetry in the staggered fermion action.

For the numerator, we need some particular set of lattice four-fermion operators which correspond to the continuum $\Delta S=2$ operator. There are two methods to transcribe the continuum $\Delta S=2$ operator on the lattice with staggered fermions: the one spin trace formalism and the two spin trace formalism [9, 10, 11]. In the continuum, there is no difference between the one spin trace form and the two spin trace form of the $\Delta S=2$ operator, since they are connected by the Fierz transformation. However, in the staggered fermion method, the one spin trace form is different from the two spin trace form due to a pure lattice artifact of 4 degenerate flavors.

In the two spin trace formalism, the four–fermion operator in the numerator in Eq. (17) is transcribed to the lattice [9, 10, 11] as a sum of four operators:

[TABLE]

where V (or A) represents the vector (or axial) spin structure, P represents the pseudoscalar-like flavor structure and the subscript $ab;ba$ (or $aa;bb$) represents the color indices of the quark fields (the details of these notations are described in Ref. [11]). The operators in Eq. (19) have the same chiral behavior in the limit of vanishing quark mass as the continuum $\Delta S=2$ operator [10, 11, 23]. In addition, the matrix elements of ${\cal O}^{Latt}_{2TR}$ show the same leading logarithmic dependence on the renormalization scale as the continuum $\Delta S=2$ operator [11, 30, 31].

In the one spin trace formalism, the four-fermion operator of the numerator in Eq. (17) is transcribed to the lattice as follows:

[TABLE]

where again the details of this notation is given in Ref. [11]. In contrast with the two spin trace formalism, the individual terms in Eq. (20) do not possess the same chiral behavior as the continuum $\Delta S=2$ operator [11]. We must add ${\cal O}^{1TR}_{\rm chiral\ partner}$ in order to preserve the correct continuum chiral behavior [11]. By imposing the correct chiral behavior on ${\cal O}^{Latt}_{1TR}$, we determine the chiral partner operator [11] as follows:

[TABLE]

This forces the resulting operator to respect the continuum chiral behavior. The next question is whether the whole operators including the additional chiral partner operators still have the continuum leading logarithmic behavior. First of all, we need to choose the basis operators such that they belong to the identity representation with respect to the 90*∘* axial rotation group (a subgroup of the exact $U_{A}(1)$ symmetry group). This particular choice of the basis operators guarantees the analytic chiral behavior of the continuum. Second, we find that an eigen-operator (Eq. (20) and Eq. (21)) possesses the same leading logarithmic behavior as the continuum $\Delta S=2$ operator [11]. For $B_{K}$ measurements, the difference between the one spin trace operators and the two spin trace operators vanishes as $a\rightarrow 0$ [11]. We have used both one spin trace and two spin trace operators to calculate $B_{K}$ on the lattice and the numerical results are compared later in this paper.

4 Simulations and Measurement Parameters

In this section we describe the simulation parameters for the gauge configurations and $B_{K}$ measurements. Our old hadron mass calculation [27] on a $16^{3}\times 32$ lattice at $\beta=5.7$ informed us that a longer time dimension allows a more precise fitting of the hadron propagators’ exponential time dependence. For this reason the volume of the configuration was chosen as $16^{3}\times 40$. The coupling constant was $\beta$ = 5.7  ($1/a\approx 2$ GeV). The dynamical effects of two degenerate flavors of staggered fermions with a mass 0.01 were incorporated into the gauge configurations, using the hybrid molecular dynamics R-algorithm [28]. The sea quark mass (0.01) has been fixed through all the measurements even though various valence quark masses were chosen for the $B_{K}$ measurement. The gauge configurations were updated by the hybrid molecular dynamics R-algorithm with molecular dynamics step size 0.0078125 [28] and a trajectory length of 0.5 time units.

Now the measurement parameters for $B_{K}$: Every 60 trajectories $B_{K}$ has been measured. The total number of the gauge configuration samples for $B_{K}$ measurements was 155. We have used both cubic wall source and conventional even-odd wall source methods to create the pseudo Goldstone boson. We used two separate wall sources to create $K^{0}$ and $\bar{K}^{0}$ mesons. The distance between these two separate wall sources was 36 lattice units. For each $B_{K}$ measurement, both wall sources were shifted by 15 lattice units in the time direction from the position used in the previous measurement, while the distance between them was fixed to 36 through all the measurements. The valence (quark, antiquark) mass pairs for the $K$ mesons were (0.01, 0.01), (0.02, 0.02), (0.03, 0.03), (0.004, 0.01), (0.004, 0.02), (0.01, 0.03) and (0.004, 0.05). The quark propagators were calculated, using the lattice doubling method described in the previous section. The stopping condition of the conjugate gradient residual for the quark propagator was set to $1.0\times 10^{-8}$. The quark propagators were gauge fixed to Landau condition with a numerical precision of $a^{4}g^{2}(\partial_{\mu}A_{\mu})^{2}<1.0\times 10^{-7}$.

5 Data Analysis

We now present the numerical results of ${\cal M}_{K}$ (the numerator of $B_{K}$) and the vacuum saturation amplitude ${\cal M}^{V}_{K}$ (the denominator of $B_{K}$) with respect to the lattice Euclidean time for Kaon with a quark anti-quark mass pair = (0.01, 0.01). Meanwhile, we explain how the central plateau region has been selected to determine $B_{K}$. We also discuss the numerical results for unrenormalized (naive, bare) $B_{K}$ with respect to the various quark and anti-quark masses. In addition, we will present the numerical measurement of the wrong flavor channel $((V+A)\otimes S)^{2TR}$ in order to see how much contamination comes from the operator mixing and the excited hadronic states which are supposed to vanish in the limit of $a\rightarrow 0$.

5.1 The Denominator: Vacuum Saturation Amplitude

We define the denominator of $B_{K}$ as the vacuum saturation amplitude:

[TABLE]

Numerical data for ${\cal M}^{V}_{K}$ is presented in Figure 1 (the even-odd wall source method) and in Figure 2 (the cubic wall source method). The data points in Figure 1, 2, and 3 are obtained by single elimination jack-knife method. Each data point in Figure 1(the even-odd wall source) has about twice larger error than that in Figure 2 (the cubic wall source). However, one needs to notice that the cubic wall source method takes four times longer time to compute than the even-odd wall source method.

In Table 1, the results of the covariance fitting of the vacuum saturation together with their $\chi^{2}$ per degree of freedom are given with respect to the various fitting ranges in the lattice Euclidean time for quark mass $m_{q}a=0.01$. From Table 1, for the cubic wall source method the minimum of the $\chi^{2}$ per degree of freedom extends to the fitting range of $11\leq t\leq 24$. For the even-odd wall source method, the minimum of the $\chi^{2}$ per degree of freedom occurs in the fitting range of $13\leq t\leq 23$. It is important to keep in mind the fact that both the one spin trace form and the two spin trace form of $B_{K}$ have the vacuum saturation in common.

5.2 The Numerator: ${\cal M}_{K}$

We define the numerator of $B_{K}$ as:

[TABLE]

Numerical data for ${\cal M}_{K}$ calculated in the two spin trace formalism with the even-odd wall source method is presented in Figure 3. Numerical data for ${\cal M}_{K}$ calculated in the two spin trace formalism with the cubic wall source method is shown in Figure 4. Numerical data for ${\cal M}_{K}$ calculated in the one spin trace formalism with the even-odd wall source method is drawn in Figure 5. All the data points in Figures 3, 4, and 5 are obtained by the single elimination jack-knife method. Each data point in Figure 3 (the even-odd wall source) has an about twice larger error than that in Figure 4 (the cubic wall source). The nearby data points in Figure 3 and Figure 5 (the even-odd wall source) have larger fluctuations than those in Figure 4 (the cubic wall source), which are reflected on $\chi^{2}/\rm d.o.f.$ in Table 2. However, this was not observed in the results for vacuum saturation amplitude. It is important to note the fact that the computational time for the cubic wall source method is four times longer than that for the even-odd wall source method.

In Table 2, the covariance matrix fitting results of ${\cal M}_{K}$ and its $\chi^{2}/(\rm d.o.f.)$ are collected with respect to the fitting range in the lattice Euclidean time. From Table 2 it is difficult to choose the optimal fitting range as all of the $\chi^{2}/(\rm d.o.f.)$ values are within $0.5\sim 1$ $\sigma$ for ${\cal M}_{K}$ calculated in the two spin trace form with the cubic wall source. This is also true for ${\cal M}_{K}$ calculated in the one spin trace form with the even-odd wall source. For these cases, we therefore choose the optimal fitting range consistent with the vacuum saturation amplitude ${\cal M}^{V}_{K}$.

For ${\cal M}_{K}$ calculated in the two spin trace form with the even-odd wall source, we notice that the optimal fitting range is $13\leq t\leq 22$. The ${\cal M}_{K}$ results in the cubic wall source method has more correlation (less fluctuation) between neighboring time slices than those in the even-odd wall source method, while this is not obvious for the measurements of vacuum saturation amplitude ${\cal M}^{V}_{K}$.

5.3 The Ratio: $B_{K}$

In the previous sections, we have discussed the covariance matrix fitting result of ${\cal M}_{K}$ and ${\cal M}^{V}_{K}$ as a function of the fitting range. From this analysis of the fitting ranges, we determine the optimal fitting range for $B_{K}$. The optimal fitting ranges we have chosen to determine $B_{K}$ with respect to various quark mass pairs are summarized in Table 3. Once the fitting range is chosen, $B_{K}$ can be determined through the covariance fitting to a constant. One of our fitting procedures is naive jack-knife analysis of the data (conventional) and the other uses the covariance fitting on the bootstrap ensembles.

The numerical results for $B_{K}$ calculated in the two spin trace form with the even-odd wall source method are shown in Figure 6. The numerical results for $B_{K}$ calculated in the two spin trace form with the cubic wall source method are drawn in Figure 7. The numerical results for $B_{K}$ calculated in the one spin trace form with the even-odd wall source method appear in Figure 8.

Table 4 is a collection of the numerical results of $K$ meson mass and $B_{K}$ with respect to the quark mass pairs calculated in the optimal fitting range by the jack-knife method. The results for $K$ meson mass ($m_{K}$) are obtained by analyzing the results for axial current matrix element with an external $K$ meson state $\langle 0\mid\bar{S}(\gamma_{0}\gamma_{5}\otimes\xi_{5})D\mid K\rangle$, which are also used to obtain the vacuum saturation amplitude. We fitted the logarithm of the axial current data to the linear function of the Euclidean time $A+m_{K}t$ in the optimal fitting range. We would like to thank Pavlos Vranas for checking these mass results using his own fitting program. Table 5 presents the results of $B_{K}$ obtained by the covariance fitting over the optimal fitting range on the bootstrap ensembles.

5.4 The Wrong Flavor Channel

Here we address two important questions on the validity of our approach to $B_{K}$. The first question comes from the fact that the higher-loop radiative correction of the four-fermion composite operators cause the violation of the continuum spin and flavor symmetries. It is important to know non-perturbatively, how large is the contribution of such symmetry violating terms to the weak matrix element measurements for finite lattice spacing. Note, such terms are supposed to vanish in the limit of $a\rightarrow 0$. The smaller the contribution from these wrong flavor channels is, the more reliable our connection between the continuum and lattice observables (either at tree level or at one loop level). The second question is how exclusively we can select the pseudo-Goldstone mode through our improved wall source method. In other words, we want to ask how reliably our wall source technique suppresses the unwanted hadronic states.

We have chosen the $((V+A)\otimes S)^{\rm 2TR}$ operator in order to address the above two questions. The matrix element of this operator with external $K$ mesons is supposed to vanish in the continuum limit ($a\rightarrow 0$) of lattice QCD, due to the vanishing flavor trace. The numerical results for this wrong flavor channel are shown in Figure 9 for the even-odd wall source and Figure 10 for the cubic wall source. From Figures 9 and 10, we notice that the value of the wrong flavor channel is extremely suppressed (less than 1% of $B_{K}$) in both cases. This implies that the unwanted operator mixing of $((V+A)\otimes S)^{\rm 2TR}$ should be at most 1% of $B_{K}$ since it is suppressed by $\alpha_{s}/(4\pi)$ as well as by the vanishing flavor trace. This fact that the unwanted operator mixing is smaller than 1% of $B_{K}$ is of great significance to our approximate matching between the continuum and lattice composite operators. So far, one has neglected those terms of wrong flavor channels which enters at order $g^{2}$, when one connects the lattice $B_{K}$ to the continuum $B_{K}$ at one loop level [10, 11, 19, 30]. The main reason was that these wrong flavor channels will not contribute to $B_{K}$ at all as $a\rightarrow 0$. Hence, the remaining difficulty was to know how large is the contribution of those wrong flavor channels at finite lattice spacing. Our non-perturbative measurements of one of wrong flavor channels shows that the contribution from these wrong flavor channels will be at most 1% of $B_{K}$ and so much smaller than the statical and other systematic errors at $\beta=5.7$ ($a^{-1}=2GeV$). This gives us a great confidence in our approximate matching at one loop level, where such terms of wrong flavor channels, which enters at order $g^{2}$, are neglected.

From Figures 9 and 10, we also observe that the nearest neighboring data points have stronger correlation (less fluctuating) in the cubic wall source than in the even-odd wall source. The four times more floating point computation in the cubic wall source than in the even-odd wall source explains why the error bars of the data in Figure 10 is about half of those in Figure 9. Hence, we conclude that the even-odd wall source is equivalent to the cubic wall source from the standpoint of statistics. We believe that it is better to use the cubic wall source in weaker coupling simulations because the unwanted contamination from the degenerate and excited hadronic states becomes even more severe there.

6 Data Interpretation

Here, we interpret the numerical results in terms of the physics. First of all, we will address the technical questions about the improved wall source methods. Second, we will compare the numerical results of the one spin trace form and the two spin trace form and discuss the meaning of the consistency between the two formalism. Third, we will discuss the chiral behavior of $B_{K}$, and the chiral behavior of the individual components making up $B_{K}$ (i.e. $B_{A}$, $B_{A1}$, $B_{A2}$, $B_{V}$, $B_{V1}$, $B_{V2}$). In addition, we will discuss the effects of non-degenerate quark anti-quark pairs on the chiral behavior of these quantities. We will also discuss the chiral behavior from the standpoint of both $1/N_{c}$ suppression and (partially quenched) chiral perturbation theory. Fourth, we will compare our numerical results with earlier works and discuss the effect of quenching on $B_{K}$ measurements. Finally, we will present our best value of $B_{K}$ in the physical limit as well as our fitting procedure.

6.1 Comparison of Wall Sources

Numerical results for unrenormalized (naive, bare, or tree-level) $B_{K}$ for various average quark mass, calculated both by the even-odd wall source and by the cubic wall source are shown in Figure 11. The values of $B_{K}$ calculated by the two wall sources agree within errors. The error bars for $B_{K}$ calculated by the cubic wall source are about half of those of the even-odd wall source. The cubic wall source method takes four times more computational time than the even-odd wall source method. We conclude that for $B_{K}$ measurements at $\beta=5.7$ the cubic wall source results are consistent with the even-odd wall source results.

6.2 Comparison of One Spin Trace Form

and Two Spin Trace Form

The numerical results for unrenormalized $B_{K}$ in both one spin trace and two spin trace methods are presented in Figure 12. The numerical results for one-loop renormalized (tadpole improved through mean field theory [11, 29, 31, 30, 32, 33]) $B_{K}$ both in the one spin trace and two spin trace methods are shown in Figure 13. The numerical results for one-loop renormalized (RG improved [11, 34]) $B_{K}$ both in the one spin trace and two spin trace form are given in Figure 14. From Figures 12, 13 and 14, we notice that the results for the renormalized $B_{K}$ calculated in the one spin trace method agree with those in the two spin trace method better than those for the unrenormalized $B_{K}$.

Let us explain how we have obtained the renormalized coupling constant for the perturbative expansion. The one-loop renormalization of the $\Delta S=2$ four-fermion operators on the lattice is explained in Ref. [11, 31, 30] in detail. The detailed explanation of matching between the continuum and lattice observables at one loop level is also given in Ref. [11, 31, 30]. We discuss here how to obtain the coupling constant for the renormalization of the composite operators. For perturbative matching at one-loop level, one needs a well-defined coupling constant to use as the perturbative expansion parameter [33, 34]. We have chosen the $\overline{MS}$ coupling constant at $\mu_{\overline{MS}}=\pi/a$ scale as our perturbative expansion parameter [19]. There are two methods to obtain the $\overline{MS}$ coupling constant from the bare lattice coupling constant. One is a non-perturbative approach using tadpole improvement by mean field theory [33] and the other is a purely perturbative approach using renormalization group improvements [34]. We will refer to these two approaches to obtaining $g^{2}_{\overline{MS}}$ as the “tadpole improved” and “RG improved” methods. We use both methods to obtain the $\overline{MS}$ coupling constant at the renormalization scale $\pi/a$.

The details of the non-perturbative approach of tadpole improvement by mean field theory is given in Ref [33, 11, 34, 31, 30]. The $\overline{MS}$ coupling constant ($g^{2}_{\overline{MS}}$) is related to the tadpole-improved coupling constant ($g^{2}_{MF}$) and the bare lattice coupling constant ($g^{2}_{0}$) through mean field theory [11, 34]:

[TABLE]

where $g^{2}_{0}(a)$ is the bare lattice coupling constant, $U_{\Box}$ is a unit gauge link plaquette, and

[TABLE]

For $\beta=5.7$ and $N_{f}=2$, we obtain the tadpole improved $\overline{MS}$ coupling constant:

[TABLE]

The coupling constant in Eq. (27) is used for tadpole improved matching between the continuum and the lattice observables.

The details about a perturbative approach of RG improvement is explained in Ref. [34]. In this approach, the $\overline{MS}$ coupling constant ($g^{2}_{\overline{MS}}$) is related to the bare lattice coupling constant ($g^{2}_{0}(a)$):

[TABLE]

For $\beta=5.7$ and $N_{f}=2$, we obtain the $\overline{MS}$ coupling constant at $\mu_{\overline{MS}}=\frac{\pi}{a}$ scale:

[TABLE]

The coupling constant in Eq. (29) is used for RG improved matching between the continuum and the lattice observables.

6.3 Chiral Logarithms and Non-degenerate Quark Mass Pairs

There have been a number of theoretical attempts to understand the chiral behavior of the $K$ meson $B$ parameters in terms of chiral perturbation theory [23, 35, 36, 37]. In order to discuss the chiral behavior of the $B$ parameters in an organized way, we need to consider a theory with four valence flavors: $S$ and $S^{\prime}$ both with mass $m_{s}$ as well as $D$ and $D^{\prime}$, both with mass $m_{d}$. Let the $K^{0}$ be the $\bar{S}\gamma_{5}D$ pion and the $K^{\prime 0}$ be the corresponding state with primed quarks (i.e. $\bar{S}^{\prime}\gamma_{5}D^{\prime}$ pion). Let us define $B_{V1}$, $B_{V2}$, $B_{A1}$ and $B_{A2}$ as:

[TABLE]

where $a$ and $b$ are color indices. The $B_{A}$, $B_{V}$ and $B_{K}$ are expressed in terms of $B_{V1}$, $B_{V2}$, $B_{A1}$ and $B_{A2}$:

[TABLE]

Let us discuss the chiral behavior of the B parameters defined above ($B_{A1}$, $B_{A2}$, $B_{V1}$, $B_{V2}$). In the vacuum saturation approximation with $1/N_{C}$ suppression, $B_{A1}$ = 0.25, $B_{A2}$ = 0.75 and $B_{V1}=B_{V2}$ = 0 (obviously wrong!). For $m_{u}=m_{d}\neq m_{s}$, chiral perturbation theory predicts the following results [10, 23]:

[TABLE]

where in the vacuum saturation approximation (tree-level chiral perturbation), $\alpha_{+}=1$, $\alpha_{-}=\frac{1}{3}$, $\delta_{1}=\frac{1}{3}$, $\delta_{2}=1$, $\beta_{1}=1$, $\beta_{2}=\frac{1}{3}$ while $\gamma_{1}$ and $\gamma_{2}$ is not determined. $A$, $C$ and $D$ are defined as

[TABLE]

where

[TABLE]

$\Lambda$ is introduced as a momentum cut-off regularization for chiral perturbation theory. $m_{ss}^{2}=2\mu m_{s}$, $m_{dd}^{2}=2\mu m_{d}$ and $m_{K}^{2}=\mu(m_{s}+m_{d})$ (i.e. $2m_{K}^{2}=m_{ss}^{2}+m_{dd}^{2}$). Note there is no quadratic divergence in $B_{V1}$, $B_{V2}$, $B_{A1}$, and $B_{A2}$ even though $I_{1}(m)$ has a quadratic divergence. The quadratic divergences in $C$ and $D$ cancel out. The quadratic divergences in $A$ is not a function of quark masses $m_{s}$ or $m_{d}$. Hence, these quadratic divergences in $A$ can be absorbed into the coefficients $\alpha_{-}$ and $\alpha_{+}$. As a summary of chiral perturbation theory, let us choose $\Lambda=4\pi f$ and rewrite the results for $B_{V1}$, $B_{V2}$, $B_{A1}$, and $B_{A2}$ as follows:

[TABLE]

where the constant terms proportional to $\beta_{i}$ are absorbed into $\delta_{i}$, and

[TABLE]

Here we keep full functional form with respect to $\epsilon$ for theose terms of order $\ln(z)$ and $z^{0}$. For those terms of order $z\ln(z)$ and $z$, we do Taylor expansion with respect to $\epsilon$ and keep only the terms of order $\epsilon^{0}$ and $\epsilon^{2}$. The coefficients $\sigma_{i}$, $\eta_{i}$, $\xi_{i}$ and $\rho_{i}$ ($i\in\{V1,V2,A1,A2\}$) are unknown. All the results in Eqn. (37-40) and Eqn. (46-49) are calculated using full QCD as the fundamental theory. Note, $B_{V1}$, $B_{V2}$, $B_{A1}$, and $B_{A2}$ have a branch point at $\epsilon=\pm 1$ which is non-singular. Let us set the domain of $\epsilon$ to $-1\leq\epsilon\leq 1$ and $\epsilon\in R$. $B_{V1}$, $B_{V2}$, $B_{A1}$, and $B_{A2}$ have no singularity on this physical domain of $\epsilon$. Also note, $B_{V1}$, $B_{V2}$, $B_{A1}$, and $B_{A2}$ are even functions of $\epsilon$ which means that although we switch $m_{s}$ with $m_{d}$, the physics does not know it at all (i.e. it does not change). Hence, the power series expansion of $B_{V1}$, $B_{V2}$, $B_{A1}$, and $B_{A2}$ with respect to $\epsilon$ should have only even powers of $\epsilon$. The one-loop corrections to $B_{V1}$, $B_{V2}$ $B_{A1}$ and $B_{A2}$ include a term of the order $I_{2}(m_{K})/m_{K}^{2}$ which is proportional to $\beta_{i}$ in Eqn. (37-40). $I_{2}(m_{K})/m_{K}^{2}$ is proportional to $\ln(m_{K})$ which is logarithmically divergent as $m_{K}\rightarrow 0$. These logarithmically divergent terms are called enhanced chiral logarithms [23], and were originally noticed by Langacker and Pagels [38]. It is important to note that in full QCD, the enhanced chiral logarithms in Eqn. (37-40) and Eqn. (46-49) are all proportional to $\beta_{i}\frac{\mu^{2}}{m_{K}^{2}}I_{2}(m_{K})$, which depends only on the average mass of the quark antiquark pair.

Let us adapt the above results to the (partially) quenched limit of QCD [23]. Here, quenched approximation means neglecting all the internal fermion loops and keeping only pure gauge interactions (i.e. sea quark mass is infinitely heavy), while the partially quenched approximation implies that the sea quark mass of the internal fermion loops is different from the valence quark mass. Both partially quenched and quenched approximations have additional infra-red problems which are absent in full QCD, since the sea quark fermion determinant can not regulate the infra-red pole singularity of the valence quark propagator. There are two important differences between (partially) quenched and full QCD. The first difference is that the meson eigenstates are not the same. The second difference comes from $\eta^{\prime}$ loops, which can not contribute in full QCD simply because the $\eta^{\prime}$ is too heavy. The diagrams of the hairpin type ($\eta^{\prime}$ loops) present in the (partially) quenched chiral perturbation for $B_{K}$ vanish in the limit of $m_{s}=m_{d}$ [23]. The logarithmically divergent terms which come from the hairpin diagrams are called quenched chiral logarithms [24, 25, 23]. Here, the key point is that in (partially) quenched QCD, the quenched chiral logarithms in the B parameters, if present, must be a function of both average quark mass ($m_{s}+m_{d}$) and mass difference ($m_{s}-m_{d}$) of the quark anti-quark pair, while the enhanced chiral logarithms common in both full and (partially) quenched QCD are not a function of quark mass difference but a function of only average quark mass.

The quantity $B_{A1}$ in Eqn. (39) has a finite constant term which is a factor of around 1/3 smaller than that in $B_{A2}$ in Eqn. (40), while the enhanced logarithmic term in $B_{A1}$ is around 3 times larger than that in $B_{A2}$. The enhanced logarithmic term in $B_{V1}$ is also around 3 times larger than that in $B_{V2}$. Therefore, we have chosen $B_{A1}$ as a useful measurement adequate to observe both the enhanced chiral logarithms and, if present, the (partially) quenched chiral logarithms.

We plot $B_{A1}$ and $B_{A2}$ with respect to the average quark masses in Figure 15 and $B_{V1}$ and $B_{V2}$ with respect to the quark masses in Figure 16. One can see a difference in the chiral behavior of $B_{A1}$ and $B_{V1}$ between the case with degenerate quark anti-quark pairs of mass: {0.01, 0.02, 0.03} and the situation with non-degenerate quark anti-quark pairs with the mass pairs: {(0.004, 0.01), (0.004, 0.02), (0.004, 0.05), (0.01, 0.03)}. In order to more precisely see the effect of non-degenerate quark anti-quark pairs on $B_{A1}$, we must interpolate in quark mass. Toward this end, we fit the data of the degenerate quark anti-quark pairs to the following function:

[TABLE]

and, second, subtract this form for the degenerate case from the non-degenerate data as follows:

[TABLE]

where $B_{A1}(m_{K},m_{s}-m_{d})$ is the numerical data for the non-degenerate quark mass pairs.

The covariance fitting results of the degenerate data on the jack-knifed ensembles are is

[TABLE]

Using chiral perturbation theory, one can estimate the value of $A_{1}$ in Eq. (52). From Eq. (48),

[TABLE]

where $\mu\equiv m_{K}^{2}/(m_{s}+m_{d})$. Using the results of our lattice QCD simulation ($\mu\cong 2.3a^{-1}$), this amplitude $A_{1}$ is written in terms of a poorly known quantity $\beta_{1}\sim 1$: $A_{1}\cong-11.6\beta_{1}$. Using the value of $K$ meson mass and strange quark mass in the particle data book ($m_{s}\cong 100\sim 300MeV$ and $m_{s}/m_{d}\cong 25$), one can also obtain $A_{1}\cong-(0.5\sim 3.5)\beta_{1}$ in terms of a poorly known quantity $\beta_{1}$, where the range of values is chosen to reflect the large uncertainty in the chiral perturbation theory prediction.

The numerical data for $\Delta B_{A1}(m_{K})$ with respect to the average quark mass are plotted in Figure 17. From Figure 17, it is obvious that there is a strong sensitivity to the non-degenerate quark mass pairs (i.e. $\Delta B_{A1}(m_{K})$ is not only a function of $(m_{s}+m_{d})$ but also a function of $(m_{s}-m_{d})$). This additional divergence may be related to (partially) quenched chiral logarithms. It could also come from finite volume dependence on the lightest mass of the quark anti-quark pair. Something else might cause this unexpected divergence. The key point is that there is an additional divergence which is visible in our numerical simulations. At any rate, we would like to raise the following questions:

• •

$\eta^{\prime}$ loop: Could the additional hairpin diagram in partially quenched chiral perturbation theory explain quantitatively the additional divergence which is sensitive to the non-degenerate quark masses? In fact, it is known that the contribution from $\eta^{\prime}$ loops vanishes as $m_{s}\rightarrow m_{d}$ [23, 10].

• •

finite volume effect: Could the finite volume effect of $B_{A1}$ or $B_{V1}$ be sensitive to the non-degenerate quark masses? The small eigenvalues and their density in lattice QCD is regulated by the finite volume. Could these small eigenvalues be set up such that they are sensitive to the lightest mass of the non-degenerate quark anti-quark pair?

• •

fermion determinant in partially quenched QCD: In partially quenched QCD, the sea quark mass is different from the valence quark mass. In partially quenched QCD, the fermion determinant of the sea quark may suppress the coefficient of the logarithmic divergence more efficiently than the case of quenched QCD, which has no fermion determinants at all. Could one see the much larger effect of the non-degenerate quark anti-quark pair on $B_{A1}$ and $B_{V1}$ in quenched QCD than in partially quenched QCD?

• •

scaling violation: The anomaly current in staggered fermion formulation $\bar{\chi}\overline{(\gamma_{\mu}\gamma_{5}\otimes I)}\chi$ is not a conserved current for finite lattice spacing. Hence, the corresponding pseudo-Goldstone pion $\bar{\chi}\overline{(\gamma_{5}\otimes I)}\chi$ has a serious contamination of finite lattice spacing, which is supposed to vanish by flavor symmetry restoration in the continuum limit of $a\rightarrow 0$. This suggests that even in quenched QCD, $\eta^{\prime}$ in staggered fermion formulation for finite lattice spacing may be much heavier than pseudo Goldstone pion $\bar{\chi}\overline{(\gamma_{5}\otimes\xi_{5})}\chi$ and that its mass may have a scaling violation term of order $a$ to make matters worse. How large is the scaling violation of $\eta^{\prime}$ mass as a function of quark mass?

• •

other possibility: Is there something else related to the systematics of the non-degenerate quark masses on the lattice?

One important thing is that this dependence on quark mass difference amplifies the enchanced chiral logarithms rather than reducing them. We see this in Figure 15, 16 and 17.

$B_{A}$ and $B_{V}$ are plotted in Figure 18, together with $B_{K}$ with respect to the average quark masses. One can notice that all these divergences related to both (partially) quenched chiral logarithms (if present) and enhanced chiral logarithms (present in $B_{A1}$, $B_{A2}$, $B_{V1}$ and $B_{V2}$) are canceled out in $B_{K}$ which is finite in the chiral limit.

6.4 Comparison with Earlier Work

We now compare our numerical results of $B_{K}$ with those of other groups as well as comparing our results at $\beta=5.7$ (full QCD with $N_{f}=2$) with the results at $\beta=6.0$ (Quenched QCD).

There have been two groups to calculate $B_{K}$ at $\beta=5.7$ with staggered fermions in full QCD. Kilcup [15] has calculated the unrenormalized (naive) $B_{K}$ at $\beta=5.7$ ($16^{3}\times 32$, full QCD with two dynamical flavors of a mass $m_{\rm sea}a$ = 0.01, 0.015, 0.025). The lattice scale was $a^{-1}=1.9\sim 2.0$ GeV. The number of independent configurations was 50 and $B_{K}$ measurements were done twice in different locations on the lattice for each individual configuration, to make the total number of configurations equivalent to 100. Quark wall propagators were fixed in Landau gauge, and periodic boundary conditions in space and Dirichlet boundary conditions in the time direction were imposed.

Ukawa et al. [19, 20] studied $B_{K}$ at $\beta=5.7$ on a lattice of size $20^{3}\times 20$ (duplicated in the time direction) with two flavors of dynamical staggered quarks of mass $m_{\rm sea}a=0.01$. The lattice scale was $a=0.085\sim 0.09$ fm ($a^{-1}=2.2\sim 2.4$ GeV). Both Landau gauge operators and gauge-invariant operators were used. Quark propagators were calculated with Dirichlet (periodic) boundary condition in the time (space) direction.

The differences between our numerical simulation and those of other groups’ are the lattice size, the boundary conditions on the quark propagators in the time direction and the color summation over the meson wall sources. We are not summing the three values of the color index for the meson wall sources on the individual configuration sample. Instead, we choose a different color index for the meson wall sources in each measurement (in other words, color indices of meson wall sources are spread over the configuration samples with equal statistical weight instead of being summed on each configuration sample). The $B_{K}$ results of this work, Kilcup and Ukawa et al. are summarized in Table 6. The $B_{A1}$, $B_{A2}$, $B_{V1}$ and $B_{V2}$ of this work and Kilcup [39] are compared in Tables 8 and 8. We believe that the differences in Tables 6, 8 and 8 can be reasonably explained by the different lattice size, the boundary conditions on the quark propagators, color summation of the meson wall source, poor statistics, the uncertainties in the lattice spacing, etc.. Hence, we conclude that all the measurements in Tables 6, 8 and 8 are consistent with one another.

It is also important to compare the results of full QCD ($N_{f}=2$) with those of quenched QCD. There have been three groups to calculate $B_{K}$ at $\beta=6.0$ ($a^{-1}=2.0\sim 2.1$ GeV) in quenched QCD. This corresponds to $\beta=5.7$ in full QCD.

The results of this work (full QCD, $\beta=5.7$, $16^{3}\times 40$), Sharpe et al. [13, 18] (quenched QCD, $\beta=6.0$, $24^{3}\times 40$), Kilcup [15, 18] (quenched QCD, $\beta=6.0$, $16^{3}\times 40$) and Ukawa et al. [19, 20] (quenched QCD, $\beta=6.0$, $24^{3}\times 40$) are compared in Table 9. From Table 9, we conclude that the dependence of $B_{K}$ on $m_{\rm sea}$ is too weak to detect in our numerical simulations (in other words, the effect of quenched approximation is less than 15 % in $B_{K}$ measurements). This has been predicted by the quenched chiral perturbation theory in the limit of $m_{s}\rightarrow m_{d}$ [10, 23].

6.5 Fitting Procedures for $B_{K}$

$B_{K}$ describes $K^{0}-\bar{K}^{0}$ mixing at the energy scale of about 500 MeV (i.e. in the low energy limit of QCD dynamics). It is not known how to calculate the dependence of $B_{K}$ on the valence quark mass directly from the QCD Lagrangian. For this reason we adopt the chiral effective Lagrangian (equivalent to current algebra), which is valid in the energy region below the $\rho$ meson mass. This chiral effective lagrangian is not a cure-all solution to the low energy dynamics of QCD. However, it gives us a reasonable guide to understand the leading chiral behavior of QCD. We will now use the predictions of the chiral effective Lagrangian to interpolate between our $B_{K}$ results to make a $B_{K}$ prediction for physical quark masses.

The corrections from chiral perturbation theory to ${\cal M}_{K}$ in full QCD were calculated in Ref. [23, 35, 36, 37]. The results of that calculation [23] ($m_{u}=m_{d}\neq m_{s}$) in full QCD were

[TABLE]

where $\epsilon$ is defined as

[TABLE]

$c_{1}$, $c_{2}$ are unknown coefficients, but can be determined by numerical simulation on the lattice. The results from quenched chiral perturbation theory [23] ($m_{u}=m_{d}\neq m_{s}$) are

[TABLE]

where quenched chiral perturbation [23] predicts

[TABLE]

$c^{\prime}_{1}$ and $c^{\prime}_{2}$ are unknown coefficients. The term proportional to $\delta$ is a contribution of $\eta^{\prime}$ loops appearing in quenched QCD, which are absent in full QCD. Because of this term, $B_{K}^{\rm quenched}$ has a singular branch point at $\epsilon=\pm 1$. Note, $B_{K}^{\rm full\ QCD}$ does not have any singular branch point. Equations (57) and (59) are the theoretical predictions for $B_{K}$ as a function of the light quark masses.

Let us choose our fitting function on the basis of the predictions of (quenched) chiral perturbation theory in Eq. (57) and (59). Parts of our numerical simulation of $B_{K}$ are classified in the category of partially quenched QCD ($N_{f}=2$ but $m_{\rm sea}\neq m_{\rm valence}$). We notice that the coefficients of those terms proportional to $\epsilon^{2}m_{K}^{2}\ln(m_{K}^{2})$ in Eqn. (57) and Eqn. (59) are different from each other, which implies that the coefficient of these terms should be determined by our numerical data.

We have tried a linear fitting function ($B_{K}(m_{q})=\alpha_{1}+\alpha_{2}m_{q}a$) even though the linear term is a next to leading order correction from the chiral perturbation (as can be seen in Eqn. (57) and (59)). The $\chi^{2}$/(d.o.f.) for the linear covariance fitting on the jack-knifed ensembles of $B_{K}$ data calculated in the two spin trace form with the cubic wall source method is 17.2(37). This implies that the fitting is poor and that we need an additional term (e.g. $m_{q}a\ln(m_{q}a)$) to fit the data.

Hence, we choose the first fitting function as

[TABLE]

where $\alpha_{i}$ ($i=1,2,3$) are unknown coefficients to be determined and $m_{q}a=\frac{1}{2}(m_{s}+m_{d})$. In this first fitting function, we neglected the effect of non-degenerate quark mass pairs (i.e. terms proportional to $\epsilon^{2}$ are neglected). From the chiral perturbation theory Eqn. (57) and (59), the predicted value of the ratio $\alpha_{2}/\alpha_{1}$ is

[TABLE]

where $0\leq\Delta\leq 1$. Using our lattice QCD simulation results for $\mu$ ($\mu\cong 2.3a^{-1}$), the ratio $\alpha_{2}/\alpha_{1}$ is written in terms of poorly known quantity $\Delta$:

[TABLE]

Using the value of $K$ meson mass and strange quark mass in the particle data book ($m_{s}\cong 100\sim 300MeV$ and $m_{s}/m_{d}\cong 25$), one can also obtain the ratio $\alpha_{2}/\alpha_{1}$

[TABLE]

in terms of a poorly known quantity $\Delta$, where the range of values is chosen to reflect the large uncertainty in the chiral perturbation theory prediction. The $B_{K}$ results of the above 3 parameter covariance fitting on the bootstrap ensembles are summarized in Table 11 (unrenormalized) and 11 (tadpole-improved renormalized).

From Table 11 and 11, we observe that the renormalization with tadpole improvements makes $B_{K}(m_{q})$ in agreement between the one spin trace form and the two spin trace form once they are obtained with the same source method and covariance-fitted in the same range. The covariance fitting results for the ratio $\alpha_{2}/\alpha_{1}$ were not consistent among the various measurements, mainly because the 3 parameters have a wide domain to fit the 7 $B_{K}$ data points as a function of the average quark mass. This gave us a motivation for the second fitting trial function, which will have only two parameters.

Our second fitting trial function is chosen such that the ratio $\alpha_{1}/\alpha_{2}$ is fixed to the case of the degenerate quark anti-quark pair ($\epsilon=0$). The reason is that this ratio is universal for the degenerate quark pair regardless of quenched approximation and that our $\epsilon$ values are small enough to take into consideration later as a perturbative expansion parameter. The second fitting trial function is

[TABLE]

where $\alpha_{i}$ ($i=1,2$) are unknown coefficients. We take the effective mass ($m_{K}$) of $f_{K}$ ($\langle 0\mid A_{\mu}\mid K^{0}\rangle$) measurements in Table 4 as a definition of $m_{K}$ in Eq. (65). In this second fitting trial function, we set the coefficient of the leading term in the chiral perturbation expansion to the theoretically expected value in Eq. (57, 59) and also we again neglect the effect of non-degenerate quark mass pairs. When we set $a^{-1}$ = 1.8 GeV from the $\rho$ meson mass and choose $f_{\pi}$ = 93 MeV, the covariance fitting results for the two spin trace form (cubic wall source) are summarized in Table 13 and 13.

From Table 13 and 13, we notice that the fitting results for the tadpole-improved renormalized $B_{K}$ in various types of the measurements are in good agreement with one another, while the fitting results for the unrenormalized $B_{K}$ are not consistent between the one spin trace form and the two spin trace form.

Let us discuss how we can detect the effect of non-degenerate quark anti-quark pairs on $B_{K}$ (i.e. the dependence of $B_{K}$ on $\epsilon$). The strategy is the following. First, we divide the numerical data in two parts: one part (we call this the degenerate part) contains only the numerical data for degenerate quark mass pairs {(0.01, 0.01), (0.02, 0.02), (0.03, 0.03)} and the other part (the non-degenerate part) contains only the data for non-degenerate quark mass pairs {(0.004, 0.01), (0.004, 0.02), (0.004, 0.05), (0.01, 0.03)}. Next, we fit the degenerate part to the second trial function in Eqn. (65) which are supposed to be exact for the degenerate quark mass pairs up to the given order in the chiral perturbative expansion. For the degenerate part,

[TABLE]

Now we define a function which represents the effect of non-degenerate quark mass pairs:

[TABLE]

where $B_{K}(m_{K},\epsilon)$ are our numerical data with non-vanishing $\epsilon^{2}=(m_{s}-m_{d})^{2}/(m_{s}+m_{d})^{2}$ and $B_{K}^{\rm degenerate}(m_{K})$ represents the fitting function for the degenerate part, given in Eq. (66). Finally, we try to find a functional form, if possible, to fit $\Delta B_{K}(m_{K})$ numerical data to. In Figure 19 (two spin trace form, cubic wall source), Figure 20 (two spin trace form, even-odd wall source) and Figure 21 (one spin trace form, even-odd wall source), we plot $\Delta B_{K}$ with respect to the average quark mass, only for the non-degenerate quark anti-quark pairs. As you can see in the figures, it is hard to find a functional form which can explain all of the data. From these figures, we notice that the dependence of $B_{K}$ on $\epsilon^{2}$ is extremely small (especially in the domain near the physical $K$ meson mass). Hence, we conclude that we could not detect any significant effect of non-degenerate quark mass pairs on $B_{K}$ within the precision of our numerical study.

As a conclusion to this section, let us present our best value of $B_{K}$. Since the data of the cubic wall source method has better statistics than that of the even-odd wall source method, we shall quote the covariance fitting results of the cubic wall source data as our best value. For the second trial function, the uncertainty in the lattice scale $a$ and the decay constant $f_{\pi}$ ($f_{K}$) produces uncertainty in $B_{K}(m_{K})$ of the same order of magnitude as the statistical uncertainty, while covariance fitting to the first trial function is not as much sensitive to the lattice scale $a$ and the decay constant $f_{\pi}$. Hence, we choose the results of the covariance fitting to the first fitting trial function with three parameters on the bootstrap ensembles as our best value. Interpolated from the first fitting trial function, the physical results are

[TABLE]

where N.D.R. implies naive dimensional regularization scheme and the errors represent purely statistical uncertainty. Here the physical results imply the $B_{K}$ values for the physical kaon mass ($m_{K}=497.7$ MeV). Here, we have completely neglected the systematic errors related to the scale ($a$) uncertainty, the coupling ($g^{2}_{\overline{\rm MS}}$) uncertainty, the contamination from the unwanted operator mixing, and the contamination from unwanted hadronic states which can couple to the operators used. We also could not control finite volume effects, finite temperature effects, or finite lattice spacing effects.

7 Conclusion

Here, we summarize what we have learned through the numerical simulation of $B_{K}$ and what needs further investigation in the future.

The results for $B_{K}$ from the improved wall source (cubic wall source) were in good agreement with those of the conventional even-odd wall source. It is shown that the cubic wall source suppresses the contamination from the wrong flavor channels efficiently. However, the cubic wall source takes four times more computational time than the even-odd wall source. In the limit of $a\rightarrow 0$, the SU(4) flavor symmetry is supposed to be recovered and so there will be a serious contamination from pseudo Goldstone pions with various flavor structures and the various $\rho$ mesons. The cubic wall source is quite promising in the weak coupling region to exclude the contamination from unwanted hadronic eigenstates.

We can transcribe the continuum $\Delta S=2$ operator to the lattice with staggered fermions in two different ways. Theoretically, both formalisms of operator transcription must be equivalent to each other in the limit $a\rightarrow 0$. The numerical results in the one spin trace formalism were consistent with those in the two spin trace formalism after the proper renormalization (with either tadpole or RG ($\overline{MS}$) improvement). We have learned how important the proper renormalization is, as well as the careful choice of the coupling constant for the perturbative expansion.

We have tried to understand the effects of the non-degenerate quark masses on $B_{K}$ and the individual components making up $B_{K}$. The effects of the non-degenerate quark antiquark pairs on $B_{K}$ were too small to observe within the precision of our numerical simulation (especially near the region of physical $K$ meson mass). Why this effect is so small needs more careful theoretical investigation. Chiral perturbation theory suggests that $B_{A1}$ (axial part of $B_{K}$ with one color loop) is the best observable to detect the enhanced chiral logarithms which are not expected to be a function of quark mass difference. We observed an additional divergence which depends on the quark mass difference. This additional divergence needs more thorough investigation and understanding. We wonder whether partially quenched chiral perturbation can explain this additional divergence, or how much finite volume effects on $B_{A1}$ depend on the lightest mass of the non-degenerate quark anti-quark pair. Qualitatively, chiral perturbation theory is consistent with our numerical work. However, more theoretical research on the (partially) quenched chiral perturbation and its quenched chiral logarithms is necessary to see whether the (partially) quenched chiral perturbation can explain the effects of the non-degenerate quark antiquark pairs on $B_{K}$ and its individual components.

We could not observe any dynamical fermion effect on $B_{K}$. It is difficult to understand why internal fermion loops are so unimportant for $B_{K}$, since the Dirac eigenvalue spectrum of quarks is supposed to be completely different between quenched QCD and full QCD. This also needs further theoretical understanding.

Through the numerical study of $B_{K}$ in this paper, we have learned that the cubic source method is promising for the weak coupling simulation and that the one spin trace formalism is consistent with the two spin trace formalism. It is true that lattice QCD results for $B_{K}$ are more solid and believable after this work. However, it is also true that there are many details which need more thorough investigation and understanding.

8 Acknowledgement

One of the authors (W. Lee) is indebted a lot to Norman H. Christ. This work could not have been done without his consistent help and encouragement. Helpful discussion with Robert D. Mawhinney at the early stage of this work is acknowledged with gratitude. One of the authors (W. Lee) would like to express his sincere gratitude to Donald Weingarten at IBM T.J. Watson Research Center for his superb lecture of bootstrap analysis. One of the authors (W. Lee) would like to thank Gregory Kilcup for his kind help in many ways during his visit to Columbia University. Helpful discussion with Stephen Sharpe and A. Ukawa during Santa Fe workshop is acknowledged with sincere gratitude. We would like to thank Zhihua Dong for his kind help in gauge fixing programming. Helpful conversation with Shailesh Chandrasekaharan, Dong Chen, and Decai Zhu is acknowledged with gratitude.