Improved lattice computation of proton decay matrix elements

TL;DR

This paper reports an improved lattice QCD calculation of proton decay matrix elements, achieving higher statistical accuracy and refining key constants, which impacts predictions of proton lifetime and dark matter models.

Contribution

The study introduces an error reduction technique in lattice computations, significantly improving the precision of proton decay matrix elements and related low energy constants.

Findings

Matrix element accuracy improved to 10-15% from 30-40%.

Form factors are 1.4 times smaller than previous estimates using BChPT.

Proton lifetime estimates increase by a factor of 2-3.

Abstract

We present an improved result of lattice computation of the proton decay matrix elements in $N_f=2+1$ QCD. In this study, the significant improvement of statistical accuracy by adopting the error reduction technique of All-mode-averaging, is achieved for relevant form factor to proton (and also neutron) decay on the gauge ensemble of $N_f=2+1$ domain-wall fermions in $m_\pi=0.34$--0.69 GeV on 2.7~fm$^3$ lattice as used in our previous work \cite{Aoki:2013yxa}. We improve total accuracy of matrix elements to 10--15\% from 30--40\% for $p\rightarrow\pi e^+$ or from 20--40\% for $p\rightarrow K \bar\nu$. The accuracy of the low energy constants $\alpha$ and $\beta$ in the leading-order baryon chiral perturbation theory (BChPT) of proton decay are also improved. The relevant form factors of $p\rightarrow \pi$ estimated through the "direct" lattice calculation from three-point function appear to be 1.4 times smaller than those from the "indirect" method using BChPT with $\alpha$ and $\beta$. It turns out that the utilization of our result will provide a factor 2--3 larger proton partial lifetime than that obtained using BChPT. We also discuss the use of these parameters in a dark matter model.

Tables

Table 1. Table 1: Lattice ensemble set and parameters. m 𝑚 m refers to the domain-wall fermion mass for the degenerate light quarks ( u 𝑢 u and d 𝑑 d ). N g subscript 𝑁 𝑔 N_{g} is the number of approximate in both light and strange quark propagators. N eig subscript 𝑁 eig N_{\rm eig} denotes the number of low-mode used in light quark propagator, and “res” is value of squared norm of residual vector for the sloppy solver in AMA. Values of hadron masses are measured with extended quark source and sink with gauge invariant Gaussian smearing.

$a^{-1}$ GeV	$m$	$m_{\pi }$ (GeV)	$m_K$ (GeV)	$m_N$ (GeV)	$N_{\mathrm{eig}}$	$N_g$	res	$N_{\mathrm{conf}}$
1.7848(6)	0.005	0.340(1)	0.594(2)	1.179(5)	300	32	0.003	91
	0.01	0.427(1)	0.626(1)	1.269(5)	300	32	0.003	55
	0.02	0.574(1)	0.688(1)	1.452(4)	200	32	0.003	39
	0.03	0.694(1)	0.744(2)	1.598(5)	200	32	0.003	44

Table 2. Table 2: Error budget of α 𝛼 \alpha and β 𝛽 \beta at μ = 2 𝜇 2 \mu=2 GeV. “ χ 𝜒 \chi ” column is systematic uncertainty due to chiral extrapolation.” a 2 superscript 𝑎 2 a^{2} ” column denotes the systematic uncertainty due to 𝒪 ​ ( a 2 ) 𝒪 superscript 𝑎 2 \mathcal{O}(a^{2}) . “ Δ Z subscript Δ 𝑍 \Delta_{Z} ” and “ Δ a subscript Δ 𝑎 \Delta_{a} ” column are systematic uncertainties from renormalization factor and lattice spacing.

LECs	statistical error	systematic error
		$\chi$	$a^2$	${\Delta }_Z$	${\Delta }_a$
$\alpha$(GeV3) = $-$0.0144(15)	0.0003	0.0005	0.0007	0.0012	0.0002
$\beta$(GeV3) = 0.0144(15)	0.0004	0.0005	0.0007	0.0012	0.0002

Table 3. Table 3: Relative error of systematic uncertainty in chiral extrapolation estimated from three fitting ranges; “light” is fitting range without the heaviest point, “heavy” is fitting range without the lightest point. “total” is total one in quadrature. For reference, we also show the value of chi-squared per degree-of-freedom in our fitting as in “ χ dof 2 subscript superscript 𝜒 2 dof \chi^{2}_{\rm dof} ” column.

Matrix element	Relative error in chiral extrapolation
	total	${\chi }_{\mathrm{dof}}^2$	light	${\chi }_{\mathrm{dof}}^2$	heavy	${\chi }_{\mathrm{dof}}^2$	$𝒪\left(q^4\right)$	${\chi }_{\mathrm{dof}}^2$	$𝒪\left(m{q}^2\right)$	${\chi }_{\mathrm{dof}}^2$
$⟨{\pi }^0\left|{\left(ud\right)}_R{u}_L\right|p⟩$	1.8%	0.6	1.6%	0.8	0.8%	0.8	0.7%	0.6	0.3%	0.2
$⟨{\pi }^0\left|{\left(ud\right)}_L{u}_L\right|p⟩$	5.7%	1.4	3.8%	2.0	4.3%	1.2	2.3%	2.2	2.6%	1.9
$⟨K^0\left|{\left(us\right)}_R{u}_L\right|p⟩$	2.8%	1.4	2.7%	1.7	0.7%	1.5	0.7%	1.6	1.1%	1.4
$⟨K^0\left|{\left(us\right)}_L{u}_L\right|p⟩$	3.1%	1.7	0.8%	1.9	3.0%	1.7	1.0%	2.0	2.1%	0.2
$⟨K^{+}\left|{\left(us\right)}_R{d}_L\right|p⟩$	3.5%	1.3	3.4%	1.5	1.0%	1.5	1.6%	1.3	2.0%	0.8
$⟨K^{+}\left|{\left(us\right)}_L{d}_L\right|p⟩$	7.5%	1.6	2.3%	2.2	7.2%	1.5	3.3%	2.1	1.9%	2.7
$⟨K^{+}\left|{\left(ud\right)}_R{s}_L\right|p⟩$	1.6%	0.9	1.0%	1.2	1.2%	1.1	1.3%	0.8	1.3%	0.1
$⟨K^{+}\left|{\left(ud\right)}_L{s}_L\right|p⟩$	3.9%	1.7	2.1%	2.4	3.3%	1.6	1.4%	2.2	1.5%	1.7
$⟨K^{+}\left|{\left(ds\right)}_R{u}_L\right|p⟩$	2.7%	1.0	2.3%	0.8	1.4%	1.1	2.3%	1.0	0.7%	0.8
$⟨K^{+}\left|{\left(ds\right)}_L{u}_L\right|p⟩$	2.1%	1.8	1.5%	2.4	1.4%	1.8	0.8%	2.2	1.6%	0.8
$⟨\eta \left|{\left(ud\right)}_R{u}_L\right|p⟩$	39.7%	1.0	31.7%	1.0	23.8%	1.4	9.4%	1.0	4.7%	1.6
$⟨\eta \left|{\left(ud\right)}_L{u}_L\right|p⟩$	2.8%	1.0	1.3%	1.2	2.5%	1.1	1.9%	1.8	1.5%	0.7

