Lattice QCD determination of neutron-antineutron matrix elements with physical quark masses

TL;DR

This paper presents the first lattice QCD calculations of neutron-antineutron transition matrix elements using physical quark masses, providing crucial data for new physics constraints from experimental searches.

Contribution

It introduces a detailed lattice QCD methodology for calculating neutron-antineutron matrix elements with physical quark masses, including systematic uncertainty analysis and non-perturbative renormalization.

Findings

Matrix elements are significantly larger than previous MIT bag model estimates.

Calculations performed with chiral-symmetric discretization and physical quark masses.

Systematic uncertainties from discretization and finite-volume effects need further study.

Abstract

Matrix elements of six-quark operators are needed to extract new physics constraints from experimental searches for neutron-antineutron oscillations. This work presents in detail the first lattice quantum chromodynamics calculations of the necessary neutron-antineutron transition matrix elements including calculation methods and discussions of systematic uncertainties. Implications of isospin and chiral symmetry on the matrix elements, power counting in the isospin limit, and renormalization of a chiral basis of six-quark operators are discussed. Calculations are performed with a chiral-symmetric discretization of the quark action and physical light quark masses in order to avoid the need for chiral extrapolation. Non-perturbative renormalization is performed, including a study of lattice cutoff effects. Excited-state effects are studied using two nucleon operators and multiple values of source-sink separation. Results for the dominant matrix elements are found to be significantly larger compared to previous results from the MIT bag model. Future calculations are needed to fully account for systematic uncertainties associated with discretization and finite-volume effects but are not expected to significantly affect this conclusion.

Tables

Table 1. Table 1: Summary of operator properties and relations to notations used in other papers. The last column shows 1-loop QCD anomalous dimensions of the operators (see Sec.  V ).

$Q_I$	Ref. Syritsyn et al. (2016)	Ref. Rao and Shrock (1982)	Ref. Grojean et al. (2018)	${(I,I_3)}_R\otimes {(I,I_3)}_L$	${\gamma }_{𝒪}^{\left(0\right)}$
$-\frac{3}{4}{Q}_1$	$\left[{\left(RRR\right)}_{\mathrm{𝟏}}\right]$	$3{𝒪}_{\left\{RR\right\}R}^3={𝒪}_{\left\{RR\right\}R}^2-{𝒪}_{\left\{RR\right\}R}^1$	$12{𝒪}_1$	${(1,-1)}_R\otimes {(0,0)}_L$	$4$
$-\frac{3}{4}{Q}_2$	$\left[{\left(RR\right)}_{\mathrm{𝟏}}{L}_{\mathrm{𝟎}}\right]$	$3{𝒪}_{\left\{LR\right\}R}^3={𝒪}_{\left\{LR\right\}R}^2-{𝒪}_{\left\{LR\right\}R}^1$	$6{𝒪}_2$	${(1,-1)}_R\otimes {(0,0)}_L$	$-4$
$-\frac{3}{4}{Q}_3$	$\left[R_{\mathrm{𝟏}}{\left(LL\right)}_{\mathrm{𝟎}}\right]$	$3{𝒪}_{\left\{LL\right\}R}^3={𝒪}_{\left\{LL\right\}R}^2-{𝒪}_{\left\{LL\right\}R}^1$	$12{𝒪}_3$	${(1,-1)}_R\otimes {(0,0)}_L$	$0$
$-\frac{5}{4}{Q}_4$	$\left[{\left(RRR\right)}_{\mathrm{𝟑}}\right]$	${𝒪}_{R\left\{RR\right\}}^1+4{𝒪}_{\left\{RR\right\}R}^2$	—	${(3,-1)}_R\otimes {(0,0)}_L$	$+24$
$-Q_5^{𝒫}$	${\left[{\left(RR\right)}_{\mathrm{𝟐}}{L}_{\mathrm{𝟏}}\right]}_{\left(1\right)}$	${𝒪}_{L\left\{RR\right\}}^1$	$-4{𝒪}_4^{𝒫}$	${(2,-2)}_R\otimes {(1,1)}_L$	$+12$
$\frac{1}{4}{Q}_6^{𝒫}$	${\left[{\left(RR\right)}_{\mathrm{𝟐}}{L}_{\mathrm{𝟏}}\right]}_{\left(2\right)}$	${𝒪}_{\left\{LR\right\}R}^2$	$-2{𝒪}_5^{𝒫}$	${(2,-1)}_R\otimes {(1,0)}_L$	$+12$
$\frac{3}{4}{Q}_7^{𝒫}$	${\left[{\left(RR\right)}_{\mathrm{𝟐}}{L}_{\mathrm{𝟏}}\right]}_{\left(3\right)}$	${𝒪}_{R\left\{RL\right\}}^1+2{𝒪}_{\left\{RR\right\}L}^2$	$-4{𝒪}_6^{𝒫}$	${(2,0)}_R\otimes {(1,-1)}_L$	$+12$

Table 2. Table 2: Results of two-point function fits from different time ranges: ground- and excited-state energies, reduced χ 2 / N dof superscript 𝜒 2 subscript 𝑁 dof \chi^{2}/N_{\text{dof}} , and optimal shrinkage parameters λ ∗ superscript 𝜆 \lambda^{*} . The uncertainties in individual fits are statistical. The last line shows “fit averages” with statistical and systematic uncertainties computed as described in Appendix  B.3 .

$t_{PP}^{min}$	$t_{PS}^{min}$	$t^{max}$	$N_{\text{dof}}$	$E_0$	$E_1$	${\chi }^2/N_{\text{dof}}$	${\lambda }^{\ast }$
6	4	13	12	0.578(23)	1.23(27)	0.50	0.14
6	6	13	10	0.556(22)	1.11(15)	0.42	0.15
6	5	13	11	0.560(24)	1.13(21)	0.40	0.14
5	5	13	12	0.566(20)	1.26(9)	0.40	0.13
7	5	13	13	0.554(69)	0.98(43)	0.42	0.15
Weighted Ave				0.565(24)(8)	1.21(15)(65)

Table 3. Table 3: Fit ranges, bare matrix element results and uncertainties in lattice units, reduced χ 2 superscript 𝜒 2 \chi^{2} showing goodness-of-fit, and optimal shrinkage parameters used for each three-point function fit for the electroweak-singlet operators Q 1 subscript 𝑄 1 Q_{1} , Q 2 subscript 𝑄 2 Q_{2} , Q 3 subscript 𝑄 3 Q_{3} , and Q 5 subscript 𝑄 5 Q_{5} .

${\tau }_P^{min}$	${\tau }_S^{min}$	$t_{\text{sep}}^{max}$	$N_{\text{dof}}$	${ℳ}_1^{\text{lat}}\times {10}^5$	${\chi }^2/N_{\text{dof}}$	${\lambda }^{\ast }$	${ℳ}_2^{\text{lat}}\times {10}^5$	${\chi }^2/N_{\text{dof}}$	${\lambda }^{\ast }$	${ℳ}_3^{\text{lat}}\times {10}^5$	${\chi }^2/N_{\text{dof}}$	${\lambda }^{\ast }$	${ℳ}_5^{\text{lat}}\times {10}^5$	${\chi }^2/N_{\text{dof}}$	${\lambda }^{\ast }$
6	2	13	70	-4.13(0.92)	0.25	0.77	8.50(1.07)	0.40	0.35	-5.01(0.76)	0.44	0.32	-0.098(39)	0.62	0.53
6	4	13	25	-3.81(1.78)	0.44	0.72	6.46(2.15)	0.31	0.31	-3.21(1.25)	0.40	0.29	-0.063(45)	0.53	0.41
6	3	13	45	-3.85(1.07)	0.30	0.76	8.24(1.52)	0.34	0.31	-4.09(0.95)	0.47	0.30	-0.068(38)	0.54	0.50
5	3	13	51	-4.09(0.92)	0.28	0.75	8.61(1.06)	0.34	0.29	-4.50(0.67)	0.44	0.29	-0.077(22)	0.54	0.47
7	3	13	40	-3.87(1.13)	0.34	0.76	8.13(1.32)	0.37	0.32	-4.05(1.00)	0.50	0.31	-0.069(32)	0.55	0.53
Weighted Ave				-3.99(1.08)(0.13)			8.28(1.29)(0.54)			-4.37(0.86)(0.52)			-0.075(32)(10)

Table 4. Table 4: Summary of renormalization constants from fits in different q 𝑞 q ranges and with 1- and 2-loop QCD running, with statistical uncertainties and uncorrelated χ 2 superscript 𝜒 2 \chi^{2} -values. The first two columns show the fit ranges and the perturbative orders of QCD matching. The last two rows show the w m subscript 𝑤 𝑚 w_{m} -weighed ( 112 ) final values with statistical and systematic uncertainties for Z I S ​ I subscript superscript 𝑍 𝑆 𝐼 𝐼 Z^{SI}_{I}  ( 74 ) and the conversion coefficients between the lattice bare values and the MS ¯ ​ ( 2 ​ GeV ) ¯ MS 2 GeV \overline{\text{MS}}(2\,\text{GeV}) scheme.