Table 4. Table 4: Table of renormalized W 0 subscript 𝑊 0 W_{0} in the physical kinematics at 2 GeV in MS ¯ ¯ MS \overline{\rm MS} NDR scheme. The fourth column presents relative error of systematic uncertainties; “ χ 𝜒 \chi ” is coming from chiral extrapolation given from three different fitting ranges as explained in the text, “ q 4 superscript 𝑞 4 q^{4} ” and “ a 2 superscript 𝑎 2 a^{2} ” column is uncertainty of higher order correction than 𝒪 ​ ( q 2 ) 𝒪 superscript 𝑞 2 \mathcal{O}(q^{2}) and lattice artifact at 𝒪 ​ ( a 2 ) 𝒪 superscript 𝑎 2 \mathcal{O}(a^{2}) respectively. “ m s subscript 𝑚 𝑠 m_{s} ” column is uncertainty for use of unphysical strange quark mass. Δ Z subscript Δ 𝑍 \Delta_{Z} and Δ a subscript Δ 𝑎 \Delta_{a} are error of renormalization factor and lattice scale estimate respectively.

Matrix element	$W_0$ GeV2	stat.[%]	Systematic error [%]
			total	$\chi$	$q^4$	$m{q}^2$	$a^2$	$m_s$	${\Delta }_a$	${\Delta }_Z$
$⟨{\pi }^0\left|{\left(ud\right)}_R{u}_L\right|p⟩$	-0.131(4)(13)	3.0	9.7	1.8	0.7	0.3	5.0	-	0.6	8.1
$⟨{\pi }^0\left|{\left(ud\right)}_L{u}_L\right|p⟩$	0.134(5)(16)	3.4	11.6	5.7	2.3	2.6
$⟨{\pi }^{+}\left|{\left(du\right)}_R{d}_L\right|p⟩$	-0.186(6)(18)	3.0	9.7	1.8	0.7	0.3
$⟨{\pi }^{+}\left|{\left(du\right)}_L{d}_L\right|p⟩$	0.189(6)(22)	3.4	11.6	5.7	2.3	2.6
$⟨K^0\left|{\left(us\right)}_R{u}_L\right|p⟩$	0.103(3)(11)	2.8	10.4	2.8	0.7	1.1	5.0	3.0	0.6	8.1
$⟨K^0\left|{\left(us\right)}_L{u}_L\right|p⟩$	0.057(2)(6)	3.5	10.7	3.1	1.0	2.1
$⟨K^{+}\left|{\left(us\right)}_R{d}_L\right|p⟩$	-0.049(2)(5)	3.7	10.9	3.5	1.6	2.0
$⟨K^{+}\left|{\left(us\right)}_L{d}_L\right|p⟩$	0.041(2)(5)	4.4	13.1	7.5	3.3	1.9
$⟨K^{+}\left|{\left(ud\right)}_R{s}_L\right|p⟩$	-0.134(4)(14)	3.2	10.3	1.6	1.3	1.3
$⟨K^{+}\left|{\left(ud\right)}_L{s}_L\right|p⟩$	0.139(4)(15)	3.0	10.9	3.9	1.4	1.5
$⟨K^{+}\left|{\left(ds\right)}_R{u}_L\right|p⟩$	-0.054(2)(6)	3.6	10.6	2.7	2.3	0.7
$⟨K^{+}\left|{\left(ds\right)}_L{u}_L\right|p⟩$	-0.098(3)(10)	2.8	10.3	2.1	0.8	1.6
$⟨\eta \left|{\left(ud\right)}_R{u}_L\right|p⟩$	0.006(2)(3)	30.0	42.1	39.7	9.4	4.7	5.0	-	0.6	8.1
$⟨\eta \left|{\left(ud\right)}_L{u}_L\right|p⟩$	0.113(3)(12)	3.1	10.2	2.8	1.9	1.5

Table 5. Table 5: Table of renormalized W μ subscript 𝑊 𝜇 W_{\mu} (Eq.  3 ), which is the form factor in the physical kinematics with final state of μ + superscript 𝜇 \mu^{+} .

Matrix element	$W_{\mu }$ GeV2
$⟨{\pi }^0\left|{\left(ud\right)}_R{u}_L\right|p⟩$	-0.118(3)(12)
$⟨{\pi }^0\left|{\left(ud\right)}_L{u}_L\right|p⟩$	0.119(4)(14)
$⟨{\pi }^{-}\left|{\left(du\right)}_R{u}_L\right|n⟩$	-0.167(4)(16)
$⟨{\pi }^{-}\left|{\left(du\right)}_L{u}_L\right|n⟩$	0.169(5)(20)
$⟨K^0\left|{\left(us\right)}_R{u}_L\right|p⟩$	0.099(2)(10)
$⟨K^0\left|{\left(us\right)}_L{u}_L\right|p⟩$	0.061(2)(7)
$⟨\eta \left|{\left(ud\right)}_R{u}_L\right|p⟩$	0.011(2)(3)
$⟨\eta \left|{\left(ud\right)}_L{u}_L\right|p⟩$	0.108(3)(11)

Table 6. Table 6: Comparison of matrix element calculation in lattice QCD from N f = 0 subscript 𝑁 𝑓 0 N_{f}=0 to N f = 3 subscript 𝑁 𝑓 3 N_{f}=3 with several groups. The errors of α 𝛼 \alpha , β 𝛽 \beta and W 0 subscript 𝑊 0 W_{0} are denoted as the total one, which is combined with statistical and systematic errors in the quadrature. We remark that ∗ denotes the result using renormalization constants with perturbative matching factors having an error as explained in the text. With the error corrected the values would increase as much as ∼ similar-to \sim 7 %.

Ref.	JLQCD	CP-PACS	RBC	QCDSF	RBC/	This
		& JLQCD			UKQCD	work
	(2000)	(2004)	(2007)	(2008)	(2008,2014)
	[29]	[31]	[30]	[46]	[32, 1]
Fermion	Wilson	Wilson	DW	Wilson	DW	DW
$N_f$	0	0	0 and 2	2	3	3
		(3.3)3	Quench	(1.68)3	(2.65)3	(2.65)3
Volume	(2.4)2		(1.6)3
(fm3)	$\times$4.1		Two-flavor
			(1.9)3
$a$ (fm)	0.09	0	Quench	0.07	0.11	0.11
			0.1
			Two-flavor
			0.12
	0.45–0.73	0.6–1.2	Quench	0.42–1.18	0.34–0.69	0.34–0.69
$m_{\pi }$			0.39–0.58
(GeV)			Two-flavor
			0.48–0.67
Renorm.	One-loop	One-loop	NPR	NPR	NPR	NPR
$\mu$	$1/a$, $\pi /a$	2 GeV	2 GeV	2 GeV	2 GeV	2 GeV
	$-0.015\left(1\right)$	$-0.0090{(}_{-21}^{+10})$	Quench	$-0.0091\left(4\right)$	$-{0.0119}^{\ast }\left(26\right)$	$-0.0144\left(15\right)$
$\alpha$			$-{0.0100}^{\ast }\left(19\right)$
(GeV3)			Two-flavor
			$-{0.0118}^{\ast }\left(21\right)$
	$0.014\left(1\right)$	$0.0096{(}_{-22}^{+11})$	Quench	$0.0090\left(4\right)$	${0.0128}^{\ast }\left(28\right)$	$0.0144\left(15\right)$
$\beta$			${0.0108}^{\ast }\left(21\right)$
(GeV3)			Two-flavor
			${0.0118}^{\ast }\left(21\right)$
$p\to {\pi }^0$
	$-0.134\left(16\right)$	-	Quench	-	$-{0.103}^{\ast }\left(41\right)$	$-0.131\left(13\right)$
$W_0^{LR}$			$-{0.060}^{\ast }\left(18\right)$
(GeV2)			Two-flavor
			-
	$0.128\left(17\right)$	-	Quench	-	${0.133}^{\ast }\left(40\right)$	$0.134\left(16\right)$
$W_0^{LL}$			${0.086}^{\ast }\left(22\right)$
(GeV2)			Two-flavor
			-

Equations

$$\Gamma(N\rightarrow P+\bar{l})=\frac{m_{N}}{32\pi}\Big{[}1-\Big{(}\frac{m_{P}}{m_{N}}\Big{)}^{2}\Big{]}^{2}\Big{|}\sum_{I}C^{I}W_{0}^{I}(N\rightarrow P)\Big{|}^{2}+O(m_{l}/m_{N}),$$

$$\langle P(p)|\mathcal{O}^{\Gamma\Gamma^{\prime}}(q)|N(k,s)\rangle=P_{\Gamma^{\prime}}\Big{[}W_{0}^{\Gamma\Gamma^{\prime}}(q^{2})-\frac{iq/}{m_{N}}W_{1}^{\Gamma\Gamma^{\prime}}(q^{2})\Big{]}u_{N}(k,s).$$

$$W_{\mu }\equiv W_0(-m_{\mu }^2)+\frac{m_{\mu }}{m_N}{W}_1(-m_{\mu }^2).$$

$$O^{\Gamma {\Gamma }^{\prime }}=U_{\Gamma {\Gamma }^{\prime }}{\epsilon }^{ijk}(q^{iT}C{P}_{\Gamma }{q}^j)P_{{\Gamma }^{\prime }}{q}^k,$$

$$⟨{\pi }^0|(ud{)}_{\Gamma }{u}_{{\Gamma }^{\prime }}|p⟩$$

$$⟨{\pi }^{+}|(ud{)}_{\Gamma }{d}_{{\Gamma }^{\prime }}|p⟩$$

$$⟨K^0|(us{)}_{\Gamma }{u}_{{\Gamma }^{\prime }}|p⟩$$

$$⟨K^{+}|(us{)}_{\Gamma }{d}_{{\Gamma }^{\prime }}|p⟩$$

$$⟨K^{+}|(ud{)}_{\Gamma }{s}_{{\Gamma }^{\prime }}|p⟩$$

$$⟨K^{+}|(ds{)}_{\Gamma }{u}_{{\Gamma }^{\prime }}|p⟩$$

$$⟨\eta |(ud{)}_{\Gamma }{u}_{{\Gamma }^{\prime }}|p⟩$$

$$⟨{\pi }^{+}|(ud{)}_{\Gamma }{d}_{{\Gamma }^{\prime }}|p⟩=\sqrt{2}⟨{\pi }^0|(ud{)}_{\Gamma }{u}_{{\Gamma }^{\prime }}|p⟩,$$

$$R_3^{\Gamma L}(t,t_1,t_0;\vec{p};P)$$

$$\underset{t_1-t,t-t_0\to \infty }{\lim }{R}_3^{\Gamma L}(t,t_1,t_0;\vec{p},P_4)$$

$$\underset{t_1-t,t-t_0\to \infty }{\lim }{R}_3^{\Gamma L}(t,t_1,t_0;\vec{p},i{P}_4{\gamma }_j)$$

$${\mathcal{N}}_5$$

$${\mathcal{N}}_{45}$$

$$U^{RL}=0.705(11)(56),U^{LL}=0.706(11)(56),$$

$$⟨0|(ud{)}_R{u}_L|p⟩=\alpha P_R{u}_p,⟨0|(ud{)}_L{u}_L|p⟩=\beta P_L{u}_p,$$

$$R_{\alpha }=\frac{C_{N\mathcal{O}}^R(t)}{C_N(t)}{Z}_N\stackrel{t\to \infty }{\to }\alpha ,R_{\beta }=\frac{C_{N\mathcal{O}}^L(t)}{C_N(t)}{Z}_N\stackrel{t\to \infty }{\to }\beta$$

$$f(\tilde{m})=c_0+c_1\tilde{m},$$

$${\tilde{m}}_{ud}^{\mathrm{phys}}=0.001382,$$

$$\alpha =-0.0144(3)(21)\mathrm{Ge}{\mathrm{V}}^3,\beta =0.0144(3)(21)\mathrm{Ge}{\mathrm{V}}^3,$$

$$F_{W_0}=A_0+A_1\tilde{m}+A_2{q}^2,$$

$$u_R{d}_R{d}_R{\Psi }_R\Phi /{\Lambda }^3+\text{hc},$$

$$⟨{\pi }^0|(ud{)}_R{u}_L|p⟩$$

$$⟨{\pi }^0|(ud{)}_L{u}_L|p⟩$$

$$⟨K^0|(us{)}_R{u}_L|p⟩$$

$$⟨K^0|(us{)}_L{u}_L|p⟩$$

$$⟨K^{+}|(us{)}_R{d}_L|p⟩$$

$$⟨K^{+}|(us{)}_L{d}_L|p⟩$$

$$⟨K^{+}|(ud{)}_R{s}_L|p⟩$$

$$⟨K^{+}|(ud{)}_L{s}_L|p⟩$$

$$⟨K^{+}|(ds{)}_R{u}_L|p⟩$$

$$⟨K^{+}|(ds{)}_L{u}_L|p⟩$$

$$⟨\eta |(ud{)}_R{u}_L|p⟩$$

$$⟨\eta |(ud{)}_L{u}_L|p⟩$$

$$⟨{\pi }^0|(ud{)}_R{u}_L|p{⟩}_{sp}=\frac{\alpha }{\sqrt{2}{f}_0}{P}_L{u}_N,⟨{\pi }^0|(ud{)}_L{u}_L|p{⟩}_{sp}=\frac{\beta }{\sqrt{2}{f}_0}{P}_L{u}_N,$$

$$⟨{\pi }^0|(ud{)}_{\Gamma }{u}_L|p{⟩}_{sp}=P_L{W}_{sp}{u}_N,$$

$$W_{\vec{p}=0}=\underset{t_1-t,t-t_0\to \infty }{\lim }{R}_3^{\Gamma }(t,t_1,t_0;0,P_4),$$

Full text

1 Introduction