$p[\text{GeV}$	$Z^{\text{pert}}$	ndof	$Q_1$	${\chi }^2$	$Q_2$	${\chi }^2$	$Q_3$	${\chi }^2$	$Q_4$	${\chi }^2$	$Q_5$	${\chi }^2$
1.6:3.5	1L	28	0.425(30)	31.7	0.378(36)	25.6	0.369(39)	33.8	0.615(55)	22.3	0.509(47)	27.4
	2L	28	0.432(31)	32.0	0.380(37)	25.5	0.372(39)	33.9	0.646(58)	23.0	0.519(48)	27.7
1.6:4.0	1L	42	0.458(11)	37.8	0.403(14)	33.3	0.421(11)	47.1	0.605(14)	28.2	0.526(11)	35.6
	2L	42	0.471(11)	38.4	0.405(14)	33.2	0.426(11)	47.9	0.650(15)	29.6	0.544(11)	37.4
1.6:4.5	1L	56	0.462(10)	69.3	0.409(08)	68.4	0.427(08)	83.8	0.622(14)	50.2	0.530(07)	62.9
	2L	56	0.477(10)	72.1	0.411(08)	68.3	0.433(08)	86.1	0.673(15)	55.0	0.551(08)	68.5
2.0:3.5	1L	22	0.508(75)	24.6	0.469(61)	19.0	0.434(68)	28.9	0.763(91)	17.1	0.602(79)	21.9
	2L	22	0.523(76)	24.6	0.471(62)	19.0	0.440(69)	28.9	0.826(98)	17.1	0.626(82)	21.7
2.0:4.0	1L	36	0.476(15)	30.7	0.433(17)	25.3	0.452(14)	38.0	0.638(32)	24.6	0.559(24)	29.2
	2L	36	0.494(16)	30.7	0.435(17)	25.3	0.460(14)	38.4	0.698(35)	24.7	0.586(25)	29.6
2.0:4.5	1L	50	0.477(12)	60.4	0.433(10)	55.7	0.451(10)	70.4	0.656(21)	43.5	0.556(13)	53.1
	2L	50	0.495(13)	62.1	0.435(10)	55.6	0.459(11)	71.7	0.720(23)	45.6	0.585(14)	55.9
$Z_I^{SI}\left({\mu }_0\right)$ $\left({\delta }^{\text{stat}}\right)\left({\delta }^{\text{sys}}\right)$			0.471(15)(15)		0.420(17)(17)		0.437(16)(20)		0.644(22)(35)		0.543(16)(19)
$C_I^{\overline{\text{MS}}\leftarrow \text{lat}}$ $\left({\delta }^{\text{stat}}\right)\left({\delta }^{\text{sys}}\right)$			0.433(14)(14)		0.429(17)(17)		0.425(15)(20)		0.520(18)(28)		0.527(15)(19)

Table 5. Table 5: Matrix element results for the chiral basis operators with independent non-zero matrix elements in the isospin limit. The second column shows the renormalized matrix elements at a scale of 2 GeV and total uncertainty including statistical and systematic uncertainties from the bare matrix elements and non-perturbative renormalization factor added in quadrature. Renormalized results use the MS ¯ ¯ MS \overline{\text{MS}} scheme with N f = 4 subscript 𝑁 𝑓 4 N_{f}=4 active quark flavors and are obtained through nonperturbative RI-MOM renormalization and perturbative matching to MS ¯ ¯ MS \overline{\text{MS}} . The third column shows the corresponding MS ¯ ¯ MS \overline{\text{MS}} renormalized matrix elements and uncertainties after renormalization group evolution from 2 GeV to a higher scale of 700 TeV. The fourth and fifth columns show comparisons with the results of the same matrix elements in the MIT bag model from Ref.  Rao and Shrock ( 1982 ) as described in the main text.

Operator	${ℳ}_I^{\overline{\text{MS}}}\left(2\text{ GeV}\right),$	${ℳ}_I^{\overline{\text{MS}}}\left(700\text{ TeV}\right),$	$\frac{{ℳ}_I^{\overline{\text{MS}}}\left(2\text{ GeV}\right)}{\text{MIT bag A}}$	$\frac{{ℳ}_I^{\overline{\text{MS}}}\left(2\text{ GeV}\right)}{\text{MIT bag B}}$
$Q_1$	$-46\left(13\right)\times {10}^{-5}{\text{ GeV}}^6$	$-26\left(7\right)\times {10}^{-5}{\text{ GeV}}^6$	4.2	5.2
$Q_2$	$95\left(17\right)\times {10}^{-5}{\text{ GeV}}^6$	$144\left(26\right)\times {10}^{-5}{\text{ GeV}}^6$	7.5	8.7
$Q_3$	$-50\left(12\right)\times {10}^{-5}{\text{ GeV}}^6$	$-47\left(11\right)\times {10}^{-5}{\text{ GeV}}^6$	5.1	6.1
$Q_5$	$-1.06\left(48\right)\times {10}^{-5}{\text{ GeV}}^6$	$-0.23\left(10\right)\times {10}^{-5}{\text{ GeV}}^6$	-0.84	1.6

Equations

$$\begin{array}{rl}{\mathcal{O}}_{{\chi }_1{\chi }_2{\chi }_3}^1& =(u_i^{T}C{P}_{{\chi }_1}{u}_j)(d_k^{T}C{P}_{{\chi }_2}{d}_l)(d_m^{T}C{P}_{{\chi }_3}{d}_n)T_{\{ij\}\{kl\}\{mn\}}^{\text{(symm)}},\\ {\mathcal{O}}_{{\chi }_1{\chi }_2{\chi }_3}^2& =(u_i^{T}C{P}_{{\chi }_1}{d}_j)(u_k^{T}C{P}_{{\chi }_2}{d}_l)(d_m^{T}C{P}_{{\chi }_3}{d}_n)T_{\{ij\}\{kl\}\{mn\}}^{\text{(symm)}},\\ {\mathcal{O}}_{{\chi }_1{\chi }_2{\chi }_3}^3& =(u_i^{T}C{P}_{{\chi }_1}{d}_j)(u_k^{T}C{P}_{{\chi }_2}{d}_l)(d_m^{T}C{P}_{{\chi }_3}{d}_n)T_{[ij][kl]\{mn\}}^{\text{(asym)}}\end{array}$$

$$T_{\{ij\}\{kl\}\{mn\}}^{\text{(symm)}}$$

$$T_{[ij][kl]\{mn\}}^{\text{(asym)}}$$

$$C={\gamma }_2{\gamma }_4=C^{\ast }=-C^T=-C^{†},$$

$$(q_1^{T}C{q}_2{)}_{\chi }^{S,A}=(q_1^{T}C{P}_{\chi }{q}_2{)}^{S,A}=\pm (q_2^{T}C{P}_{\chi }{q}_1{)}^{S,A}$$

$$(\tilde{\psi }\psi {)}_{\chi }^A,(\tilde{\psi }{\tau }^a\psi {)}_{\chi }^S,$$

$$\begin{array}{rl}{Q}_1& =(\tilde{\psi }\psi {)}_R^{A_1}(\tilde{\psi }\psi {)}_R^{A_2}(\tilde{\psi }{\tau }^{+}\psi {)}_R^{S_3}{T}^{A_1{A}_2{S}_3}=-4{\mathcal{O}}_{RRR}^3,\\ Q_2& =(\tilde{\psi }\psi {)}_L^{A_1}(\tilde{\psi }\psi {)}_R^{A_2}(\tilde{\psi }{\tau }^{+}\psi {)}_R^{S_3}{T}^{A_1{A}_2{S}_3}=-4{\mathcal{O}}_{LRR}^3,\\ Q_3& =(\tilde{\psi }\psi {)}_L^{A_1}(\tilde{\psi }\psi {)}_L^{A_2}(\tilde{\psi }{\tau }^{+}\psi {)}_R^{S_3}{T}^{A_1{A}_2{S}_3}=-4{\mathcal{O}}_{LLR}^3,\end{array}$$

$$\begin{split}Q_{4}&=\Big{[}(\tilde{\psi}\tau^{3}\psi)^{S_{1}}_{R}\,(\tilde{\psi}\tau^{3}\psi)^{S_{2}}_{R}-\frac{1}{5}(\tilde{\psi}\tau^{a}\psi)^{S_{1}}_{R}\,(\tilde{\psi}\tau^{a}\psi)^{S_{2}}_{R}\Big{]}\,(\tilde{\psi}\tau^{+}\psi)^{S_{3}}_{R}\,T^{S_{1}S_{2}S_{3}}=-\frac{4}{5}{\mathcal{O}}^{1}_{RRR}-\frac{16}{5}{\mathcal{O}}^{2}_{RRR}\,,\end{split}$$

$$\begin{split}Q_{5}&=(\tilde{\psi}\tau^{-}\psi)^{S_{1}}_{R}\,(\tilde{\psi}\tau^{+}\psi)^{S_{2}}_{L}\,(\tilde{\psi}\tau^{+}\psi)^{S_{3}}_{L}\,T^{S_{1}S_{2}S_{3}}={\mathcal{O}}^{1}_{RLL}\,,\\ Q_{6}&=(\tilde{\psi}\tau^{3}\psi)^{S_{1}}_{R}\,(\tilde{\psi}\tau^{3}\psi)^{S_{2}}_{L}\,(\tilde{\psi}\tau^{+}\psi)^{S_{3}}_{L}\,T^{S_{1}S_{2}S_{3}}=-4{\mathcal{O}}^{2}_{RLL}\,,\\ Q_{7}&=\Big{[}(\tilde{\psi}\tau^{3}\psi)^{S_{1}}_{L}\,(\tilde{\psi}\tau^{3}\psi)^{S_{2}}_{L}-\frac{1}{3}(\tilde{\psi}\tau^{a}\psi)^{S_{1}}_{L}\,(\tilde{\psi}\tau^{a}\psi)^{S_{2}}_{L}\Big{]}\,(\tilde{\psi}\tau^{+}\psi)^{S_{3}}_{R}\,T^{S_{1}S_{2}S_{3}}=-\frac{4}{3}{\mathcal{O}}^{1}_{LLR}-\frac{8}{3}{\mathcal{O}}^{2}_{LLR}\,.\end{split}$$

$$\begin{array}{r}{Q}_I^R(\mu )=Z_{IJ}^R(\mu )Q_I^{\text{lat}}.\end{array}$$

$$(I,I_3{)}_L\otimes (I,I_3{)}_R\to (I,I_3{)}_{L+R}=(1,-1).$$

$$⟨\bar{n}|Q_4|n⟩{|}_{m_u=m_d}=0.$$

$$Q_5\sim (2,-2{)}_L\otimes (1,1{)}_R,Q_6\sim (2,-1{)}_L\otimes (1,0{)}_R,Q_7\sim (2,0{)}_L\otimes (1,-1{)}_R,$$

$$⟨I,I_3|{\hat{I}}_{+}|I,I_3-1⟩=\sqrt{(I+I_3)(I-I_3+1)}=⟨I,I_3-1|{\hat{I}}_{-}|I,I_3⟩$$

$$\displaystyle\mathbf{3}_{R}:\;\left(\begin{array}[]{r}(u^{T}Cu)^{S}_{R}\\ \sqrt{2}(u^{T}Cd)^{S}_{R}\\ (d^{T}Cd)^{S}_{R}\end{array}\right)\sim\left(\begin{array}[]{rr}1,&+1\\ 1,&0\\ 1,&-1\end{array}\right)_{R}\,,$$

$$\displaystyle\mathbf{5}_{L}:\;\left(\begin{array}[]{c}(u^{T}Cu)^{S}_{L}(u^{T}Cu)^{S}_{L}\\ 2(u^{T}Cu)^{S}_{L}(u^{T}Cd)^{S}_{L}\\ \sqrt{\frac{2}{3}}\big{[}(u^{T}Cu)^{S}_{L}(d^{T}Cd)^{S}_{L}+2(u^{T}Cd)^{S}_{L}(u^{T}Cd)^{S}_{L}\big{]}\\ 2(u^{T}Cd)^{S}_{L}(d^{T}Cd)^{S}_{L}\\ (d^{T}Cd)^{S}_{L}(d^{T}Cd)^{S}_{L}\end{array}\right)\sim\left(\begin{array}[]{rr}2,&+2\\ 2,&+1\\ 2,&0\\ 2,&-1\\ 2,&-2\end{array}\right)_{L}\,.$$

$$⟨\bar{n}|Q_5|n⟩=⟨\bar{n}|Q_6|n⟩=-\frac{3}{2}⟨\bar{n}|Q_7|n⟩,$$

$$\langle\bar{n}|{\mathcal{O}}^{1}_{RLL}|n\rangle=(-4)\langle\bar{n}|{\mathcal{O}}^{2}_{RLL}|n\rangle=(+2)\langle\bar{n}|\big{(}{\mathcal{O}}^{1}_{LLR}+2{\mathcal{O}}^{2}_{LLR}\big{)}|n\rangle\,.$$

$$Q_I^{\mathcal{P}}$$

$$Q_I^{\mathcal{C}}$$

$$Q_I^{\mathcal{T}}$$

$$Q_{I}^{\dagger}=[Q_{I}]_{\psi\leftrightarrow\bar{\psi},\,L\leftrightarrow R}=-\eta_{C}^{*6}\eta_{P}^{*6}\big{(}{\mathcal{C}}{\mathcal{P}}\big{)}\,Q_{I}\,\big{(}{\mathcal{C}}{\mathcal{P}})^{-1}\,.$$

$${\left(⟨{\bar{n}}_{+1/2}^{(-)}|Q_I|n_{+1/2}^{(+)}⟩\right)}^{\ast }=⟨n_{+1/2}^{(+)}|Q_I^{†}|{\bar{n}}_{+1/2}^{(-)}⟩$$

$$\displaystyle\quad=-\eta_{C}^{*6}\eta_{P}^{*6}\left<{n^{(+)}_{+1/2}}\right|\big{(}{\mathcal{C}}{\mathcal{P}}\big{)}\,Q_{I}\,\big{(}{\mathcal{C}}{\mathcal{P}})^{-1}\left|{\bar{n}^{(-)}_{+1/2}}\right>=\left<{\bar{n}^{(-)}_{+1/2}}\right|Q_{I}\left|{n^{(+)}_{+1/2}}\right>\,.$$

$$\begin{array}{rl}⟨{\bar{n}}_{+1/2}^{(-)}|Q_I|n_{+1/2}^{(+)}⟩& =⟨{\bar{n}}_{+1/2}^{(-)}|{\mathcal{P}}^{-1}{Q}_I^{\mathcal{P}}\mathcal{P}|n_{+1/2}^{(+)}⟩\\ & =-{\eta }^{\ast 6}⟨{\bar{n}}_{+1/2}^{(-)}|Q_I^{\mathcal{P}}|n_{+1/2}^{(+)}⟩,\end{array}$$

$$⟨{\bar{n}}_{\pm 1/2}^{(-)}|Q_I|n_{\pm 1/2}^{(+)}⟩=⟨{\bar{n}}_{\pm 1/2}^{(-)}|{\mathcal{T}}^{-1}{Q}_I^{\mathcal{T}}\mathcal{T}|n_{\pm 1/2}^{(+)}⟩$$

$$=(\mp {\eta }_T^{\ast 3})⟨{\bar{n}}_{\mp 1/2}^{(-)}|Q_I^{\mathcal{T}}|n_{\mp 1/2}^{(+)}⟩(\mp {\eta }_T^{\ast 3})=⟨{\bar{n}}_{\mp 1/2}^{(-)}|Q_I|n_{\mp 1/2}^{(+)}⟩.$$

$${\mathcal{M}}_I=⟨{\bar{n}}_{+1/2}^{(-)}|Q_I|n_{+1/2}^{(+)}⟩=⟨{\bar{n}}_{-1/2}^{(-)}|Q_I|n_{-1/2}^{(+)}⟩,$$

$${\mathcal{M}}_I^{\ast }={\mathcal{M}}_I,{\mathcal{M}}_I^P=-{\mathcal{M}}_I.$$

$${\mathcal{M}}_4=0,{\mathcal{M}}_5={\mathcal{M}}_6=-\frac{3}{2}{\mathcal{M}}_7,$$

$${\mathcal{L}}_{n\text{-}\overline{n}}=\sum _{I(\mathcal{P})}\frac{{\tilde{C}}_I^{(\mathcal{P})}(\mu )}{{\Lambda }_{\text{BSM}}^5}{\left(\frac{v^2}{{\Lambda }_{\text{BSM}}^2}\right)}^{I_L[Q_I^{(\mathcal{P})}]}{Q}_I^{(\mathcal{P})}(\mu )+\text{h.c.},$$

$${\mathcal{M}}_{n\text{-}\overline{n}}=\frac{1}{\Lambda_{\text{BSM}}^{5}}\Big{[}\sum_{I=1,2,3}\widetilde{C}_{I}{\mathcal{M}}_{I}+\frac{v^{2}}{\Lambda_{\text{BSM}}^{2}}\sum_{I=1,2,3,5}\widetilde{C}_{I}^{\mathcal{P}}{\mathcal{M}}_{I}^{\mathcal{P}}+\frac{v^{4}}{\Lambda_{\text{BSM}}^{4}}\widetilde{C}_{5}{\mathcal{M}}_{5}\Big{]}\,,$$

$$\begin{array}{rl}{G}_{nn(\sigma )}^{J{J}^{\prime }}(t>0)& =\sum _{\mathbf{x}}⟨\text{vac}|n_{\sigma }^{(+)J^{\prime }}(\mathbf{x},t)n_{\sigma }^{(+)J†}(0)|\text{vac}⟩={\Gamma }_{\alpha {\alpha }^{\prime }}^{\sigma (+)}\sum _{\mathbf{x}}⟨n_{{\alpha }^{\prime }}^{J^{\prime }}(\mathbf{x},t){\bar{n}}_{\alpha }^J(0)⟩\end{array},$$

Full text

I Introduction

In the contemporary theory of particles and fields, there is no fundamental reason for baryon number $B$ to be conserved. Quantum effects in the Standard Model (SM) can lead to $B$ violation, and at temperatures above the electroweak phase transition sphaleron processes can efficiently convert baryons into antileptons while preserving $(B-L)$, where $L$ is lepton number. Low-temperature $B$-violating effects have not been observed experimentally, and their existence would have significant implications for the stability of nuclear matter. However, the observed baryon-antibaryon asymmetry of the universe cannot be explained within the SM, which fulfills Sakharov’s conditions for baryogenesis Sakharov (1967) but does not contain enough baryon number and CP violation to reproduce the observed baryon asymmetry of the universe Cohen et al. (1993); Rubakov and Shaposhnikov (1996); Dine and Kusenko (2003); Canetti et al. (2012). Moreover, while $(B-L)$ symmetry is preserved in the SM, it likely has to be violated in its extensions (BSM theories) aimed at explaining baryogenesis, since electroweak sphaleron transitions would otherwise “wash out” any net baryon number generated by $(B-L)$-conserving interactions in the early universe.

Baryon number violation might be experimentally observed in proton decays Miura (2016) or neutron-antineutron oscillations Baldo-Ceolin et al. (1994); Chung et al. (2002); Abe et al. (2015); Bergevin (2011). The implications of these two hypothetical processes are fundamentally different: proton decay changes baryon number by $|\Delta B|=1$ unit and involves (anti)leptons, while neutron-antineutron oscillations change baryon number by $|\Delta B|=2$ units and do not involve leptons. Proton decay, even if observed, does not necessarily violate $(B-L)$ and may be insufficient to explain baryogenesis.

Despite decades of searches, neither process has been observed, constraining the strength of $B$-violating interactions. In particular models of baryogenesis, this may require higher level of $CP$ violation, which is in turn constrained by searches for the electric dipole moments of neutrons, nuclei, and atoms. However, excluding theories of baryogenesis using results from these experiments requires knowledge of nucleon matrix elements of $B$- and $CP$-violating effective interactions expressed in terms of fundamental fields, quarks and gluons. For neutron-antineutron transitions, these calculations have previously been performed using nucleon models Rao and Shrock (1982). Modern lattice QCD methods permit model-independent calculation of these matrix elements. This paper reports the first completely nonperturbative calculation of the neutron-antineutron transition matrix elements computed in lattice QCD with physical quark masses and chiral symmetry. In particular, we find that lattice QCD calculations result in substantially larger $n\text{-}\overline{n}$ matrix elements compared to nucleon model calculations. Our findings imply that $n\text{-}\overline{n}$ oscillation experiments should observe 1-2 orders of magnitude more oscillation events than was previously expected for the same BSM physics parameters.

This paper describes in detail our methodology for computing neutron-antineutron matrix elements of operators changing baryon number by $|\Delta B|=2$ units, which have already been reported in a short publication Rinaldi et al. (2019). In particular, the operator definitions, symmetry properties of their matrix elements, and their impact on phenomenology within $SU(2)_{L}\times U(1)$-symmetric extensions are discussed in Sec. II. The setup for our calculation of these matrix elements on a lattice is described in Sec. III. Extraction of ground-state matrix elements from lattice correlation functions and analysis of potential excited state contaminations are performed in Sec. IV. Nonperturbative renormalization and matching to the $\overline{\text{MS}}$ scheme are described in Sec. V. The final results for $n\text{-}\overline{n}$ matrix elements and their uncertainties are provided in Sec. VI. In Section VII, we discuss briefly the impact of our results in light of other potential sources of systematic uncertainties that are not controlled in our present calculation.

II Effective $n\text{-}\overline{n}$ interactions

II.1 Chiral basis of $n\text{-}\overline{n}$ operators

A complete basis of color-singlet, electrically-neutral six-quark operators with $uudddd$ flavor structure can be constructed from operators of the form Chang and Chang (1980); Kuo and Love (1980); Rao and Shrock (1982, 1984); Caswell et al. (1983)

[TABLE]

where quark spinor indices are implicitly contracted in the parentheses, the $P_{L,R}=\frac{1}{2}(1\mp\gamma_{5})$ are chiral projectors, and the quark color tensors $T$ are

[TABLE]

with $S_{i}$, $A_{i}$ standing for the symmetrized and antisymmetrized pairs of color indices, respectively. These operators are identical in Euclidean and Minkowski spaces with the charge-conjugation spin matrix $C$,111 To avoid confusion, throughout the paper we use Euclidean $\gamma$-matrices $(\vec{\gamma},\gamma_{4})_{\mathcal{E}uc}=(\vec{\gamma},\gamma_{4})_{\mathcal{E}uc}^{\dagger}=(-i\vec{\gamma},\gamma_{0})_{{\mathcal{M}}}$ satisfying $\gamma_{\mu}^{\dagger}=\gamma_{\mu}$.

[TABLE]

that satisfies the usual condition $C\gamma_{\mu}C^{\dagger}=-\gamma_{\mu}^{T}$. Operators involving vector diquarks $(q^{T}CP_{\chi}\gamma_{\mu}q)$ or tensor diquarks $(qCP_{\chi}\sigma_{\mu\nu}q)$ are redundant and can be related to linear combinations of the operators in Eq. (1) by spin Fierz relations. The two choices of chirality for each $(qCP_{\chi}q)$ diquark above in ${\mathcal{O}}^{1,2,3}$ provide an overcomplete basis of 18 operators. Fierz relations ${\mathcal{O}}^{2}_{\chi\chi\chi^{\prime}}-{\mathcal{O}}^{1}_{\chi\chi\chi^{\prime}}=3{\mathcal{O}}^{3}_{\chi\chi\chi^{\prime}}$ reduce the number of independent operators to 14.222 These Fierz relations are valid in four spacetime dimensions but are violated in dimensional regularization at two-loop order Buchoff and Wagman (2016). The $\overline{\text{MS}}$ scheme defined in Ref. Buchoff and Wagman (2016) includes evanescent operator counterterms that ensure that renormalized matrix elements obey these Fierz relations. Provided that matching between BSM theory and SM effective operators is consistently performed in this $\overline{\text{MS}}$ scheme or is performed at a high enough scale that one-loop QCD corrections are negligible, these Fierz relations can be assumed for $\overline{\text{MS}}$ renormalized matrix elements.

All 14 independent effective six-quarks operators are electrically neutral and change the baryon number by $\Delta B=-2$ units. However, they are not independent under isospin symmetry transformations. The electroweak (EW) symmetry $SU(2)_{L}\times U(1)_{Y}$ requires that all interactions are $SU(2)_{L}$-singlet, which may be achieved with with additional factors of the Higgs field (see Sec.II.5). Furthermore, since the chiral symmetry $SU(2)_{L}\otimes SU(2)_{R}$ is preserved exactly in the massless perturbation theory, and preserved with good precision on a lattice with chiral fermions, it is more convenient to use a basis made of operators having definite values of chiral $L,R$-isospin.

The operators in Eq. (1) are built from color-symmetric ($\mathbf{6}_{c}$) and antisymmetric ($\mathbf{\bar{3}}_{c}$) chiral diquarks, which can be denoted as

[TABLE]

where $q_{1,2}=u$ or $d$ and the relative signs upon quark permutation come from their anticommutation, $C^{T}=-C$, and color (anti)symmetry. Using the isospin doublet $\psi=(u,d)$ and its conjugate $\tilde{\psi}=(\psi^{T}C\,i\tau^{2})$, the chiral isoscalar and isovector diquarks can be written as

[TABLE]

where $\tau^{a}$ are the Pauli matrices, $[\tau^{a},\tau^{b}]=2i\epsilon^{abc}\tau^{c}$. The details of isospin classification were given in Ref. Buchoff and Wagman (2016), and here we list only the chiral-basis operators and their relation to the conventional basis (1). All the $SU(2)_{L}$-singlet operators can be constructed from some $R$-diquarks and $L$-isoscalar diquarks, resulting in three operators belonging to the $(\mathbf{1}_{L},\mathbf{3}_{R})$ irreducible representation of the chiral isospin,

[TABLE]

and one $(\mathbf{1}_{L},\mathbf{7}_{R})$ operator

[TABLE]

where $\tau^{\pm}=\frac{1}{2}(\tau^{1}\pm i\tau^{2})$. The remaining 10 independent $n\text{-}\overline{n}$ transition operators are not $SU(2)_{L}$ singlets. Of these additional operators, three belong to the $(\mathbf{5}_{L},\mathbf{3}_{R})$ irreducible representation,

[TABLE]

The remaining seven independent operators $Q_{1}^{\mathcal{P}},\cdots,Q_{7}^{\mathcal{P}}$ are obtained from $Q_{1},\cdots,Q_{7}$ by parity transformation discussed below (35) and belong to the $(\mathbf{3}_{L},\mathbf{1}_{R})$, $(\mathbf{7}_{L},\mathbf{1}_{R})$, and $(\mathbf{3}_{L},\mathbf{5}_{R})$ irreducible representations. The operators $Q_{1},\cdots,Q_{7},Q_{1}^{P},\cdots,Q_{7}^{P}$ form a complete basis of 14 linearly independent $SU(3)_{C}\times U(1)_{EM}$-invariant dimension-9 operators with baryon number $\Delta B=-2$ and isospin $\Delta I_{3}=-1$333 The isospin of operators $\Delta I_{Q}$ is defined here as $[Q,\vec{I}]=\Delta\vec{I}_{Q}Q$, leading to the selection rule $I_{i}-I_{f}=\Delta I_{Q}$ for the isospins of initial and final states. .

Isospin properties of the $n\text{-}\overline{n}$ operators are summarized in Tab. 1, together with relations to notations used in other papers. In the following sections, we will discuss nucleon matrix elements only of the operators $Q_{1,2,3,5}$, and the other matrix elements can be easily obtained using symmetries discussed below.

II.2 Operator mixing

In this work, we study lattice regularized operators that have to be nonperturbatively renormalized and then perturbatively converted to the $\overline{\text{MS}}$ scheme using the one-loop matching results of Ref. Buchoff and Wagman (2016), as described in Sec. V,

[TABLE]

The renormalization matrix $Z_{IJ}^{R}$ takes especially simple form in the “chiral basis” consisting of elements $Q_{I=1\ldots 7}^{({\mathcal{P}})}$, because they belong to different chiral multiplets and cannot mix with each other due to chiral $SU(2)_{L}\times SU(2)_{R}$ symmetry of massless QCD. Although some chiral representations appear in Tab. 1 more than once, they are actually also prevented from mixing. Specifically, operators $Q^{({\mathcal{P}})}_{5,6,7}$ consist of different components (“rows”) of chiral $\mathbf{3}$- and $\mathbf{5}$-multiplets, and transform differently under $SU(2)_{L}\times SU(2)_{R}$. Operators $Q^{({\mathcal{P}})}_{1,2,3}$ cannot mix with each other for a more subtle reason. Even though they belong to the same chiral representation, they contain different numbers of left- and right-handed diquarks. While the $U(1)_{A}$ symmetry is violated in QCD by the ABJ anomaly, operators $Q_{1,2,3}$ do not mix in perturbative QCD because perturbative gluon exchanges preserve the $U(1)_{A}$ transformation properties of external quark fields in their respective Green’s functions. At the diagram level, there are only quark (and no antiquark) external fields, which cannot be contracted into closed loops, and thus penguin-like diagrams do not appear. This point is discussed and illustrated by an explicit two-loop perturbative calculation in Ref. Buchoff and Wagman (2016).

In order to avoid mixing of renormalized operators, one has to define renormalization matrix $Z_{IJ}^{R}$ in a scheme respecting chiral symmetry, such as $\overline{\text{MS}}$, and perform perturbative matching calculations in massless QCD. Likewise, to avoid mixing of bare lattice operators $Q_{I}^{\text{lat}}$, chiral symmetry must be preserved in lattice QCD regularization, which requires [Möbius] domain wall ([M]DWF) or overlap fermion discretization. The MDWF action that we use in this work has been shown to have good chiral properties Blum et al. (2016) (see Sec. III), and our lattice results may be safely matched to perturbative QCD in the UV regime. Finally, nonperturbative effects such as spontaneous chiral symmetry breaking and $U(1)_{A}$-violating topological fluctuations (instantons) in the QCD vacuum could lead to operator mixing in nonperturbative renormalization (NPR). Mixing can be also induced by the light quark masses and residual chiral symmetry violation. However, as we study NPR numerically in Sec. V, we find that this mixing is negligible ($\approx O(10^{-3})$) and can be safely neglected at our level of precision.

II.3 Isospin relations between matrix elements

Since the chiral symmetry of QCD is spontaneously broken $SU(2)_{L}\otimes SU(2)_{R}\to SU(2)_{L+R}$, the isospin selection rules for $n\text{-}\overline{n}$ matrix elements constrain only the total isospin $I_{L+R}$ of the effective operators $Q^{({\mathcal{P}})}_{1,\dotsm,7}$. The $n\text{-}\overline{n}$ transition changes the isospin by $\Delta I_{3}=-1$, therefore the $L,R$ isospins must add as

[TABLE]

The operator in the $(\mathbf{1}_{L},\mathbf{7}_{R})$ representation (8) with the total isospin $I_{L+R}=3$ cannot couple a neutron to an antineutron ($I_{L+R}=\pm\frac{1}{2}$) in our calculation that is performed with $SU(2)_{f}$-symmetric QCD with $m_{u}=m_{d}$, therefore

[TABLE]

Even if the isospin-breaking effects $\sim(m_{u}-m_{d})\neq 0$ are included, such $SU(2)_{f}$-violating matrix elements will be suppressed with powers of $(m_{u}-m_{d})/\Lambda_{\text{QCD}}$ relative to those of other operators.

Similarly, while the $Q_{5,6,7}$ operators introduced in Eqs. (9) are linearly independent, isospin symmetry leads to additional relations between their $n\text{-}\overline{n}$ matrix elements that make two of them redundant. The relations between them are determined by the $(I,I_{3})_{L+R}=(1,-1)$ component in the product of their chiral factors,

[TABLE]

as well as their normalization. To find the latter, one can use $SU(2)$ ladder operators

[TABLE]

to construct the full $\mathbf{3}_{R}$ and $\mathbf{5}_{L}$ isospin multiplets starting from $(u^{T}Cu)^{S}_{R}\sim(1,+1)_{R}$ and $(u^{T}Cu)^{\{S_{1}}_{L}(u^{T}Cu)^{S_{2}\}}_{L}\sim(2,+2)_{L}$, respectively:

[TABLE]

Combining these components to construct $Q_{5,6,7}$ according to Eq. (13) yields their relative normalizations. Taking into account the Clebsch-Gordan coefficients for the projection (11), one obtains the relations between matrix elements

[TABLE]

which are also fulfilled in lattice contractions up to the machine precision. Additionally, one can check that these relations hold, e.g., for the results of the Bag-model calculation Rao and Shrock (1982) in the form

[TABLE]

II.4 ${\mathcal{C}}$, ${\mathcal{P}}$, and ${\mathcal{T}}$ relations

The discrete symmetries ${\mathcal{C}}$, ${\mathcal{P}}$, and ${\mathcal{T}}$, which are conserved in QCD, imply further relations for $n\text{-}\overline{n}$ transition matrix elements. Since the form of $n\text{-}\overline{n}$ operators is identical in Minkowski and Euclidean space, we study the relations between their matrix elements in Minkowski space but using Euclidean $\gamma$-matrix conventions. From the usual transformations for the fermion fields, we obtain ${\mathcal{C}}$,${\mathcal{P}}$,${\mathcal{T}}$-transformation properties for quark bilinears and the 6-quark operators, which are summarized in Appendix A,

[TABLE]

where $\eta_{{\mathcal{C}},{\mathcal{P}},{\mathcal{T}}}$ are arbitrary complex phases accompanying the ${\mathcal{C}},{\mathcal{P}},{\mathcal{T}}$ transformations of fermion fields. These factors and the relevance of Eq. (37) for $CP$-violating processes are discussed further in Refs. McKeen and Nelson (2016); Berezhiani and Vainshtein (2015); Gardner and Yan (2016). The conjugated operators $Q_{I}^{\dagger}$ are related to $Q_{I}$ by the ${\mathcal{C}}{\mathcal{P}}$ transformation,

[TABLE]

The ${\mathcal{C}}{\mathcal{P}}$ transformation also relates the transition matrix elements $\bar{n}\to n$ and $n\to\bar{n}$ , which can be shown to be real. For that, one has to use the transformation properties of the neutron and antineutron states (see Appendix A for the details):

[TABLE]

Parity relates $n\bar{n}$ transition matrix elements of $Q_{I}^{P}$ and $Q_{I}$,

[TABLE]

where the phase factor is complementary to that in Eq. (35). For the conventional choice $\eta_{P}=1$, it is clear that only the pseudoscalar combination $(Q_{I}-Q_{I}^{\mathcal{P}})$ has nonzero matrix elements, since $n\to\bar{n}$ transition changes parity. Note that in all the cases, the arbitrary phase factors $\eta_{C,P,T}$ arising from the transformations of $Q_{I}$ cancel with the phase factors arising from the transformations of the states.

Finally, with the help of the ${\mathcal{T}}$-reflection, one can also show that the matrix elements do not depend on the direction of the (anti)neutron spin. Using the transformation properties of the neutron and antineutron states,

[TABLE]

All spin-flip matrix elements of $Q_{I}$ are trivially zero because $Q_{I}$ are (pseudo)scalars.

Denoting the ground-state $n\bar{n}$ transition matrix elements for each $Q_{I}$ by

[TABLE]

the matrix element results derived above can be summarized as

[TABLE]

In conjunction with the results from Sec. II.3,

[TABLE]

this implies that in the isospin limit where Eq. (44) is valid, $n\bar{n}$ transition rates involving the 14 operators $Q_{I}^{({\mathcal{P}})}$ are given in terms of 4 real $n\bar{n}$ transition matrix elements ${\mathcal{M}}_{1,2,3,5}$.

II.5 $n\bar{n}$ Effective Field Theory

The $|\Delta B|=2$ effective interactions discussed above must be generated by some extension of the Standard model at yet unknown scale $\Lambda_{\text{BSM}}$. It is generally assumed that such extensions have higher symmetry, which is broken at scales below $\Lambda_{\text{BSM}}$ to the electroweak symmetry $SU(2)_{L}\times U(1)_{Y}$, and thus the effective interactions must be EW-symmetric. From the discussion above it follows that only $Q_{1,2,3}$ are $SU(2)_{L}\times U(1)_{Y}$-singlets, while $Q_{4}$, $Q_{5(67)}$, and all $Q_{I}^{{\mathcal{P}}}$ operators are not. These latter operators require additional EW-charged factors to make them EW-symmetric, which affect the power counting and result in higher suppression by the $\Lambda_{\text{BSM}}$ scale.

Such factors can be easily constructed from the Higgs field doublet $\phi$ and its conjugate $i\tau_{2}\phi^{*}$ to compensate for the $SU(2)_{L}$- and hypercharge of the operators $Q_{I}^{({\mathcal{P}})}$. The Higgs v.e.v. $v$ in unitary gauge leads to nonzero effective $n\text{-}\overline{n}$ interaction in the form

[TABLE]

where $\widetilde{C}_{I}(\mu)$ are dimensionless Wilson coefficients and $I_{L}[Q_{I}^{({\mathcal{P}})}]$ is the left-handed isospin of the operator $Q_{I}^{({\mathcal{P}})}$. In addition to Eq. (45), the full $|\Delta B|=2$ Lagrangian must also include combinations of electrically charged $|\Delta B|=2$ operators with oppositely charged Higgs fields to assure the EW symmetry above the EW scale. Such interactions can lead to $n\leftrightarrow\bar{p}$ and $p\leftrightarrow\bar{n},\bar{p}$ transitions, and the emitted charged Higgs bosons (e.g., decaying into leptons) would compensate for the change in the electric charge. These transitions are suppressed by at least one factor of $(v^{2}/\Lambda_{\text{BSM}}^{2})$. 444 Isospin breaking effects in QCD may result in suppression in powers of $(m_{u}-m_{d})/\Lambda_{\text{QCD}}$ instead of $v/\Lambda_{\text{BSM}}$, which is beyond the scope of the present paper. Also, while additional higher-dimensional operators suppressed at the same level may be constructed using field derivatives, they are less relevant and not considered here.

Using the effective Lagrangian (45) and the relations derived in the previous sections, the full $n\text{-}\overline{n}$ matrix element can be written as

[TABLE]

where ${\mathcal{M}}^{({\mathcal{P}})}_{I}$ are the nucleon matrix elements of operators $Q^{({\mathcal{P}})}_{I}$. The dimensionless low-energy constants $\widetilde{C}^{({\mathcal{P}})}_{I}(\mu)$ depend on the scale $\mu$ only logarithmically and can be computed perturbatively by using $C_{I}(\Lambda_{\text{BSM}})\sim O(1)$ given by a particular BSM scenario as an initial condition for renormalization group evolution. A nonperturbative calculation of the matrix elements ${\mathcal{M}}_{I}$ is presented in the following sections.

III Lattice setup

In this section, we fist recount the details of the lattice QCD gauge configurations and propagators used in this study, and then describe the construction of (anti)neutron correlation functions with the $n\text{-}\overline{n}$ operators. The QCD gauge field configurations were generated with the Iwasaki gauge action on a $48^{3}\times 96$ lattice and $N_{f}=2+1$ flavors of dynamical Möbius Domain Wall fermions. The fermion masses are tuned to be almost exactly at the physical point Blum et al. (2016), such that the pion mass is approximately $m_{\pi}=139.2(4)\text{MeV}$ and the scale (the lattice spacing) is $a=0.1141(3)\text{ fm}$. The residual mass $m_{res}$, which encapsulates the residual violation of chiral symmetry, is smaller than 50% of the input quark mass. The physical lattice size $L\approx 5.45\text{ fm}$ and $m_{\pi}L=3.86$ should be sufficient to suppress finite volume effects of the $n\text{-}\overline{n}$ matrix elements to a level below our target precision. In particular, according to chiral perturbation theory, these finite size effects are expected to be $\lesssim 1\%$ Bijnens and Kofoed (2017).

The three-point functions needed to evaluate the matrix elements of the operators $Q^{({\mathcal{P}})}_{I}$ require six quark propagators for the $u$ and $d$ quarks flavors; they result from Wick contractions of the six-quark operators with the (anti)neutron interpolating fields. There are no disconnected quark-loop diagrams because the operators $Q_{I}$($Q_{I}^{\dagger}$) contain only quarks (antiquarks). For both two- and three-point lattice correlation functions we compute propagators on 30 independent gauge field configurations separated by 40 molecular dynamics time steps. All the quark propagators required for a single sample are computed from a point source located at the operator insertion point, which is identified in the analysis with the origin $x_{0}=(0,0,0,0)$ using translational invariance. To reduce stochastic uncertainty, sampling of the neutron correlation functions is enhanced by all-mode-averaging Shintani et al. (2015), in which we compute 1 exact and 81 low-precision samples evenly distributed over the 4D volume on each gauge configuration. The low-precision quark propagators are computed with low-mode deflation and the conjugate gradient algorithm truncated at 250 iterations.

The propagators are contracted at the sink into intermediate baryon blocks Doi and Endres (2013); Detmold and Orginos (2013) with polarized nucleon and antinucleon quantum numbers to minimize the time spent in the contraction step of the calculation. (Anti)neutron source and sink interpolating operators are constructed with either point or Gaussian-smeared (anti)quarks and are denoted with $n^{J=P,S}$, respectively. Final contraction at the propagator source yields an (anti)neutron two-point correlation function sample with a point source at $x_{0}$. Thus, the polarized neutron two-point correlation function with zero spatial momentum for positive time $t>0$ is

[TABLE]

and, similarly, for the polarized antineutron

[TABLE]

where the polarization matrix $\Gamma^{\sigma(\pm)}=\frac{1\pm\gamma_{4}}{2}\frac{1+\sigma\gamma_{3}\gamma_{5}}{2}$ projects on the selected parity ($\pm$) and spin $\sigma=\pm\frac{1}{2}$, and the interpolating operator at the source is $J=P$ and the one at the sink can be either $J^{\prime}=P$ or $J^{\prime}=S$. Neutron/antineutron two-point functions have the spectral representation

[TABLE]

where the overlap factors $Z_{m}^{J}$ are identical for neutrons and antineutrons in either spin orientations.

The three-point functions involve two neutron or two antineutron fields to create and annihilate states with opposite baryon numbers. Using the (anti)neutron states defined in Eq. (98), one can express the three-point correlation function containing, for example, the $n\leftarrow\bar{n}$ transition matrix element (see Fig. 1):

[TABLE]

where $n^{(-)}_{-\sigma}$ is nucleon interpolating field that creates an antineutron with spin $\sigma$555 Note that spin-flip of a spinor incorporates nontrivial signs in order to satisfy $n_{-(-\sigma)}=-n_{\sigma}$ similar to ${\mathcal{T}}$ transformation, which is responsible for the relative signs of the neutron and antineutron states. , both the (anti)neutron operators are summed over the spatial coordinate to project on zero momentum. By calculating quark propagators with point sources located at the operator insertion point and momentum-projected $P$ and $S$ sinks located on all time slices, the correlation functions $G_{nQ_{I}^{\dagger}\bar{n}}^{JJ^{\prime}}$ can be accessed for all smearing combinations $PP$, $PS$, $SP$, and $SS$, any temporal separation between the source and the sink $t_{\text{sep}}=t_{1}+t_{2}$, and any operator separation from the source $\tau=t_{1}$. The same propagators are used to calculate $PP$ and $PS$ two-point correlation functions. The spectral representation for Eq. (50) analogous to Eq. (49) is given by

[TABLE]

where $(\mathcal{M}_{I})_{m^{\prime}m}=\left<m^{\prime}\vphantom{Q_{I}m}\right|Q_{I}\left|m\vphantom{m^{\prime}Q_{I}}\right>$, the ground-state matrix element of interest is $\mathcal{M}_{I}=(\mathcal{M}_{I})_{00}$, and the overlap factors $Z_{m}^{J}$, $Z_{m^{\prime}}^{J^{\prime}}$ are the same as in Eq. 49. We perform contractions for all combinations of point and smeared sources and sinks in the three-point functions to enhance the analysis of the ground and excited state matrix elements in the next section. To reduce stochastic uncertainties, we also average lattice matrix elements over the spins of the neutron and antineutron states. The specific combinations of (anti)neutron 4-spinor components in the three-point functions that give matrix elements ${\mathcal{M}}_{I}$ are

[TABLE]

where the signs correspond to the conventions listed in Appendix A. As shown in Sec. II.4, these matrix elements are real, and combining them with the conjugated ones is also used to enhance statistics following Eq. (43).

IV Analysis of matrix elements

To account for excited-state contributions, we perform two-state fits to a truncation of Eq. (50),

[TABLE]

where $\mathcal{A}_{I}^{JJ^{\prime}}$ and $\mathcal{B}_{I}^{JJ^{\prime}}$ are products of overlap factors and matrix elements involving only excited states, which are discarded in our calculation. The ground-state overlap factors $Z_{0}^{P}$ and $Z_{0}^{S}$ required to extract matrix elements of $G_{nQ_{I}^{\dagger}\bar{n}}^{JJ^{\prime}}$ can be obtained independently from fits of two-point functions $G_{nn}^{PP}$ and $G_{nn}^{PS}$ to an analogous two-state model

[TABLE]

The energies $E_{0}$ and $E_{1}$ appear in both Eq. (53) and Eq. (54), therefore fits of $G_{nQ_{I}^{\dagger}\bar{n}}$ may be simplified by fixing the state energies $E_{0}$, $E_{1}$ to values determined from fits of two-point functions $G_{nn}^{JJ^{\prime}}$. In principle, the overlaps with excited neutron states $Z_{1}^{J}$ are also determined from two-point function fits, thus the number of parameters in Eq. (53) can be reduced by factoring $\mathcal{A}_{I}^{JJ^{\prime}}$, $\mathcal{B}_{I}^{JJ^{\prime}}$ into excited-state matrix elements and overlap factors $Z_{0,1}^{J}$, of which only the latter would depend on the neutron interpolating operators. It would be possible if the two- and three-point functions were saturated by contributions only from the ground and the first excited states, or their contributions could be reliably distinguished from higher-energy states omitted from Eqs. (53,54). However, as our two-point function fits in Fig. 2 show, there are higher excited-state contributions to $G_{nn}$; in particular, there is large systematic uncertainty on $E_{1}$ in (see Tab. 2).

These considerations lead us to adopt the following fit strategy: first, a combined fit of $G_{nn}^{PP}$ and $G_{nn}^{PS}$ to Eq. (54) is used to determine the four parameters $E_{0,1}$ and $Z_{0}^{P,S}$ as summarized in Fig. 2 and Tab. 2; then, a combined fit of $G_{nQ_{I}^{\dagger}\bar{n}}^{PS}$, $G_{nQ_{I}^{\dagger}\bar{n}}^{SP}$, and $G_{nQ_{I}^{\dagger}\bar{n}}^{SS}$ to Eq. (53) is used to determine the six parameters ${\mathcal{M}}_{I}$, $\mathcal{A}_{I}^{PS}$, $\mathcal{A}_{I}^{SP}$, $\mathcal{A}_{I}^{SS}$, $\mathcal{B}_{I}^{PS}=\mathcal{B}_{I}^{SP}$, and $\mathcal{B}_{I}^{SS}$. Also, since $PP$ three-point functions would have even large excited-state contamination and $PP$ three-point/two-point ratios are not close to their plateau region for the $t_{\text{sep}}$ used here (not shown), we do not include $G_{nQ_{I}^{\dagger}\bar{n}}^{PP}$ in our analysis.

With all-mode-averaging described in Sec. III, we obtain one unbiased sample of the two- and three-point functions per gauge field configuration. The number of gauge field configurations used in this calculation $N_{\text{conf}}=30$ is not large enough to obtain nondegenerate determination of a covariance matrix for the required number of data points $31\leq K\leq 76$ included in the three-point correlator fits. Therefore, spin and parity symmetries are used to increase the effective number of unbiased samples of correlation functions. Thus, $G_{nQ^{\dagger}\bar{n}}$ and $G_{n(Q_{I}^{\mathcal{P}})^{\dagger}\bar{n}}$ with two polarizations are treated as four samples per gauge-field configuration, resulting in $N=120$ samples for each data point after all-mode-averaging bias correction. Polarized two-point functions $G_{nn(\pm 1/2)}$, $G_{\bar{n}\bar{n}(\pm 1/2)}$ are similarly combined to obtain a statistical ensemble of $N=120$ two-point functions. Although this yields an “ensemble” with $N>K$ samples, it is still not sufficient for reliable determination of covariance matrix, which typically requires $N\gtrsim K^{2}$.

For both two-point and three-point functions, finite sample-size fluctuations may make the sample covariance matrix ill-determined and lead to a numerically unstable inverse covariance matrix required for least-squares fitting. Shrinkage Stein (1956); Ledoit and Wolf (2004) has been proposed as a method of improving the condition number of covariance matrix estimates. Denoting the sample covariance matrix by $S$, the corresponding covariance matrix estimate with shrinkage is given by

[TABLE]

where $0\leq\lambda\leq 1$ is the shrinkage parameter, and $\text{diag}(S)$ is a particular “shrinkage target”. Taking any $\lambda>0$ “shrinks” the spectrum of the covariance matrix by reducing the relative size of off-diagonal correlations compared to the diagonal covariance matrix elements. This leads to a better-conditioned covariance matrix and a more robust estimate of the inverse covariance matrix used for $\chi^{2}$-minimization. Trivial $\lambda=0$ corresponds to no shrinkage, while $\lambda=1$ removes off-diagonal correlations completely, which is equivalent to an uncorrelated fit. Therefore, varying the parameter $0\leq\lambda\leq 1$ interpolates continuously between correlated (albeit with potentially poorly-determined covariance matrix) and uncorrelated fits. A standard prescription for choosing the optimal shrinkage parameter is to minimize the rms difference between $\Sigma(\lambda)$ and the true covariance matrix. A sample estimator for the optimal shrinkage parameter $\lambda^{*}$ is suggested in Ref. Ledoit and Wolf (2004) and summarized in Appendix B.1. Bootstrap covariance matrices with optimal shrinkage666 $\lambda^{*}$ are chosen to provide optimal shrinkage for the normalized sample correlation matrix as described in Appendix B.1 rather than the bootstrap covariance matrix. It is possible that finite-sample-size bias will lead to differences between the optimal shrinkage parameters for the two matrices. Since $\lambda^{*}$ defined by this prescription vanishes in the infinite-statistics limit, the bootstrap covariance matrix obtained with this choice of shrinkage parameter will provide an unbiased (but not necessarily optimal) estimate of the true covariance matrix in the infinite-statistics limit.

$\Sigma^{*}=\Sigma(\lambda^{*})$ are obtained by inserting $\lambda^{*}$ from Eq. (105) into Eq. (55) with $S$ the bootstrap covariance matrix obtained from $N_{boot}=10,000$ samples of two-point and three-point correlation functions. The effects of shrinkage on the central values, uncertainties, and goodness-of-fit of the matrix element fits described below are explored by varying $\lambda$, and the results for one choice of fit range are shown in Figs. 3-4. For all operators, the central values and the statistical uncertainties are relatively insensitive to the value of the shrinkage parameter once $\lambda>0$. The $\chi^{2}/N_{\text{dof}}$ values decrease sharply in a small region around $\lambda=0$, however they are much less sensitive for larger $\lambda$ values. In all cases, the optimal values $\lambda^{*}$ for the shrinkage parameter are found outside of the region of strong dependence of $\chi^{2}$ on $\lambda$.

The average two-point function and the corresponding bootstrap covariance matrix with optimal shrinkage are used for nonlinear $\chi^{2}$-minimization to determine $E_{0}$, $E_{1}$, $\sqrt{Z_{0}^{P}}$, and $\sqrt{Z_{0}^{S}}$. $\chi^{2}$-minimization is reduced to a two-parameter optimization problem by variable projection(VarPro) technique Golub and Pereyra (2003); Dianne P. O’Leary (2013) detailed in Appendix B.2. In VarPro, the products of overlap factors in Eq. (54) are found from a linear $\chi^{2}$-fit for particular values of $E_{0,1}$ and the solution is substituted back into $\chi^{2}$ in order to obtain a two-parameter function $\chi^{2}_{VP}(E_{0},E_{1})$, which is then minimized using nonlinear numerical methods. With these four parameters held fixed, the remaining six free parameters in the three-point function fit (53) can also be found from a linear $\chi^{2}$ fit. The parameter covariance matrix for all 10 parameters is subsequently estimated using an additional correlated bootstrap resampling and fitting of two- and three-point function data. The original bootstrap covariance matrices are used in all these fits in order to avoid the possibility of ill-conditioned covariance matrices. The bootstrap parameter covariance matrix is obtained from $N_{boot}^{\prime}=200$ correlated resampling draws. The parameter covariance matrix is diagonalized, the eigenvalues are resampled $N_{boot}^{\prime\prime}=200$ times, and the resampled eigenvalues are transformed back to the original parameter basis. Finally, the standard deviation of the resulting resampled values of ${\mathcal{M}}_{I}$ is used to define the marginalized uncertainty of ${\mathcal{M}}_{I}$ for each fit range shown in Tab. 3. An analogous procedure is used to obtain the uncertainties of $E_{0}$ and $E_{1}$ shown in Tab. 2. Results for the matrix elements from the fits that include our smallest $t_{\text{sep}}$ value are compared to the ratios of three-point to two-point functions (adjusting for proper overlap factors) in Figs. 5. In addition, the two-point function fits are compared to the corresponding effective masses in Fig. 2.

Systematic uncertainties of our analysis procedure are studied by varying the time ranges of data included in the two-state fits. Results of fits of two-point function data $G^{PP}_{nn(\bar{n}\bar{n})}(t^{min}_{PP}\leq t\leq t_{\text{sep}}^{max})$ and $G^{PS}_{nn(\bar{n}\bar{n})}(t^{min}_{PS}\leq t\leq t_{\text{sep}}^{max})$ for a variety of $t^{min}_{PP}$ and $t^{min}_{PS}$ are shown in Tab. 2. Results of corresponding fits of three-point function data $G_{nQ_{I}^{\dagger}\bar{n}}^{SS}(\tau_{min}^{P}+\tau_{min}^{S}\leq t_{\text{sep}}\leq t_{\text{sep}}^{max},\tau_{min}^{S}\leq\tau\leq t_{\text{sep}}-\tau_{min}^{S})$ and $G_{nQ_{I}^{\dagger}\bar{n}}^{PS}(2\tau_{min}^{S}\leq t_{\text{sep}}\leq t_{\text{sep}}^{max},\tau_{min}^{P}\leq\tau\leq t_{\text{sep}}-\tau_{min}^{S})$ for a variety of $\tau_{min}^{S}$ and $\tau_{min}^{P}$ are shown in Tab. 3. The data for $SP$ and $PS$ three-point correlation functions are averaged using relation $G_{nQ_{I}^{\dagger}\bar{n}}^{SP}(t_{\text{sep}},\tau)=G_{nQ_{I}^{\dagger}\bar{n}}^{PS}(t_{\text{sep}},t_{\text{sep}}-\tau)$ to reduce the number of data points in the fits. Results for bare ground-state matrix elements from different fits are in very good agreement with each other, as shown in Fig.6.

The five fit range choices shown in Tabs. 2-3 result in acceptable correlated $\chi^{2}/N_{\text{dof}}$ values in fits of two- and three-point function data. These results are combined into final estimates of ${\mathcal{M}}_{I}$ and estimates of their statistical and systematic uncertainties. Since the various fits have different $N_{\text{dof}}$, we use a weighted-averaging procedure defined in Appendix B.3. For a particular fit, the weight is a combination of the likelihood that the fit describes the data (we use its $p$-value as the likelihood proxy) and its statistical precision, to penalize both fits that fail to describe data and fits that do not constrain the relevant parameters. The same weights are used to determine the average statistical uncertainty, which ensures that including multiple similar fits will not lead to a spurious reduction in the final statistical uncertainty. The weighted mean-square difference between each fit result and the weighted average is used to define the systematic uncertainty due to arbitrariness of choice of a fit window. Applying this weighted averaging procedure to the ground-state energy $E_{0}$ of the two-point function yields the result for the nucleon mass that agrees well with the physical value,

[TABLE]

where we have used the scale-setting result $a=0.1141(3)\,\text{fm}$ from Ref. Blum et al. (2016), which has negligible uncertainty for our purposes as it is much smaller compared to other uncertainties in our calculation. Applying the same procedure to the fit results in Tab. 3 provides our final estimate of the bare matrix elements including statistical and fitting systematic uncertainties,

[TABLE]

These lattice regularized matrix elements can be related to renormalized matrix elements through NPR as described in the next section.

V Renormalization of lattice operators

Since the matrix elements of the 6-quark operators are computed on a lattice, they have to be converted to some perturbative scheme, e.g., $\overline{\text{MS}}$, before they can be used in BSM phenomenology. We calculate conversion factors between lattice-regularized operators and their perturbative definitions nonperturbatively, by computing their Green’s functions on a lattice and matching them to perturbative calculations. The operators $Q^{({\mathcal{P}})}_{I}$ are the lowest-dimension operators with $\Delta B=-2$, therefore they can only either mix with each other, or get discretization corrections from higher-dimensional operators that vanish in the continuum limit. In the chiral basis, all the 14 operators transform differently under $U(2)_{L}\otimes U(2)_{R}$ flavor symmetry, so they can mix only due to the spontaneous chiral symmetry breaking (S$\chi$SB) in QCD, non-perturbative $U(1)_{A}$ violation, or chiral symmetry violations by quark masses and discretization of the fermion action. Mixing due to quark masses and non-perturbative effects should be small if renormalization is carried out in the UV region $|p|\approx\mu\gg\{\Lambda_{QCD},m_{q}\}$ where perturbative matching is applicable. Furthermore, effects of the explicit chiral symmetry violation by the (M)DWF fermion action on a lattice are suppressed as the “residual mass” $m_{res}\lesssim m_{q}$ Blum et al. (2016), and thus are also negligible. Therefore, we don’t expect that renormalization of our results will be affected by mixing between the chiral-basis operators 777 This holds even without taking the continuum limit, since the continuum and the chiral limits can be taken separately in calculations with the (M)DWF lattice fermions, see Ref. Blum et al. (2016) and references within. .

V.1 RI-MOM amplitudes on a lattice

The lattice renormalization constants for the 6-quark operators are defined as

[TABLE]

but, as will be shown below, in the chiral-diagonal basis $|Z_{I\neq J}|\ll Z_{II}\equiv Z_{I}$, so $Q_{I}^{R}(\mu)=Z_{I}Q_{I}^{\text{lat}}$ both on a lattice and in continuum perturbation theory. The nonperturbative renormalization and mixing of the six-quark operators is computed using a variant of the RI-MOM scheme Martinelli et al. (1995) with a specific choice of momenta of the external quark states. Since the external states are not color-singlets, the gauge is fixed to the Landau gauge using the Fourier-accelerated conjugate gradient algorithm Hudspith (2015). All the operators of interest with $\Delta B=2$ and $\Delta I=1$ 888 Instead of the 6-quark, we study renormalization of the 6-antiquark operators, which is equivalent but is more natural on a lattice since it does not require conjugating quark propagators. can be represented in the generic form

[TABLE]

where $A_{i}=(\alpha_{i},a_{i})$ are the spin$\times$color indices. Then their Green’s functions with external plane-wave quark states,

[TABLE]

are computed on a lattice contracting six quark propagators computed with a point source at the operator location. The same propagators are used as for the $n\text{-}\overline{n}$ three-point correlators (see Sec. III), with the only difference that prior to the contraction the propagators are Fourier-transformed at the sink. The six-quark vertex functions are obtained by “amputating” the Green’s functions (60)

[TABLE]

where contraction in $\{B_{i}\}$ is implied, and the momentum-projected quark propagators are

[TABLE]

Note that the amputated Green’s functions (61) are not symmetric with respect to permutation of the spin$\times$color indices $A_{i}$, unlike the tree-level vertex function $\Gamma_{i}^{[A_{1}A_{2}][A_{3}\ldots A_{6}]}$ in Eq. (59). This is due to the fact that $G_{I}(\{p_{i}\})$ and $\Lambda_{I}(\{p_{i}\})$ depend on the non-equal momenta $p_{i}$ of the external fields. Such dependency would break the isospin symmetry and thus may mix operators from different chiral representations.

One must choose specific momenta for external quark fields in order to preserve the chiral isospin symmetry. The simplest choice $p_{i}=p$ would result in a large momentum $p_{\mathcal{O}}=6p$ at the operator insertion leading to large perturbative corrections in conversion to the $\overline{\text{MS}}$ scheme. To avoid that, the external quark momenta are arranged so that $\sum_{i}p_{i}=0$ and, specifically, $p_{i}=\pm p$ (see Fig. 7), where $p^{2}=\mu^{2}$ determines the scale for perturbative RI-MOM $\to\,\overline{\text{MS}}$ matching. In addition, the amputated amplitudes (61) must be averaged over permutations of the $\pm p$ momenta to enforce the symmetry with respect to the external quark states Syritsyn et al. (2016),

[TABLE]

where the factors are determined by combinatorics. All possible permutations of momenta are implicitly included by Wick contractions, and the symmetries of the color$\times$spin indices are restored automatically. Perturbative matching at the one-loop level for this particular scheme has been computed in Ref. Buchoff and Wagman (2016).

The lattice renormalization factors $Z^{\text{lat}}(p^{2})$ (58) can be computed by imposing the condition

[TABLE]

where $Z_{q}$ is the lattice quark field renormalization factor

[TABLE]

The renormalization factors $Z_{IJ}^{\text{lat}}$ can be expressed in terms of the amputated and symmetrized vertex functions $\Lambda_{I}^{\{A_{i}\}}(p)$ that are projected onto the original tree-level structures $\Gamma_{J}^{\{A_{i}\}}$,

[TABLE]

where the “metric tensor” $g_{JK}$ is diagonal in the chiral basis $\overline{Q}_{I}^{({\mathcal{P}})}$. (Approximate) chiral symmetry on a lattice is important for ensuring that $Z^{IJ}$ and $\Lambda_{IJ}$ are also (predominantly) diagonal in this basis. Deviations from the diagonal form are due to the nonzero quark mass and residual chiral symmetry breaking of the DWF discretization. The effect of symmetrization (63) is evident from the magnitude of the off-diagonal components, which is shown in the log scale as the matrix

[TABLE]

in Fig. 8 comparing the momentum permutation-averaged amplitude (63) to the one with a specific choice of momentum $p_{1}=p_{3}=p_{4}=-p_{2}=-p_{4}=-p_{6}=p$ as in Fig. 7. These data are shown for the momentum p=\frac{2\pi}{a}\big{(}\frac{11}{48},\frac{11}{48},\frac{11}{48},\frac{22.5}{96}\big{)}, which is close to a 4d diagonal direction (up to $(\pi/L)$ along the time axis due to the antiperiodic boundary conditions) and $p^{2}\approx(5\,\text{GeV})^{2}$. Therefore, we conclude that in the chiral basis the renormalization matrix $Z_{IJ}$ is diagonal, $|Z_{IJ}|/\sqrt{Z_{II}Z_{IJ}}\lesssim O(10^{-3})$, which is definitely within our target precision, and the operators $Q^{(P)}_{I}$ may be renormalized multiplicatively in our lattice calculation. Additionally, we observe that the mixing between 6-quark operators containing different numbers of $L,R$-diquarks is negligible, indicating that nonperturbative chirality-changing effects due to fluctuations of topology of the QCD vacuum do not lead to mixing in excess of the $10^{-3}$ level.

We define lattice renormalization factors in the RI-MOM scheme for the $n\text{-}\overline{n}$ operators in the chiral basis as

[TABLE]

Finally, to get rid of the quark field renormalization, we use the renormalization constant $Z_{A}$ for the local axial-vector current $A_{\mu}=\bar{q}\gamma_{\mu}\gamma_{5}q$. Using the value of $Z_{A}$ computed in Ref. Blum et al. (2016), we can compute $Z_{q}(p)$ in the RI-MOM scheme from the condition

[TABLE]

where $\langle A_{\mu}q\bar{q}\rangle$ is the amputated Green’s function for the axial current computed analogously to Eq. (61). “Scale-independent” lattice renormalization factors $Z_{\Gamma}^{SI}=Z^{\text{lat}}_{\Gamma}(p)/Z^{RI,\text{pert}}(p)$ for the vector, tensor, and scalar vertices are shown in Fig. 9.

The value of the lattice renormalization constants $Z_{I}(p)$ may depend on the orientation of the momentum $p$ with respect to the lattice axes due to discretization effects. We compute the lattice vertex functions (67) for various orientations of lattice momenta interpolating between 3d-diagonal and 4d-diagonal orientations to study these effects in the following sections.

V.2 Perturbative running

In order to convert operator normalization from the RI-MOM scheme discussed above to $\overline{\text{MS}}$, perturbative matching calculations are required. To extract lattice renormalization factors independent from the momentum subtraction point $p$, the lattice factors (70) are compared to the perturbative predictions for the RI-MOM scheme in some window $p_{\text{min}}\leq|p|\leq p_{\text{max}}$ where lattice artifacts are believed to be under control. In this section, details of relevant perturbative results are summarized.

The one-loop anomalous dimensions of the operators (7-9) were computed in Ref. Caswell et al. (1983), and the $\overline{\text{MS}}$ anomalous dimensions to the $O(\alpha_{S}^{2})$ precision together with $O(\alpha_{S})$ conversion factors were computed in Ref. Buchoff and Wagman (2016). In the chiral basis, the perturbative renormalization of the operators is diagonal (no mixing), and their independent anomalous dimensions are

[TABLE]

with the coefficients $\gamma_{I}^{(0)}$ given in Tab. 1. These anomalous dimensions are substantially different, which would complicate operator renormalization if chiral symmetry was violated by a lattice fermion action and mixing was allowed. We integrate the equations (72) together with an RG equation for the coupling constant $\alpha_{S}(\mu)$ using the 4-loop $\beta(\alpha_{S})$-function. Since our lattice QCD action has $N_{f}=2+1$ dynamical flavors, the lattice factors (70) are matched to $Z^{RI}(\mu)$ factors computed in $N_{f}=3$ perturbative QCD and the coupling constant $\alpha_{S}^{N_{f}=3}$ is matched to its physical value at $\mu\leq m_{c}$. The latter is obtained from a global fit Bethke (2009) and matched at the $m_{b,c}$ quark mass thresholds. For the reference point $\mu_{0}=2\,\text{GeV}$, its values999 The coupling constant in the RI-MOM scheme is conveniently defined to be equal to the $\overline{\text{MS}}$ coupling constant.

are $\alpha_{S}^{N_{f}=3}=0.2827$ and $\alpha_{S}^{N_{f}=4}=0.2948$.

The final results are converted to $N_{f}=4$ QCD at $\mu_{0}=2\,\text{GeV}$, again matching at the $m_{c}$ threshold. The final conversion factors from lattice to the $\overline{\text{MS}}$ scheme at scale $\mu_{0}$ are

[TABLE]

where $Z^{SI}_{I}$ is a “scale-independent” lattice renormalization factor with a reference point $\mu_{0}$ defined in the next section. The perturbative scale dependence in both the $\overline{\text{MS}}$ and RI-MOM schemes with $N_{f}=3$ and 4 flavors is shown in Fig. 10.

V.3 Fits of nonperturbative and discretization effects

With known perturbative running, we can separate scale-independent renormalization from lattice artifacts and nonperturbative effects. Correlation functions computed on a lattice are subject to discretization effects that may break rotational symmetry at short distances, which are relevant for the large momenta used in the nonperturbative renormalization. In addition, they may have nonperturbative contributions that complicate matching with perturbative calculations. Below we follow closely the analysis performed in Ref. Blossier et al. (2015) and extract the scale-invariant renormalization constants $Z_{I}^{\text{SI}}$ from a fit

[TABLE]

where $Z^{SI}(\mu_{0},a)$ is the momentum-independent lattice renormalization constant, $Z_{I}^{RI,pert}(\mu)$ is the perturbative running of $Q_{I}$ in the RI-MOM scheme and $\Delta Z_{I}^{\text{disc,NP}}$ encapsulates discretization and nonperturbative corrections. In our calculation with $O(a)$-improved action, the discretization effects must scale as O\big{(}(ap)^{2}\big{)},

[TABLE]

where we also include the hypercubic invariant $\propto p^{[4]}$ (see Fig. 11)

[TABLE]

that breaks the rotational symmetry $O(4)\to H(4)$ for $k\neq 0,1$.

Although the vertex functions (61) are computed with “exceptional” kinematics $p_{\mathcal{O}}=0$ (see Fig. 7), they do not have “pole” contributions $\propto 1/p^{2}$ because, unlike the pseudoscalar density operator that can couple to pions, the 6-quark operators $Q_{I}$ can couple only to 2-baryon ($B=2$) states with masses $M\geq 2m_{N}$. However, the nonperturbative contributions are added to Eq. (74)

[TABLE]

to account for effects of the dimension-2 gluon condensate Boucaud et al. (2000, 2001); Dudal et al. (2003); Ruiz Arriola et al. (2004); Megias et al. (2006); Dudal et al. (2010); Kondo (2001) that may be present in the quark propagators used to amputate the Green’s functions. Contributions of condensates to correlation functions are scale-dependent and should be evaluated using OPE as in, e.g., Ref. Blossier et al. (2011). Such analysis has not been performed yet, and the correction in Eq. (77) should be regarded as a phenomenological assumption. Another potential source of $\propto 1/p^{2}$ effects are nonperturbative infrared contributions due to potential low-momentum subdiagrams, which may appear due to the same arguments as in Ref. Aoki et al. (2008).

We perform uncorrelated fit (74) with five parameters ($Z_{I}^{\text{SI}}$, and $A,B_{1,2},C$) to the lattice data $Z_{I}^{\text{lat}}(p)$ for varying sets of momenta $p$, and two examples are shown in Fig. 12. To keep discretization errors omitted from Eq. (75) as small as possible, we include only momenta $p$ that interpolate between the 3d- and 4d-diagonals,

[TABLE]

The lowest rotational symmetry-breaking contribution $\propto p^{[4]}/p^{4}$ to Eq. (75) is shown in Fig. 11. Values $Z_{I}^{\text{lat}}(p)$ at $H(4)$-equivalent momenta $p$ are averaged. The fit range $p_{\text{min}}^{2}\leq p^{2}\leq p_{\text{max}}^{2}$ is varied with $p_{\text{min}}=1.6,\,2.0\,\text{GeV}$ and $p_{\text{max}}=3.5,\,4.0,\,4.5\,\text{GeV}$, resulting in $27\leq n_{\text{mom}}\leq 61$ lattice momentum data points that are distinct with respect to $H(4)$ transformations. We use uncorrelated $\chi^{2}$ values to evaluate goodness-of-fit and estimate systematic uncertainties from variation of the results with the fit range and the order of the perturbation theory. Although correlated fits would be preferred, we resort to uncorrelated fits, because with a small number of independent configurations $N_{\text{cfg}}=30$, it is difficult to ensure that covariance matrices of sizes $n_{\text{mom}}\sim N_{\text{cfg}}$ are estimated with uniform reliability.

The results of the fits for all fit windows using $Z^{RI,pert}$ from 1- and 2-loop perturbative calculations $Z^{RI,pert}$ are collected in Tab. 4, together with the resulting uncorrelated $\chi^{2}$ values. In order to obtain the final value, we average the central values over all the fitting methods as described in Sec. B.3. In the last row of Table 4, we show the final conversion coefficients between the lattice bare and $\overline{\text{MS}}$-renormalized operators $Q_{I}$ that take into account the difference between $N_{f}=3$ and $N_{f}=4$ QCD perturbative running (see Sec. V.2 and Eq. (73)).

Our lattice vertex functions (61) are computed with nonzero quark masses but matched to massless perturbation theory. Since we analyze only one ensemble, we cannot take the chiral limit $m_{u/d,s}\to 0$ and our renormalization can potentially have systematic bias due to quark mass dependence. While the light quark masses are small and are unlikely to have significant effect, the strange quark mass is larger and it may bias our results. Although we cannot directly assess this quark mass dependence with data at only one combination of quark masses, we can make a rough estimate of its magnitude from the quark mass dependence of the axial-vector renormalization constant $Z_{A}$. Since the operator renormalization constants in our analysis are multiplied by factors $\propto Z_{q}^{3}\propto Z_{A}^{3}$, we can estimate their corresponding quark-mass correction as

[TABLE]

where the correction $\delta Z_{A}$ due to the $m_{u/d,s}\to 0$ limit may be conservatively estimated as

[TABLE]

Using the data from Ref. Aoki et al. (2011) obtained with a very similar fermion action and lattice spacing, we find from Eq. (79) that $\delta Z_{I}/Z_{I}\approx 1.7\%$. Considering that this (likely overestimated) correction is small compared to the uncertainties quoted in Tab. 4, we neglect it in the present analysis.

VI Results

The four indepedent non-vanishing $n\overline{n}$ matrix elements in the isospin limit are given in terms of the above bare matrix elements and renormalization factors as

[TABLE]

Combining the uncertainties from ${\mathcal{M}}_{I}^{\text{lat}}$ and $Z_{I}^{SI}$ in quadrature gives the result

[TABLE]

where the first uncertainty is the combined statistical uncertainty in ${\mathcal{M}}_{I}^{\text{lat}}$ and $Z_{I}$ and the second uncertainty is the combined systematic uncertainty associated with variation in fit window described in Sec. IV, V and Appendix B.3. Quark mass effects lead to negligible systematic uncertainties because of the nearly physical pion mass used Blum et al. (2016)101010 The renormalization constants require taking the limit $m_{u/d,s}\to 0$ for matching to their exact perturbative counterparts. The associated uncertainty is estimated in Sec. V to be small and is neglected in the present study.. Uncertainties in the determination of the lattice spacing in Ref. Blum et al. (2016) are negligible compared to the fitting uncertainties in Eq. (82). Finite-volume effects have been estimated in chiral perturbation theory to be $\lesssim 1\%$ effects for the volume used for this study Bijnens and Kofoed (2017). Discretization effects are expected to be the largest unquantified systematic uncertainty that are neglected in this work. Chiral symmetry leads to $O(a)$ improvement of the fermion action, and discretization effects on meson observables for these configurations have been seen to be percent-level Blum et al. (2016). Discretization effects will be studied and removed from future calculations with multiple lattice spacing.

Final results for the $n\text{-}\overline{n}$ transition matrix elements with statistical and systematic uncertainties added in quadrature and given in Tab. 5. These results can be directly compared with MIT bag model results previously used to relate experimental results to BSM couplings Rao and Shrock (1982) as shown in Tab. 5. Two different sets of MIT bag model parameters are used to calculate $n\text{-}\overline{n}$ transition matrix elements in Ref. Rao and Shrock (1982): in fit A the up and down quark masses are set to zero, while in fit B the up and down quark masses are set to 108 MeV and a different value is used for the “bag radius” parameter. MIT bag model results for both fit A and fit B are compared to LQCD results in Tab. 5. In LQCD, the electroweak-nonsinglet matrix element ${\mathcal{M}}_{5}$ is more than an order of magnitude smaller than the electroweak-singlet matrix elements. This feature is captured by the MIT bag model, although the sign of ${\mathcal{M}}_{5}$ differs between the two bag model parametrizations. LQCD results for the electroweak-singlet operator matrix elements ${\mathcal{M}}_{1}$, ${\mathcal{M}}_{2}$, and ${\mathcal{M}}_{3}$ are larger than MIT bag model results with both parametrizations by factors of 4-8. This difference between LQCD and MIT bag model results is significantly larger than the differences between MIT bag model results with different parameter values.

The effective Lagrangian for $n\text{-}\overline{n}$ oscillations given in Eq. (46) can be used to parameterize the $n\text{-}\overline{n}$ vacuum transition rate for a generic BSM theory as

[TABLE]

where $\eta=v^{2}/\Lambda_{\text{BSM}}^{2}$ is the ratio of the Higgs v.e.v. and the BSM scale squared. Both the matrix elements ${\mathcal{M}}$ and the Wilson coefficients $\widetilde{C}^{({\mathcal{P}})}$ are scheme- and scale-dependent, and these dependencies must cancel in $\tau_{n\text{-}\overline{n}}$. Below we present results with coefficients $\widetilde{C}$ defined in $\overline{\text{MS}}$ scheme. The Wilson coefficients in Eq. (83) are predicted to be non-zero in various BSM theories, see Refs. Mohapatra (2009); Babu et al. (2013); Phillips et al. (2016) for reviews and further references, and are calculable at tree-level in QCD at BSM scales $\mu=\Lambda_{\text{BSM}}$. The $n\text{-}\overline{n}$ vacuum transition rate is given in terms of the above results by

[TABLE]

To make the prefactor dimensionless, we use the “reference” normalization scale of $700\,\text{TeV}$. Estimates based on Eq. (84) put BSM theories with scales of $\Lambda_{\text{BSM}}\sim 700$ TeV and $O(1)$ matching coefficients within reach of next-generation experiments that will be able to detect baryon number violation with $\tau_{n\text{-}\overline{n}}^{-1}\geq 10^{9}$ s Milstead (2015); Frost (2017); Fomin et al. (2018); Hewes (2017). To more precisely assess the expected signatures of theories with $B$-violation at $\Lambda_{\text{BSM}}\sim 700$ TeV, the operators can be evolved to $\mu=\Lambda_{\text{BSM}}$ using the results of Refs. Caswell et al. (1983); Buchoff and Wagman (2016),

[TABLE]

Leading-order one-loop running from 2 GeV to 700 TeV modifies the dominant matrix elements $M_{1,3}$ by up to 59%. At next-to-leading-order (NLO), two-loop running modifies $M_{1,2,3}$ by $4-9\%$ and one-loop scheme matching modifies them by $2-8\%$. Neglected next-to-next-to-leading-order perturbative renormalization effects lead to unknown systematic uncertainties estimated to be at the level of $1-3\%$ (as the square of the relative NLO effects). The $n\text{-}\overline{n}$ transition rate can be expressed in terms of the matrix elements at this scale as

[TABLE]

This result can be combined with tree-level BSM matching results for $C_{I}^{\overline{\text{MS}}}(700\text{ TeV})$ to extract constraints on BSM theory parameters from experimental constraints on $n\text{-}\overline{n}$ oscillations.

VII Conclusion

We have performed the first lattice QCD calculation of the renormalized neutron-antineutron transition matrix elements needed to extract BSM physics constraints from $n\text{-}\overline{n}$ oscillation experiments. The precision of our final results including statistical and most systematic uncertainties is $15-30\%$ for the electroweak-singlet matrix elements ${\mathcal{M}}_{1}$, ${\mathcal{M}}_{2}$, and ${\mathcal{M}}_{3}$, which can be straightforwardly improved in future calculations. Several important sources of systematic uncertainty are under control for the first time, most importantly non-perturbative renormalization, chiral symmetry violations, excited state contamination, and quark mass dependence. The two sources of systematic uncertainty that are not completely controlled in this pioneering calculation are finite volume and discretization effects. To summarize our control of common systematic uncertainties in lattice calculations:

• •

The (nearly exact) physical pion mass $m_{\pi}=139.2(4)\text{MeV}$ in our calculation eliminates the need for chiral extrapolation, which would otherwise introduce systematic uncertainties associated with low-energy effective theory. In addition, the large difference of our results from the MIT bag model may have a similar origin as the strong suppression of proton decay matrix elements found in the chiral bag model Martin and Stavenga (2012), therefore using the realistic light quark masses in our calculation is arguably the most important systematic effect we have under control.

• •

The chirally symmetric Möbius domain wall fermion action used to generate these gauge field ensembles by the RBC/UKQCD collaborations Blum et al. (2016) and compute neutron-antineutron matrix elements in this work ensures that the 14 distinct $|\Delta B|=2$ operators do not mix with each other and renormalization and conversion of lattice operators to $\overline{\text{MS}}$ scheme is free from associated uncertainties. In particular, the nonperturbatively computed operator mixing matrix in RI-MOM scheme is diagonal up to $O(10^{-3})$ corrections, which are two orders of magnitude below other uncertainties and can be safely neglected. The identical action is used for valence quarks, so this is a fully unitary calculation.

• •

Excited-state effects are accounted for using correlated two-state fits with 10 different values of $t_{\text{sep}}$ and different combinations of nucleon source and sink smearing. The energy gaps are extracted from correlated fits to nucleon two-point functions. Since we have limited statistics for such a large number of $\tau,t_{\text{sep}}$ points included in correlated fits, we use “shrinkage” estimators to obtain well-conditioned covariance matrices. We obtain systematic errors by varying the fit ranges and averaging their results weighted by the quality-of-fit figure.

• •

Renormalization effects are included through NPR in an RI-MOM scheme as described in Sec. V and one-loop matching to $\overline{\text{MS}}$ using the results of Ref. Buchoff and Wagman (2016). Some discretization effects in NPR results such as rotational symmetry breaking and $(ap)^{2}$ dependence are studied and removed by fitting the lattice data with different quark momentum scales and orientations, varying scale ranges, and comparing to 1- and 2-loop perturbative QCD running. One presently uncontrolled systematic uncertainty in our renormalization procedure is the quark mass dependence. However, a rough estimate in Sec. V suggests that this uncertainty should not exceed $1.7\%$, which is below our current level of precision; this uncertainty will be studied in the future.

• •

Although we do not control finite-volume effects directly in this study on a single ensemble, we expect them to be small. First, finite-volume effects are suppressed with $e^{-m_{\pi}L}$ where $m_{\pi}L\sim 3.9$ for the volume used for this study, which is generally considered sufficiently large for nucleon structure calculations Aoki et al. (2017). Second, chiral perturbation theory calculations in Ref. Bijnens and Kofoed (2017) estimate that finite-volume effects lead to corrections below $1\%$ to ${\mathcal{M}}_{I}$ for the volume used in this study. Future lattice calculations at additional volumes could be used to test this prediction and perform an infinite-volume extrapolation.

• •

Discretization effects are the least-controlled systematic uncertainty in our current work. Lattice QCD calculations with finer lattice spacing(s) in the immediate future will be used to fully quantify and remove discretization effects that are not controlled in this calculation. However, it is reasonable to assume that discretization effects are small compared to our current combined uncertainty from other sources. First, the chirally-symmetric fermion action that we use is automatically $O(a)$-improved. Second, the meson decay constants computed on this ensemble (before finite volume and discretization corrections are applied) are within 0.6% of the physical values ($f_{\pi}=131.1(4),\,f_{K}=156.4(4)\,\text{GeV}$ Blum et al. (2016) compared to PDG values $f_{\pi}=130.4(2),\,f_{K}=156.2(7)\,\text{GeV}$ Olive et al. (2014)). Finally, the nucleon effective mass and energy dependence on the momentum is in close agreement with the continuum limit Syritsyn et al. (2019).

Our renormalized lattice QCD results for $n\text{-}\overline{n}$ transition matrix elements provide a significant step forward in accuracy and reliability compared to previous results from quark models and preliminary lattice studies. The matrix elements predicted by QCD are found to be 4-8 times larger than the predictions of the MIT bag model for the dominant electroweak-singlet operators. This difference between our lattice results and previously available bag model results is much larger than the statistical or systematic uncertainties present in this calculation and is also much larger than the expected size of finite-volume effects that have not yet been studied directly. There is less certainty about the size of discretization artifacts; however, the automatic $O(a)$ improvement due to the chiral symmetry as well as minuscule discretization corrections in the meson decay constants, nucleon mass and dispersion relation make large discretization effects in the $n\text{-}\overline{n}$ matrix elements very unlikely.

The difference in ${\mathcal{M}}_{I}$ between the bag model and our lattice results leads to increased experimental sensitivity to baryon-number violating interactions that may cause $n\text{-}\overline{n}$ oscillations. Numbers of events that can be observed both in quasi-free neutron oscillation experiments and underground nuclear decay experiments are proportional to $\tau_{n\text{-}\overline{n}}^{-2}\propto|{\mathcal{M}}_{I}|^{2}$, therefore the $\times(4\ldots 8)$ larger values of the $n\text{-}\overline{n}$ matrix elements found in our work lead to $\times(16\ldots 64)$ increase in the event rates. Since our results are obtained from ab initio QCD calculations in a model-independent way, they must be used for more precise assessments of the potential of planned $n\text{-}\overline{n}$ oscillation searches as well as stronger constraints on theories of baryon-number violation and baryogenesis in the future.

Acknowledgements.

The authors would like to express gratitude to Yuri Kamyshkov, Rabi Mohapatra, Martin Savage, Steve Sharpe, Robert Shrock, Mike Snow, Brian Tiburzi for multiple illuminating discussions. We are indebted to Norman Christ, Bob Mawhinney, Taku Izubuchi, Oliver Witzel, and the rest of the RBC/UKQCD collaboration for access to the physical point domain-wall lattice gauge configurations used in this work. This work has been supported by the U. S. Department of Energy under grant contract no. DE-FG02-00ER41132 (INT). This work was performed also under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory (Lawrence Livermore National Security, LLC) under contract DE-AC52-07NA27344; and Brookhaven National Laboratory supported by the U. S. Department of Energy under contract DE-SC0012704. ER is supported by the RIKEN Special Postdoctoral Researcher fellowship. SS is supported by the RHIC Physics Fellow Program of the RIKEN BNL Research Center. MLW was supported by a MIT Pappalardo Fellowship and acknowledges support by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under grant Contract Number DE-SC0011090. Quark propagators and contractions were computed using the USQCD computing resources at Fermilab funded by the Office of Science of the US Department of Energy, as well as computing resources at the Lawrence Livermore National Laboratory made available through the Institutional Computing Grand Challenge program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357 The calculations were performed using Chroma Edwards and Joo (2005) and Qlua Pochinsky (senta) software packages, in particular with the efficient MDWF inverter developed by A. Pochinsky Pochinsky (sentb).

Appendix A ${\mathcal{C}}$-, ${\mathcal{P}}$-, ${\mathcal{T}}$-symmetries and nucleon states

Fermion field transformations under ${\mathcal{C}}$, ${\mathcal{P}}$, and ${\mathcal{T}}$ are given by

[TABLE]

where $C$ is given in Eq. (4) and the spin matrix $T$ is

[TABLE]

and has the property

[TABLE]

Both the color-symmetric and antisymmetric quark bilinears transform as

[TABLE]

from which the transformation properties for the 6-quark operators (35-37) follow.

In order to find the effect of these symmetries on the nucleon states and matrix elements of the operators, we spell out explicitly the neutron interpolating operators,

[TABLE]

which transform as

[TABLE]

which are used to construct (anti)neutron states on a lattice. This construction is more natural in the standard (Dirac-Pauli) basis in which the $\gamma_{4}$ matrix is diagonal. It is related to the de Grand–Rossi basis commonly used in lattice calculations by the transformation

[TABLE]

The operators (95) create the neutron and antineutron states with definite $\hat{z}$-spin as

[TABLE]

which can be found to transform as

[TABLE]

These states are used to determine the properties of the $n\text{-}\overline{n}$ matrix elements in Sec.II.4 and define them in terms of three-point functions in Sec. III.

Appendix B Statistical analysis

This work uses techniques such as bootstrap resampling that are common to lattice calculations as well as tools that are less common and detailed below: shrinkage estimation of covariance matrices, VarPro $\chi^{2}$-minimization, and weighted averaging of multiple fits with different numbers of degrees of freedom.

B.1 Shrinkage estimation of covariance matrices

Correlated $\chi^{2}$-fits require sample covariance matrices that are difficult to estimate when the number of data samples $N$ is limited compared to the number of data points $K$, as in our case. In order to estimate covariance matrices that can be safely inverted, we use the “optimal shrinkage estimator” described in Refs. Ledoit and Wolf (2004). Shrinkage involves replacing the covariance matrix with a linear combination of a well-conditioned “shrinkage target” and the original covariance matrix. It has been shown that expectation values of “shrunk” covariance matrices are closer to the true covariance matrix than the sample covariance matrix Stein (1956). The condition number of the covariance matrix is also improved by shrinkage and estimates of $\chi^{2}$ relying on the inverse covariance matrix are more robust. Shrinkage targets that better approximate the true covariance matrix naturally lead to better estimates of the true covariance matrix from a finite sample, but any prescription for defining the “shrinkage parameter” introduced below that leads to zero shrinkage in the infinite statistics limit will provide a consistent estimator for the true covariance matrix.

The estimator in Ref. Ledoit and Wolf (2004) uses a shrinkage target proportional to the $K\times K$ identity matrix $I$ where $K$ is the number of data points. However, correlation functions in lattice calculations vary over orders of magnitude if a wide range of $t_{\text{sep}}$ are used for fitting. To transform the covariance matrix into a form where the shrinkage target of Ref. Ledoit and Wolf (2004) more closely resembles the true covariance matrix, we normalize the data by subtracting the mean and diving by the square root of the variance. For data points $x_{\alpha}^{i}$ where $i=1,\cdots,N$ labels decorrelated statistical samples and $\alpha=1,\cdots,K$ labels data points (i.e. $t$ in two-point function fits and $\tau$, $t_{\text{sep}}$ in three-point function fits), define normalized data points $y_{\alpha}^{i}$ and a normalized sample correlation matrix $\rho_{\alpha\beta}$ as

[TABLE]

where the sample mean and covariance are defined as

[TABLE]

The correlation matrix with optimal shrinkage is given by

[TABLE]

where $\mu=\frac{1}{K}\mathrm{Tr}[\rho]=1$ is the mean of the spectrum of $\rho$ and the optimal shrinkage parameter $\lambda^{*}$ is defined to minimize the expected Frobenius norm $||X||=\sqrt{\mathrm{Tr}[X\,X^{T}]}$ of the difference E\big{\{}\underset{\lambda}{\mathrm{min}}||\rho(\lambda)-\varrho||^{2}\big{\}} between the estimator $\rho^{*}$ and the true correlation matrix $\varrho$. A sample estimator for the optimal shrinkage parameter is given in Ref. Ledoit and Wolf (2004)

[TABLE]

The quantity $d^{2}$ estimates the dispersion of the eigenvalues of the sample correlation matrix $\rho$, which typically has a wider spectrum and correspondingly larger (worse) condition number compared to the true correlation matrix $\varrho$. The optimal estimator (102) “shrinks” the spectrum by emphasizing the diagonal elements and makes the matrix $\Sigma^{*}$ better-conditioned, resulting in more statistically stable $\chi^{2}$ values in correlated fits. Multiplying both sides of Eq. (102) by the normalization factor in Eq. (100) yields the corresponding estimator for the covariance matrix

[TABLE]

This shrinkage prescription is therefore equivalent to an interpolation between a fully correlated fit with $\lambda^{*}=0$ (no shrinkage), and an uncorrelated fit with $\lambda^{*}=1$ (full shrinkage). Although this prescription does not provide the strictly optimal $\lambda^{*}$ minimizing the distance between $\Sigma^{*}$ and $\Sigma$, it gives a simple practical prescription for a stable and consistent choice of the shrinkage parameter. Optimal closeness between $\rho^{*}$ and $\varrho$ suggests that $\Sigma^{*}$ should provide an acceptable approximation of $\Sigma$ that is better-conditioned than $S$.

B.2 VarPro $\chi^{2}$-minimization

Fluctuations of the sample mean $G(t)=\frac{1}{N}\sum_{i}G_{i}(t)$ of an ensemble of $i=1,\dots,N$ correlation functions become Gaussian distributed as $N\rightarrow\infty$ by central limit theorems. The $\chi^{2}$ value associated with the log-likelihood of the mean correlation function becomes

[TABLE]

where $f_{n}(E_{n},t)=e^{-E_{n}t}$ for a two-point correlation function and more complicated correlation functions differ only in that $f_{n}$ has more parameters. The term quadratic in $Z_{n}Z_{m}$ can be turned into a sum of squares by transforming to the eigenbasis of

[TABLE]

Since $V_{mn}$ is a symmetric positive-definite matrix, it can be diagonalized an orthogonal transformation, which is equivalent to a change of variables in the likelihood function with a trivial Jacobian, making the eigenvalues of $V_{mn}$ new independent fit variables. The minimum of $\chi^{2}$ is determined by vanishing derivatives with respect to these eigenvalues. This provides a system of constraints that can be solved to determine the overlap factors

[TABLE]

This solution can now be substituted back into Eq. (107) to give the VarPro $\chi^{2}$ function of the energies $E$ only,

[TABLE]

The $E_{n}$ and $Z_{n}$ minimizing $\chi^{2}$ in Eq. (107) can be obtaining by determining the $E_{n}$ that minimize $\chi_{VP}^{2}$ in Eq. (110) and then solving for $Z_{n}$ from Eq. (109). More details and general discussion can be found in Refs. Golub and Pereyra (2003); Dianne P. O’Leary (2013).

B.3 Averaging over fits

Our analysis of nonperturbative renormalization and the ground-state matrix elements involves fits over different ranges of data points. The $\chi^{2}$ values of these fits cannot be directly compared due to different numbers of degrees of freedom $N_{\text{dof}}$; instead, the quality of each fit $Q$ can be assessed with its $p$-value,

[TABLE]

In order to compare and average values from a family of fits as well as estimate their stochastic and systematic uncertainties in a “blind” fashion, we use the $p$-value above as a proxy for the likelihood that these fits describe data. Thus, for any parameter $x$ extracted from a particular fit $m$ as $x_{m}\pm\delta x_{m}^{\text{stat}}$, we use the combination of its statistical uncertainty (estimated with bootstrap or jackknife) and the $p$-value of the fit $p_{m}$ as the weight

[TABLE]

to compute the “global” average value $\hat{x}$ and its statistical fluctuation $\delta\hat{x}^{\text{stat}}$,

[TABLE]

while the weighted deviation from the total average serves as the estimate of the systematic uncertainty $\delta x^{\text{sys}}$,

[TABLE]

The rationale for using the weight (112) is that it penalizes both bad fits (small $p_{m}$) and unconstraining fits (large statistical uncertainty $\delta x^{\text{stat}}$ typical of overfitting). The $(\delta x^{\text{stat}})^{-2}$ factor is motivated by similarity to a weighted average of independent data. However, since all these fits are performed on the same data set, the “global” stochastic uncertainty is computed as the inverse-squared average (instead of the sum) of individual fit uncertainties. The resemblance is especially evident if the “likelihoods” of all the fits are the same, $p_{m}=\text{const}$. Thus, the estimator (114) is also somewhat robust with respect to having similar or (nearly-)duplicated fits in the set.