Although proton decay has not been observed in the experiment yet, it is an important key observable for search of the new physics beyond the Standard Model (SM). The observed proton lifetime, i.e. $\tau_{p}>8.2\times 10^{33}$ year for $p\rightarrow\pi^{0}e^{+}$ [2] (recently $\tau_{p}>1.4\times 10^{34}$ year has been reported in [3]) or $\tau_{p}>5.9\times 10^{33}$ year for $p\rightarrow K^{+}\bar{\nu}$ [4], impose tight constraints to the parameter space of the Grand Unified Theories (GUTs) and supersymmetric one (SUSY-GUTs). Currently such an experimental bound might exclude minimal SU(5) GUTs, besides SUSY-GUTs has been attractive models for the solution of hierarchy problem and the coupling unification of the SM in the GUT scale ($\sim 10^{16}$ GeV). SUSY-GUTs favor $p\rightarrow K^{+}\bar{\nu}$ decay channel within the detectable region of proton decay in the future experiment (e.g. Hyper-Kamiokande [5]).

The main mode of the proton decay through GUTs is those where a proton decays into a pseudoscalar meson and an anti-lepton. The operator product expansion (OPE) leads to the decay amplitude of such processes written in terms of the Wilson coefficients which contain all the details of the high energy part of a GUT, and the low energy QCD matrix elements of the proton and pseudoscalar states with three-quark operators. Each QCD matrix element is further decomposed into two form factors, named relevant and irrelevant form factor. Denoting the relevant form factor $W_{0}$, the partial decay width reads

[TABLE]

with $m_{N}$, $m_{P}$ and $m_{l}$ being the mass of the nucleon, pseudoscalar meson and anti-lepton, and $C^{I}$ being the Wilson coefficient of the operator of type $I$ (distinguishing flavor and chiral structure), which entering also in $W_{0}^{I}$. Parameters in a given GUT model are encoded in the Wilson coefficients $C^{I}$. The knowledge of the left hand side (experiment) and that of $W_{0}$ reported in this work will be transcribed into the knowledge of the GUT parameters. Namely the proton lifetime bound gives rise to restricting the GUT parameters [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24].

The relevant form factors are evaluated in $\overline{\rm MS}$ scheme in the naive dimensional regularization at a typical hadronic scale $\mu=2$ GeV. The matching Wilson coefficients need to be calculated in the same way.

Lattice computation of the proton decay matrix elements has rather a long history. It is started with the calculation of low energy constants (LEC) $\alpha$ and $\beta$ in quenched approximation $N_{f}=0$, where one needs to use a baryon chiral perturbation theory [25] to obtain $W_{0}$. The uncertainties of these initial computations [26, 27, 28] have been successively reduced with systematic improvements by employing the direct method [29, 30, 1], continuum limit of LEC’s in $N_{f}=0$ [31], non-perturbative renormalization with chirally invariant lattice formulation [30, 32, 1], LEC’s computed in $N_{f}=2+1$ [32], and finally $W_{0}$ computation in $N_{f}=2+1$ (see Table 6).

In the previous report [1] $W_{0}$ is calculated using the direct method with proper dynamical fermion computation using the $N_{f}=2+1$ domain-wall fermion formulation. This paper reported the result of $W_{0}$ with all the relevant systematic uncertainties removed or properly estimated. The precision, though, was not satisfactory because $W_{0}$ for the pion and kaon final state matrix elements have 20–40 % errors. Noticing the fractional error gets doubled in the partial decay width as it enters quadratically in Eq. (1), it is necessary to have more precise results for $W_{0}$ in order to make them more useful. Reduction of the statistical error into sub-dominance is essential task for this study, since it has been dominated as a half and more in total error.

A recent development of the algorithm to speed-up the measurement of matrix element in lattice QCD, called as all-mode-averaging (AMA) [33, 34, 35], enables us to further improve the statistical precision of the proton decay matrix element for the pion and kaon final states. This paper shows the update of lattice calculation of proton decay matrix element, in both “direct” and “indirect” measurements, for all decay modes on the same gauge ensembles as used in [1]. As a consequence of increasing statistical accuracy, more reliable estimate of the systematic error can be realized.

In addition, previously we have not taken into account the muon mass effect for the case of muon final state because the effect is sub-dominant compared to other uncertainty. However, with increased statistical accuracy, the effect of non-zero muon mass (106 MeV) is visible. We provide the form factors for the muon final state separately from those with positron or neutrino final states.

Here we also attempt to use our lattice computation as an input of proton decay matrix element for the model of dark matter [36, 37, 38]. Although the kinematics of our setup is not optimal for those needed for the dark matter scattering there, we provide information for them as useful bi-products.

This paper is organized as follows; after showing the notation (Section 2) and simulation parameters (Section 3), we show the updated result of lattice evaluation of the low energy constants $\alpha$ and $\beta$ in BChPT for the indirect method in Section 4, and relevant form factor $W_{0}$ of proton decay in Section 5. In Section 5 we also make an assessment of the unestimated systematic error in the indirect method and make alert to the use of them in the estimate of the proton lifetime. A description of how our results can be used in a dark matter model [36, 37, 38] is given in Section 6. A test of the soft-pion theorem of our lattice results is shown in Appendix B. Throughout the paper dimension-full quantities are expressed in the lattice unit and the lattice spacing “$a$” is suppressed in equations.

2 Proton decay matrix element

Lattice calculation is concentrating on the QCD matrix element of $N\rightarrow P$ transition, in which $N$ denotes the nucleon (proton or neutron) , and $P$ is one of pseudoscalar from $\pi^{0}$, $\pi^{\pm}$, $K^{0}$, $K^{+}$ or $\eta$ mesons. At the hadronic energy scale, only lowest dimensional operators with baryon number violation are relevant. They are the dimension-six four-fermi (three quarks and one lepton) operators [39, 40, 41]. Since the on-shell lepton can be omitted from QCD matrix element, the transition form factor from a nucleon $N(k)$ state (source field) to meson $P(p)$ state (sink field) with a momentum transfer $q=k-p$ is represented as

[TABLE]

On the physical kinematics $-q^{2}=m_{l}^{2}$, the contribution of the $W_{1}$ term is relatively small compared to the $W_{0}$ term, because for the suppression prefactor $m_{l}/m_{N}$. In $m_{l}=e^{+}$ or $\bar{\nu}$, $W_{0}$ is only relevant to proton decay matrix element since the second term is $m_{e}/m_{N}\sim\mathcal{O}(10^{-3})$, while in $m_{l}=m_{\mu}$ case, because of $m_{\mu}/m_{N}\sim\mathcal{O}(10^{-1})$, the contribution of $W_{1}$ term to the matrix element is not negligible for our target precision (below 10% precision), namely we define

[TABLE]

We estimate $W_{0}$ and $W_{1}$ and provide $W_{0}(q^{2}=0)$ for the positron and neutrino final states and $W_{\mu}$ for the anti-muon final states. Baryon number violating three-quark operator $O^{\Gamma\Gamma^{\prime}}$ reads

[TABLE]

with chiral projection $P_{\Gamma}=(1\pm\gamma_{5})/2$, where “+” is for $\Gamma=R$ and “$-$“ is for $\Gamma=L$. $q^{i}$ is quark flavor of up, down and strange with color index $i$. $U_{\Gamma\Gamma^{\prime}}$ denotes the renormalization factor which has been already computed by non-perturbative method [32] 222Recently we found an error in our one-loop perturbative formula which is used for matching between $\overline{\rm MS}$ with naive dimensional regularization and RI-SMOM renormalization scheme, Eq. (C.8) in [30] and Eq. (46) in [32] (we thank Michael Buchoff and Michael Wagman for pointing out that mistake). After correction, the impact of all matrix element calculations is about 6–7% increase for $\alpha$, $\beta$ and $W_{0}$ [30, 32, 42]. In this paper we use the corrected value as presented in Eq. (18).. Using the symmetry of parity transformation between different chirality combinations $RL\Leftrightarrow LR$ or $LL\Leftrightarrow RR$ [30] enables us to reduce four chirality combinations to two combinations, $\Gamma\Gamma^{\prime}=LL$ and $RL$. Applying the exchange symmetry of $u$ and $d$, we have the equivalence of matrix elements between proton and neutron,

[TABLE]

and furthermore, for $p\rightarrow\pi$ channel, there is a relation in the SU(2) isospin limit, which is good for our target precision,

[TABLE]

and therefore the total number of matrix element ends up to twelve. In this paper we show twelve principal matrix elements of $\langle P|\mathcal{O}^{\Gamma L}|p\rangle$, for $\Gamma=R$ and $L$ in lattice QCD.

Our target matrix element can be extracted from the computation of ratio given from three-point function and two-point function. We use the same combination used in Eq. (21) of [42],

[TABLE]

with nucleon two-point function $C_{N}(t,t_{0})$ without momentum, and pseudoscalar two-point function $C_{P}(t,t_{0};p)$ with spatial momentum $\vec{p}$. Here we also use the two projection matrices, $P=P_{4}\equiv(1+\gamma_{4})/2$ and $iP_{4}\gamma_{j}$. Three-point function $C_{\mathcal{O}^{\Gamma L}}(t_{1},t,t_{0};p)$ depends on $t-t_{0}$, the time-slice of operator, and, $t_{s}=t_{1}-t_{0}$, source-sink separation, and also injected momentum $\vec{p}$ in the operator. The factors $\sqrt{Z_{P}}$ and $\sqrt{Z_{N}}$ are overlap factors of the pseudoscalar and nucleon states to their interpolating operators. Asymptotic form of this ratio taking the trace with two projection matrices $P_{4}$ and $iP_{4}\gamma_{j}$ can be expressed as the combination of $W_{0}$ and $W_{1}$

[TABLE]

and solving the linear algebra we derive $W_{0}^{\Gamma L}$ and $W_{1}^{\Gamma L}$ simultaneously.

Calculating the three-point function $C_{\mathcal{O}^{\Gamma L}}(t_{1},t,t_{0};p)$ in Eq. (15) involves several steps. First we compute the forward quark propagator with the nucleon source located at $t=t_{0}$ with a smeared source. Then using the propagator at the meson sink position $t=t_{1}$ the sequential source computation is applied with an injection of momentum $\vec{p}$. Then the obtained backward propagator is contracted at the operator position $t=t$ with two forward propagators from the nucleon source. This process needs $1+N_{p}\times 2$ solver computations for each gauge configuration, with $N_{p}$ being the number of different meson momenta. “1” is for the forward propagator and the distinction of valence mass for the $ud$ and $s$ quarks makes the factor “2”. For the good constraint on the fitting parameters, we need to have good lever arm for $m$ (different ensembles) and variation of $\vec{p}$, which tends to sum up a large computational cost. The all-mode-averaging (AMA) technique [33, 34, 35, 43] is useful to reduce the computational cost of quark propagator by using this method. It enables us to carry out the high statistical measurement even using several momenta.

3 Lattice setup

We use the same lattice gauge ensembles as used in [42, 44, 45], which are generated with $N_{f}=2+1$ the domain-wall fermion (DWF) and Iwasaki gauge action at $\beta=2.13$, corresponding to $a^{-1}=1.7848(6)$ GeV [45], in 24${}^{3}\times$64 lattice size ($\simeq 2.65$ fm3). The four different quark masses, $m=0.005$, 0.01, 0.02 and 0.03 are used in unitary point and the corresponding pion, kaon and nucleon masses are given in Table 1. In the measurement of two-point function of pseudoscalar meson and nucleon, we use the gauge invariant Gaussian smeared source and sink function with interpolation operator on APE smeared link variable, whose parameter is same as [42]. For “indirect” method, two-point function including baryon number violating operator Eq.(4) is computed using two nucleon source operators

[TABLE]

and averaged them. On the other hand, in “direct” method for the two-point function in the ratio in Eq. (13), we use only $\mathcal{N}_{5}$, as proton interpolation operator.

The renormalization factor to make the operators to ones in $\overline{\rm MS}$ naive dimensional regularization (NDR) scheme at $\mu=2$ GeV is calculated by RI/MOM non-perturbative renormalization method combined with the RI/MOM $\to$ $\overline{\rm MS}$ matching factor calculated to the next-to-leading order in perturbation theory. The renormalization factors of $\mathcal{O}^{\Gamma}$ with $\overline{\rm MS}$ NDR at $\mu=2$ GeV are given as

[TABLE]

where the first error is statistical, and second is systematic. The systematic error is dominated by the truncation error in perturbative matching, which is done to next-to-leading order. The estimate is from the size of $\alpha_{s}^{2}(\mu=2\mbox{GeV})$ obtained by the RGE running starting from $\alpha_{s}(M_{Z})=0.1176(2)$ (see [30]).

In the computation of three-point function, we use $t_{s}=18$ ($\simeq 1.98$ fm) for source-sink separation which is shorter than $t_{s}=22$ ($\simeq 2.43$ fm) used in [1]. Although the usage of short source-sink separation will make a suppression of the statistical noise, we need to make sure the excited state contamination is negligible. As the contamination would be more serious for smaller quark mass, we test such a contamination effect by comparing the ratio $R_{3}^{\Gamma L}$ for $t_{s}=18$ and $22$ at the lightest quark mass in Section 4.

We use three non-zero spatial momenta for the mesons: $\vec{p}=(1,0,0)$, (1,1,0) and (1,1,1), where the last one is a new addition from the previous study [1]. This will provide a good lever arm for the $p^{2}$ direction as well as the good reach for the momentum range in the different kinematics (see Section 6).

The AMA technique is applied to the measurement of the three-point and two-point functions 333For the two-point function $C_{P}(t_{1},t;\vec{p})$ of pion (or eta) in the denominator of Eq. (13) we have used Kuramashi-wall source as in [42] for heavier “light” quark mass, $m=0.02$ and $m=0.03$, since there is less gain for cost per precision of signal, and thus AMA was not applied in this case.. The low-mode deflation is used in solving the even-odd decomposed Dirac kernel with conjugate gradient method. Corresponding low-mode is computed by Lanczos algorithm with Chebyshev polynomial acceleration as performed in [35]. The number of low-modes $N_{\rm eig}$ we computed for each quark mass are given in Table 1. Approximation used in AMA is also constructed by the sloppy solver using relaxed stopping condition (0.003 for the squared norm, which is compared with the $10^{-8}$ for the “exact solve” done once every configuration). $N_{g}$ shown in Table 1 presents the number of such an approximation we use in AMA. Note that in the strange quark mass we use sloppy solver without deflation to avoid the additional computation of low-mode. Even without low-mode of strange quark, AMA is also working well. Actually we check that correlation between exact and approximation is smaller than $1/N_{g}$.

4 Improved result of low energy constants

First we update the “indirect” measurement from computation of low energy constants (LECs) for baryon number violating interaction in chiral Lagrangian [25] following the method in [29, 31, 30, 32, 46]. In the “indirect” measurement, once corresponding LECs are obtained by lattice QCD, through baryon chiral perturbation theory (BChPT) together with nucleon mass, couplings to axial current (axial charge), pion decay constant and its mass, the proton decay amplitude can be evaluated. Each matrix element is proportional to LECs depending on chirality; $\alpha$ (for $RL$) and $\beta$ (for $LL$) [25, 29]. Those are defined through the nucleon to vacuum matrix elements. Writing the quark flavor explicitly

[TABLE]

with proton spinor field $u_{p}$. The above matrix elements can be extracted from the ratio of two-point function at large time-slice separation,

[TABLE]

with the nucleon decay operator ${\mathcal{O}}$ and the nucleon interpolating operator having the same flavor content. The nucleon overlap factor $Z_{N}$ is also calculated from the nucleon two-point function. From the practical point of view, this method is much cheaper than the “direct” method, since $\alpha$ and $\beta$ are obtained with single computation of quark propagator at each quark mass. Whereas, the direct method needs at least additional two propagators for each momentum values in the computation of three-point function.

Figure 1 plots our result of $R_{\alpha}(t)$ and $R_{\beta}(t)$ obtained after averaging those with two different nucleon interpolating operators $\mathcal{N}_{4}$ and $\mathcal{N}_{45}$. Fitting to plateau is done to the range $t\in[8,18]$ for all cases as we have shown the straight bar in Figure 1, where the statistical error is included.

Figure 2 shows the quark mass dependence of bare value of $\alpha$ and $\beta$. We observe that lattice data behaves as a linear function in our quark mass region, and the chiral extrapolation to physical quark mass is carried out with linear function of quark mass,

[TABLE]

with $\tilde{m}=m+m_{\rm res}$, where the residual mass has been estimated as $m_{\rm res}=0.003152(43)$ [42]. The (bare) physical quark mass has been obtained as

[TABLE]

from renormalized one $m^{\rm phys}_{ud}=-0.001770(79)$ [45].

To estimate the uncertainty due to using the linear extrapolation, we use three different fitting ranges: (i) $m\in[0.005,0.03]$ (3.4,2.7), (ii) $m\in[0.005,0.02]$ (2.5,1.1) and (iii) $m\in[0.01,0.03]$ (2.5,2.5), where the number in the brackets show the resultant $\chi^{2}$/dof for $\alpha$ and $\beta$ respectively. The systematic error due to the assumption of linear behavior is evaluated from the maximum difference of central value between (i) and (ii), and (i) and (iii). The results are tabulated in Table 2. Compared to the previous work [32], the statistical error has been improved to 2% from 10%, and systematic error of the chiral extrapolation is improved to 3% from 20%. Although the error estimation procedure is same as the previous work, because of the highly statistical precision we can use, the systematic error is properly estimated.

We estimate the systematic error of lattice artifact as 5%, which is evaluated from comparison with different lattice spacing for hadron spectrum (see [42]). The uncertainty in renormalization factor, which is dominated by the truncation error of the perturbative matching and running beyond the next-to-leading order. This error turns out to be the most dominant error.

The final value at $\mu=2$ GeV in $\overline{\rm MS}$ NDR scheme 18 extrapolated to physical quark Eq. (22) is

[TABLE]

where the first error is statistical and the second is systematic obtained by the quadrature. The total error is around 15%, which is improved from 22% [32]. The superficial relation $\alpha+\beta=0$ is observed. The relation should hold for the non-relativistic limit and the approximate relation is known to hold at least numerically in the quenched case [30]. Here we have confirmed that it holds in the $N_{f}=2+1$ case with an improved precision. Using these low energy constants, the relevant form factor can be computed via BChPT formula (see Appendix A).

5 Improved result of relevant form factor

In this section, we show our improved result of the form factors from “direct” measurement in which we compute the three-point function of $N\rightarrow P$ including baryon number violating operator. Compared to our previous study [1], the results are improved by the use of the AMA technique. We also add one larger meson momentum point, $n_{p}=(1,1,1)$ to the two non-zero momentum we had, $n_{p}=(1,0,0)$ and $n_{p}=(1,1,0)$. By that we now are able to estimate the systematic error from $\mathcal{O}(q^{4})$ term.

In Figure 3, we plot the effective mass of nucleon, pion and kaon with momentum we use in the construction of ratio Eq. (13). One clearly sees the plateau starting from $t=6$ in those hadrons, so from here we regard that the ground state is dominant from $t=6$.

Figure 4, in which we plot the time-slice dependence of the matrix element extracted from the ratio Eq. (13) for the $p\to\pi^{0}$ mode at our lightest point $m=0.005$, shows the comparison with two different source-sink separations, $t_{s}=18$ and 22 corresponding to $t_{s}=1.98$ and 2.43 fm respectively. The time separation $t_{s}=18$ is new in this study with four-time slice shorter than original $t_{s}=22$ [1]. We observe plateau for shorter separation at $t\in[12,16]$, where the denominators are also dominated with the ground state (see Fig. 3 and note that the source is located at $t=5$ here). The plateaus from two separations are consistent and the shorter separation yields significantly better statistical accuracy. Let us note that the clear plateau is observed even in largest momentum case $n_{p}=(1,1,1)$. Hereafter we only use the result in $t_{s}=18$, and further test the effect of excited state contamination by changing the fitting range below.

Figure 5 shows the result of plateau fit for $W_{0}$ using variation of fitting ranges as $t\in[11,17]$ and $t\in[12,16]$ to study the effect of excited state contamination into the signal. We observe those values are consistent within 1 sigma error in each $q^{2}$, while the central value has slight tension, especially for the lowest momentum in $\Gamma=L$. In order for estimate of systematic uncertainties including the effect of excited state contamination, we compare the results using those fitting ranges. We will back to such a discussion later.

5.1 Global fitting

To perform the extrapolation to kinematic point and physical pion mass simultaneously, we globally fit all lattice data with the linear ansatz for quark mass and $q^{2}$ dependence as

[TABLE]

where $\tilde{m}$ is the same definition as in eq.(21). Figures 6 and 7 plot renormalized $W_{0}(q^{2})$ for every decay channel in each quark mass. We observe that lattice data in each quark mass, which denotes the same symbols in Figures 6 and 7, is behaved as linear $q^{2}$ dependence. For mass dependence, we also observe the monotonic decreasing or increasing when $m$ is increasing. Even using the linear ansatz $\chi^{2}$/dof is reasonably small (note that we use uncorrelated fits) as presented in the first “$\chi^{2}_{\rm dof}$” column of Table 3.

We next study the uncertainties in the fitting related with the mass dependence, following the method used in Ref. [1]. The estimated errors are attributed to the higher order correction than $\mathcal{O}(m)$ and a part (at least) of the finite volume effect. Table 3 presents the errors estimated with the discrepancy from central value, which is obtained by full range, $m\in[0.005,0.03]$, and two fitting ranges; $m\in[0.005,0.02]$ for “light” region and $m\in[0.01,0.03]$ for “heavy” region. The error with“light” region can be an estimate of $\mathcal{O}(m^{2})$ correction since exclusion of heavy mass makes less $\mathcal{O}(m^{2})$ correction. On the other hand, the error with “heavy” region can (at least partly) be due to the finite volume effect, since the lightest point suffers most from the effect with the fixed volume. In each range, $\chi^{2}$/dof is not significantly large. The values presented in the table are taken as the maximum error compared with the result obtained in two $t$ fitting ranges $t\in[12,16]$ and $t\in[11,17]$. The “total” error of the chiral extrapolation in the table is calculated by adding two errors, “light” and “heavy”, in quadrature.

In similar manner as “light” error, the $\mathcal{O}(q^{4})$ error is estimated from the difference of the results obtained with the full range of $q^{2}$ with all the non-zero meson momentum and the shorter range where largest $|q^{2}|$ ($n_{p}=(1,1,1)$) is neglected. The result is shown in the column labeled as $\mathcal{O}(q^{4})$ in Table 3. This error turns out to be smaller than that of the “chiral” extrapolation.

As shown in Figure 6 and 7, the $q^{2}$ dependence obtained by extrapolation of data in“direct” method does not largely differ from BChPT including $\alpha$ and $\beta$ obtained in section 4, especially for that $p\rightarrow\pi$ channel has a tendency to be close to each other when increasing $q^{2}>0$. There is a discrepancy up to about a factor of 2 around the kinematics point. Such a comparison will be discussed later.

5.2 Sequential fitting

In the global fitting we estimated a part of the systematic errors due to omitting the higher order terms in the expansion of the light quark mass $\tilde{m}$ and squared momentum transfer $q^{2}$. Those estimated are of $\mathcal{O}(\tilde{m}^{2})$ and $\mathcal{O}(q^{4})$. Remaining error is of $\mathcal{O}(\tilde{m}q^{2})$. For the estimate we use the same method as in Ref. [1]. The procedure is that first the $q^{2}\to 0$ extrapolation is carried out for each fixed quark mass with linear function, then chiral extrapolation is performed (see Figure 8 and 9). By doing that we are taking into account the $q^{2}$ dependence in prefactor of the linear quark mass term, $A_{1}$ in Eq. (24). If the result is different, it is attributed as the $m$ effect in $A_{1}$, thus is of $\mathcal{O}(\tilde{m}q^{2})$. The last column of Table 3 shows the $\chi^{2}$ per degree of freedom for the final $\tilde{m}$ linear fit. The second last column represents the systematic error estimated in this analysis. It turns out to be sub-dominant in the fitting errors.

5.3 The final results

Table 4 presents the summary of the nucleon decay form factor for each operator and final state with the statistical and systematic errors. The statistical error is significantly reduced to 1/4–1/6 from our previous study [1] and now is sub-dominant. The systematic errors for the extrapolation discussed above are combined and shown in the “$(\tilde{m}q^{2})$-fit” column. Since we use a single lattice cutoff in this study, the lattice artifact, which is $\mathcal{O}(a^{2})$ correction, is estimated from the scaling study of hadron spectrum as done in [42]. The mass of valence strange quark which participate in the matrix elements of kaon final state is set equal to its sea-quark mass $m=0.04$. There is a mismatch to the physical strange mass. The associated systematic error is estimated using a subset of the $m=0.005$ ensemble by setting $m=0.343$ 444the value comes from the physical strange quark mass used in the previous study. The latest estimate of physical strange mass [45] turns out to be 0.03224 which is not far enough to change the systematic error estimate.. The difference of central value is by 3 % at most. We conservatively take 3 % as the systematic error of the form factors for the process with the kaon final state due to the use of the mismatched strange sea and valence quark masses. On the other hand, the mismatch effect of the strange sea quark is expected to be much less than that of the valence quark, thus, it is negligible in pion and eta final state. The largest uncertainty comes from that of the renormalization factor, which is dominated by the systematic error due to the truncation of the perturbative matching (Eq. (18). The total error summing up all in quadrature amounts to 10–15 % for the form factors with pion and kaon final state.

Additionally, Table 5 presents the matrix element with muon final state, $m_{l}=m_{\mu}$. $W_{\mu}$ in Eq.(3) is made from two form factors, $W_{0}$ and $W_{1}$. As one sees in Figure 10 and 11, the magnitude of $W_{1}$ in each matrix element is similar to $W_{0}$, and hence $W_{1}$ term multiplied with factor $m_{\mu}/m_{N}\sim 0.1$ in $W_{\mu}$ affects around 10% effect to matrix element in the kinematics with muon final state.

Note that for matrix element with eta final state we are ignoring the disconnected diagram, which means there remains additional uncertainty. However, the contribution of disconnected diagram expects to be small from OZI suppression. Detailed study in eta sector including disconnected diagram is beyond the scope of this paper.

All the final results of the relevant form factors of proton decay $W_{0}$ and $W_{\mu}$ with the “direct” method are summarized in Figure 12. The results are also compared with those with the “indirect” method through BChPT using lattice LECs (denoted as $W^{\alpha,\beta}_{0}$ and $W^{\alpha,\beta}_{\mu}$). The “indirect” method always overestimates the form factor. The amount is 25% or more except for two cases ($\langle K^{+}|(us)_{R/L}d_{L}|p\rangle$). In contrast to previous study [1], each error becomes lot smaller, and now we clearly see the discrepancy between $W_{0,\mu}$ and $W^{\alpha,\beta}_{0,\mu}$ for most cases.

The fact that the indirect method which uses BChPT works poorly is understandable as the physical kinematical point for the outgoing pion is far from the soft pion limit, where the ChPT description becomes arbitrary precise. We tested the soft pion theorem for the form factors of the pion final state, which is found in the appendix B. There the results from the indirect and direct method appears to be consistent with each other in the soft pion limit.

6 Application to the kinematics of dark matter model

In this section, we present a demonstration of the interesting applications to the model using other kinematics, in which energetic pion is emitted from proton and dark matter (DM) appears instead of lepton. According to [36, 37, 38], the so-called “induced nucleon decay (IND)” scenario, the proton should decay to DM particles, $(\Psi,\Phi)$, having the anti-baryon number with mass $m_{\Phi,\Psi}\sim 2$–3 GeV. This model, motivated by hypothesis of asymmetric DM model [47], assumes the net baryon number in the Universe is symmetric, in which the SM particle sector has a baryon number asymmetry while the particle $X_{1}$ in hidden sector has opposite asymmetry, and DM $(\Psi,\Phi)$ has been generated from $X_{1}$ decay in the early Universe. Under a consistency with Sakharov condition, such decay should have baryon number violation and CP violation in non-thermal circumstance. In IND model, nucleon and pseudoscalar are interacting with DM through $X_{1}$, and thus scattering process $\Psi N\rightarrow\Phi^{\dagger}P$ and $\Phi N\rightarrow\bar{\Psi}P$ occur. The interesting feature of this model is that the QCD matrix element is same as that of the standard nucleon decay, since the operator related to DM scattering is composed of effective three-quark interaction,

[TABLE]

and only difference is its kinematics of which $q^{2}$ is different from on-shell lepton. In principle lattice calculation is accessible to the matrix element at $q^{2}$ values relevant to this model, and so that we can also provide more accurate value for the prediction of this model.

The DM mass $m_{\Phi,\Psi}\sim$2–3 GeV is predicted from cosmological observation and DM stability, which is much heavier than lepton mass, so that under momentum conservation pion has finite momentum, which is a shifted region to $-q^{2}<0$, (right direction from zero in Figure 6, 7, 10 and 11). Recalling the formula of transition form factor in Eq. (2), relevant form factor is both $W_{0}$ and $W_{1}$, since DM mass is heavy, $q^{2}\sim 4m^{2}_{\Phi,\Psi}>m_{p}^{2}$. Typical meson momentum in IND model is $|\vec{p}|=1$ GeV, in which the kinematics of IND model is $q^{2}\simeq 1$ GeV2. Figure 13 plots $W_{0,1}$ and $W_{0,1}^{\alpha,\beta}$ extrapolated to $q^{2}=1$ GeV2 using lattice results. Focusing on the pion channel, one sees that “direct” lattice calculation provides 25–50% value of $W_{0,1}^{\alpha,\beta}$ used for an estimate of proton lifetime in IND model [36, 37, 38]. Concerning the convergence issue of BChPT at tree-level applying to energetic meson arises in this kinematics, our lattice result indicates such a difference from an evaluation based on tree-level BChPT may not be negligible. One sees that possible effect to proton decay amplitude when using $W_{0,1}$ in our results may be factor 4 and more suppression to the results obtained with BChPT. This potentially large systematic error needs to be considered when one use the BChPT for this purpose.

7 Summary and discussion

In this paper we present improved computation of proton decay matrix element using the all-mode-averaging technique on the $N_{f}=2+1$ domain-wall fermion configurations. Compared to previous work [1] (also see Table 6), the statistical error has been significantly reduced for both low-energy constant in baryon chiral perturbation theory (BChPT) and matrix element extracted from three-point function. Our analysis using the precise lattice data with three variation of momentum, by which we add one more higher $q^{2}$, can evaluate higher order correction than $\mathcal{O}(q^{2})$. The systematic uncertainty for the chiral extrapolation due to using unphysical pion around $m_{\pi}=0.33$ GeV, is still large rather than $\mathcal{O}(q^{4})$, $\mathcal{O}(mq^{2})$ correction, while its magnitude strongly depends on the chirality of baryon number violating operator. This uncertainty can be reduced by using larger volume than 3 fm3 in physical pion mass generated by RBC-UKQCD collaboration [45] in future work. Currently the dominated error is coming from the uncertainty of renormalization factor and lattice artifact correction, and those may be reduced by the further effort of renormalization scheme and comparison with finer lattice [45]. Final result of $W_{0,\mu}$ is presented in Table 4 and 5, in which the total error in pion channel for both $e,\nu$ and $\mu$ final state is 10–14%, and of kaon sector is also of similar precision. Compared to $W_{0,\mu}^{\alpha,\beta}$ via “indirect” method using the improved lattice calculation of LECs, $W_{0,\mu}$ from “direct” method is 1.3–1.4 time small for proton decay amplitude. This means, if our result of $W_{0,\mu}$ is incorporated into GUT model prediction instead of $W_{0,\mu}^{\alpha,\beta}$, the proton lifetime prediction may become about 2 times larger. Finally we note that our calculation is also applicable to the other kinematics corresponding to a dark matter model, and pointing out that there will be higher order correction than NLO BChPT. Our lattice calculation of $W_{0,\mu}$ can provide more reliable value for such a model.

Acknowledgments

We thank members of RIKEN-BNL-Columbia (RBC) and UKQCD collaboration for sharing USQCD resources for part of our calculation. ES thanks Hooman Davoudiasl for a useful discussion. Numerical calculations were performed using the RICC at RIKEN and the Ds cluster at FNAL. This work was supported by the JSPS KAKENHI Grant, Nos. JP22540301 (TI), JP22224003 (YA), JP16K05320 (YA), MEXT KAKENHI Grant, Nos. JP23105714 (ES), JP23105715 (TI), and U.S. DOE grants DE-SC0012704 (TI and AS). We are grateful to BNL, the RIKEN BNL Research Center, RIKEN Advanced Center for Computing and Communication (ACCC), and USQCD for providing resources necessary for completion of this work. ES also thanks the INT and organizers of Program INT-15-3 “Intersections of BSM Phenomenology and QCD for New Physics Searches”, September 14 - October 23, 2015, and Ryuichiro Kitano for his support from MEXT Grant-in-Aid for Scientific Research on Innovative Areas (No. JP25105011).

Appendix A Leading formula of proton decay matrix element in BChPT

According to BChPT [25, 29], the relevant matrix element, $W^{\alpha,\beta}_{0}$, can be represented as

[TABLE]

with low-energy parameters $D=0.80$, $F=0.47$. In this paper, we use $f=0.131$ GeV, $m_{N}=0.94$ GeV and $m_{B}=1.15$ GeV [32].

Appendix B Test of soft-pion theorem

In this section, we present the analysis of matrix element in the soft-pion limit. In this limit, each matrix element is described in term of the leading order of BChPT. In order to test the lattice calculation can make a consistent value with BChPT in the soft-pion limit, we calculate matrix element with two ways; one is BChPT using LECs $\alpha$ and $\beta$ obtained by “indirect” lattice calculation and the second is matrix element obtained by “direct” lattice calculation. Using LECs the matrix element is

[TABLE]

with subscription $sp$ denoting the soft-pion limit, which corresponds to $p_{\mu}\rightarrow 0$ and chiral limit. On the other hand, the left-hand-side of the above equation is also represented as,

[TABLE]

in which $W_{sp}$ is obtained from the form factor at $\vec{p}=(0,0,0)$ for Eq. (2) in the chiral limit. We define such a form factor as

[TABLE]

and taking the extrapolation into zero quark mass with the linear ansatz,

[TABLE]

Linear ansatz is under the assumption of negligibly small $\sqrt{\tilde{m}}\sim m_{\pi}$ term even in $\tilde{m}\simeq 0$.

Figure 14 shows the lattice result of $W_{\vec{p}=0}$ and $W_{sp}$. We also show the lines of chiral extrapolation and the extrapolated values with linear ansatz in the chiral limit. Here we estimate the systematic uncertainties due to chiral extrapolation by comparing the “light” and “heavy” region as well as in Table 3. We also add the uncertainties of renormalization factor and lattice artifact same as in Table 4. One sees that the lattice data is close to linear function and the extrapolated value is consistent with BChPT within 1 sigma error. We notice that there is no visible curvature as the square-root of quark mass. It indicates that a coefficient of square-root of quark mass may not be significantly large.