Confirming the Existence of the strong CP Problem in Lattice QCD with the Gradient Flow

TL;DR

This study uses lattice QCD with gradient flow to calculate the nucleon electric dipole moment induced by the QCD theta term, providing insights into the strong CP problem and setting bounds on the theta parameter.

Contribution

First lattice QCD calculation of nucleon EDM with gradient flow controlling topological charge, extrapolated to physical pion mass and continuum limit, confirming the strong CP problem.

Findings

Nucleon EDM proportional to theta with small coefficients

Experimental bounds constrain theta to less than 2×10^{-10}

Mild discretization effects suggest continuum-like behavior

Abstract

We calculate the electric dipole moment of the nucleon induced by the QCD theta term. We use the gradient flow to define the topological charge and use $N_f = 2+1$ flavors of dynamical quarks corresponding to pion masses of $700$, $570$, and $410$ MeV, and perform an extrapolation to the physical point based on chiral perturbation theory. We perform calculations at $3$ different lattice spacings in the range of $0.07~{\rm fm} < a < 0.11$ fm at a single value of the pion mass, to enable control on discretization effects. We also investigate finite size effects using $2$ different volumes. A novel technique is applied to improve the signal-to-noise ratio in the form factor calculations. The very mild discretization effects observed suggest a continuum-like behavior of the nucleon EDM towards the chiral limit. Under this assumption our results read $d_{n}=-0.00152(71)\ \bar\theta\ e~\text{fm}$ and $d_{p}=0.0011(10)\ \bar\theta\ e~\text{fm}$. Assuming the theta term is the only source of CP violation, the experimental bound on the neutron electric dipole moment limits $\left|\bar\theta\right| < 1.98\times 10^{-10}$ ($90\%$ CL). A first attempt at calculating the nucleon Schiff moment in the continuum resulted in $S_{p} = 0.50(59)\times 10^{-4}\ \bar\theta\ e~\text{fm}^3$ and $S_{n} = -0.10(43)\times 10^{-4}\ \bar\theta\ e~\text{fm}^3$.

Tables

Table 1. Table 1 : Summary of the lattice bare parameters for the ensembles used. N G subscript 𝑁 𝐺 N_{G} is the number of gauge configurations and N corr subscript 𝑁 corr N_{\rm corr} is the number of correlation functions calculated using many stochoastically located sources for the same gauge configuration.

	$\beta$	${\kappa }_l$	${\kappa }_s$	L/a	T/a	csw	$N_G$	$N_{\mathrm{corr}}$
M1	1.90	0.13700	0.1364	32	64	1.715	322	30094
M2	1.90	0.13727	0.1364	32	64	1.715	400	20000
M3	1.90	0.13754	0.1364	32	64	1.715	444	17834
A1	1.83	0.13825	0.1371	16	32	1.761	800	15220
A2	1.90	0.13700	0.1364	20	40	1.715	789	15407
A3	2.05	0.13560	0.1351	28	56	1.628	650	12867

Table 2. Table 3 : Fit ranges [ t m ​ i ​ n , t m ​ a ​ x ] superscript 𝑡 𝑚 𝑖 𝑛 superscript 𝑡 𝑚 𝑎 𝑥 [t^{min},t^{max}] over euclidean source-sink separation t 𝑡 t used to extract the nucleon mixing angle α N subscript 𝛼 𝑁 \alpha_{N} , along with the resulting value.

ensemble	M3	M2	M1	A1	A2	A3
fit range	[10,20]	[10,20]	[10,20]	[5,11]	[7,17]	[14,21]
fitr [fm]	[0.9,1.8]	[0.9,1.8]	[0.9,1.8]	[0.6,1.3]	[0.7,1.7]	[0.96,1.43]
${\alpha }_N$	-0.040(21)	-0.190(27)	-0.142(24)	-0.099(11)	-0.103(10)	-0.105(11)

Table 3. Table 4 : The selected starting value t s m ​ i ​ n superscript subscript 𝑡 𝑠 𝑚 𝑖 𝑛 t_{s}^{min} for the fit ranges [ t s m ​ i ​ n , T / 2 ] superscript subscript 𝑡 𝑠 𝑚 𝑖 𝑛 𝑇 2 [t_{s}^{min},T/2] for the summed Q ¯ ¯ 𝑄 \mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu , the Euclidean source-sink separation fit range [ t m ​ i ​ n , t m ​ a ​ x ] superscript 𝑡 𝑚 𝑖 𝑛 superscript 𝑡 𝑚 𝑎 𝑥 [t^{min},t^{max}] and the resulting nucleon mixing angle α N subscript 𝛼 𝑁 \alpha_{N} from the selected parameters. A comparison between computing α N subscript 𝛼 𝑁 \alpha_{N} from t s m ​ i ​ n superscript subscript 𝑡 𝑠 𝑚 𝑖 𝑛 t_{s}^{min} and fitting from t s m ​ i ​ n superscript subscript 𝑡 𝑠 𝑚 𝑖 𝑛 t_{s}^{min} onwards showed a negligible difference on all ensembles. We add the values determined in a standard way for comparison.

ensemble	M3	M2	M1
$t_s^{min}/a$	7	10	6
$t_s^{min}$ [fm]	0.63	0.9	0.54
$t/a$ fit range	[10,20]	[10,20]	[10,20]
$t$ [fm] fit range	[0.9,1.8]	[0.9,1.8]	[0.9,1.8]
${\alpha }_N$ improved	-0.098(13)	-0.201(17)	-0.0822(97)
${\alpha }_N$ standard	-0.040(21)	-0.190(27)	-0.142(24)
ensemble	A3	A2	A1
$t_s^{min}/a$	10	10	10
$t_s^{min}$ [fm]	1.21	0.98	0.69
$t/a$ fit range	[5,11]	[7,17]	[14,21]
$t$ [fm] fit range	[0.61,1.34]	[0.69,1.67]	[0.96,1.44]
${\alpha }_N$ improved	-0.1016(89)	-0.1012(75)	-0.1212(67)
${\alpha }_N$ standard	-0.099(11)	-0.103(10)	-0.105(11)

Table 4. Table 5 : F 3 p / n ​ ( Q 2 → 0 ) 2 ​ M N = d p / n subscript superscript 𝐹 𝑝 𝑛 3 → superscript 𝑄 2 0 2 subscript 𝑀 𝑁 subscript 𝑑 𝑝 𝑛 \frac{F^{p/n}_{3}(Q^{2}\rightarrow 0)}{2M_{N}}=d_{p/n} fit results over M-ensembles, taken from fig.  14 .

ensemble	$m_{\pi }=410$ MeV	$m_{\pi }=570$ MeV	$m_{\pi }=700$ MeV
$d_p$ [$e$ fm]	0.0043(99)	0.0017(83)	0.0016(64)
$d_n$ [$e$ fm]	-0.0035(66)	-0.0060(53)	-0.0009(47)

Table 5. Table 6 : F 3 p / n ​ ( Q 2 → 0 ) 2 ​ M N = d p / n subscript superscript 𝐹 𝑝 𝑛 3 → superscript 𝑄 2 0 2 subscript 𝑀 𝑁 subscript 𝑑 𝑝 𝑛 \frac{F^{p/n}_{3}(Q^{2}\rightarrow 0)}{2M_{N}}=d_{p/n} fit results over A-ensembles, taken from fig.  14 .

ensemble	$a=0.1095$ fm	$a=0.0936$ fm	$a=0.0684$ fm
$d_p$ [$e$ fm]	0.0060(30)	0.0026(25)	0.0008(18)
$d_n$ [$e$ fm]	-0.0043(20)	-0.0063(20)	-0.0023(13)

Table 6. Table 7 : Fit ranges [ t s m ​ i ​ n , T 2 ] superscript subscript 𝑡 𝑠 𝑚 𝑖 𝑛 𝑇 2 [t_{s}^{min},\frac{T}{2}] , over the symmetrically summed Q ¯ ¯ 𝑄 \mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu time t s subscript 𝑡 𝑠 t_{s} and resulting in the improved EDM determination F 3 n ​ ( Q 2 → 0 ) 2 ​ M N ≡ d ¯ n subscript superscript 𝐹 𝑛 3 → superscript 𝑄 2 0 2 subscript 𝑀 𝑁 subscript ¯ 𝑑 𝑛 \frac{F^{n}_{3}(Q^{2}\rightarrow 0)}{2M_{N}}\equiv\mkern 1.5mu\overline{\mkern-1.5mud\mkern-1.5mu}\mkern 1.5mu_{n} , over the M-ensembles, taken from fig.  21 . The unimproved results d n subscript 𝑑 𝑛 d_{n} from tab.  5 are included for comparison. The values determined at t s m ​ i ​ n superscript subscript 𝑡 𝑠 𝑚 𝑖 𝑛 t_{s}^{min} differ by the fit results at most by 10 % percent 10 10\% of the error associated.

ensemble	$m_{\pi }=410$ MeV	$m_{\pi }=570$ MeV	$m_{\pi }=700$ MeV
fit range	[6,32]	[7,32]	[4,32]
fitr [fm]	[0.54,2.9]	[0.63,2.9]	[0.63,2.9]
$d_n$ [$e$ fm]	-0.0045(26)	-0.0090(27)	-0.0027(20)
${\overline{d}}_n$ [$e$ fm]	-0.0035(66)	-0.0060(53)	-0.0009(47)

Table 7. Table 8 : Fit ranges [ t s m ​ i ​ n , T 2 ] superscript subscript 𝑡 𝑠 𝑚 𝑖 𝑛 𝑇 2 [t_{s}^{min},\frac{T}{2}] , over the symmetrically summed Q ¯ ¯ 𝑄 \mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu time t s subscript 𝑡 𝑠 t_{s} and resulting in the improved EDM determination F 3 n ​ ( Q 2 → 0 ) 2 ​ M N ≡ d ¯ n subscript superscript 𝐹 𝑛 3 → superscript 𝑄 2 0 2 subscript 𝑀 𝑁 subscript ¯ 𝑑 𝑛 \frac{F^{n}_{3}(Q^{2}\rightarrow 0)}{2M_{N}}\equiv\mkern 1.5mu\overline{\mkern-1.5mud\mkern-1.5mu}\mkern 1.5mu_{n} , over the A-ensembles, taken from fig.  21 . The unimproved results d n subscript 𝑑 𝑛 d_{n} from tab.  6 are included for comparison. The values determined at t s m ​ i ​ n superscript subscript 𝑡 𝑠 𝑚 𝑖 𝑛 t_{s}^{min} differ by the fit results at most by 10 % percent 10 10\% of the error associated.

ensemble	$a=0.1095$ fm	$a=0.0936$ fm	$a=0.0684$ fm
fit range	[3,16]	[4,20]	[10,28]
fitr [fm]	[0.36,1.9]	[0.39,2.0]	[0.69,1.9]
$d_n$ [$e$ fm]	-0.0048(13)	-.00393(97)	-0.0044(10)
${\overline{d}}_n$ [$e$ fm]	-0.0043(20)	-0.0063(20)	-0.0023(13)

Table 8. Table 9 : Neutron and proton EDM fit parameters C 1 subscript 𝐶 1 C_{1} , C 2 subscript 𝐶 2 C_{2} , C 3 subscript 𝐶 3 C_{3} extracted from the combine fits to all 6 ensembles using eq. ( 41 ), as well as the resulting χ P ​ D ​ F 2 subscript superscript 𝜒 2 𝑃 𝐷 𝐹 \chi^{2}_{PDF} . We also estimate g ¯ 0 subscript ¯ 𝑔 0 \bar{g}_{0} using C 2 subscript 𝐶 2 C_{2} and eq. ( 2 ).

	$C_1$ $\left[\overline{\theta }e{\text{ fm}}^3\right]$	$C_2$ $\left[\overline{\theta }e{\text{ fm}}^3\right]$	$C_3$ $\left[\frac{\overline{\theta }e\text{ fm}}{{\text{fm}}^2}\right]$	${\chi }_{PDF}^2$	${\overline{g}}_0^{\overline{\theta }}$ $\left[\overline{\theta }\right]$
proton	$-3.6\left(5.3\right)\times {10}^{-4}$	$-6.8\left(6.6\right)\times {10}^{-4}$	$0.20\left(31\right)$	$2.0\left(1.4\right)$	$-9.9\left(9.6\right)\times {10}^{-3}$
neutron	$3.1\left(3.2\right)\times {10}^{-4}$	$8.8\left(4.4\right)\times {10}^{-4}$	$-0.16\left(23\right)$	$1.8\left(1.5\right)$	$-12.8\left(6.4\right)\times {10}^{-3}$

Table 9. Table 10 : Neutron and proton Schiff fit parameters C 4 subscript 𝐶 4 C_{4} , C 5 subscript 𝐶 5 C_{5} extracted from the combine fits to all 6 ensembles using eq. ( 55 ), as well as the resulting χ P ​ D ​ F 2 subscript superscript 𝜒 2 𝑃 𝐷 𝐹 \chi^{2}_{PDF} . We additionally include the ratio C 2 C 4 subscript 𝐶 2 subscript 𝐶 4 \frac{C_{2}}{C_{4}} , where C 2 subscript 𝐶 2 C_{2} is the second fit parameter result from tab.  9 .

	$C_4$ $\left[\overline{\theta }e{\text{ fm}}^3\right]$	$C_5$ $\left[\frac{\overline{\theta }e{\text{ fm}}^3}{{\text{fm}}^2}\right]$	${\chi }_{PDF}^2$	$\frac{C_2}{C_4}\approx -6$
proton	$0.50\left(59\right)\times {10}^{-4}$	$-0.0022\left(73\right)$	$1.25\left(80\right)$	$-20\left(200\right)$
neutron	$-0.10\left(43\right)\times {10}^{-4}$	$-0.0057\left(51\right)$	$1.37\left(97\right)$	$70\left(970\right)$

Equations

$${\mathcal{L}}_{\mathrm{QCD}}=-\frac{1}{4}{G}_{\mu \nu }^a{G}^{a,\mu \nu }+\bar{q}(iDto0.0pt\text{hss}/-M)q-\bar{\theta }\frac{g^2}{64{\pi }^2}{ϵ}^{\mu \nu \alpha \beta }{G}_{\mu \nu }^a{G}_{\alpha \beta }^a,$$

$${\mathcal{L}}_{QCD}=-\frac{1}{4}{G}_{\mu \nu }^a{G}^{a,\mu \nu }+\bar{q}iDto0.0pt\text{hss}/q-\bar{m}\bar{q}q-\epsilon \bar{m}\bar{q}{\tau }_{3}q+m_{\ast }\bar{\theta }\bar{q}i{\gamma }^{5}q,$$

$$d_n(\theta )$$

$$d_p(\theta )$$

$${\mathcal{L}}_{\pi N}(\theta )={\bar{g}}_0\bar{N}\vec{\pi }\cdot \vec{\tau }N$$

$${\bar{g}}_0=\frac{(m_n-m_p{)}^{\mathrm{strong}}(1-{\epsilon }^2)}{4F_{\pi }\epsilon }\bar{\theta }=-14.7(2.3)\times 10^{-3}\bar{\theta },$$

$$d_n(\theta )$$

$$d_p(\theta )$$

$${⟨O⟩}_{\bar{\theta }}=⟨O⟩+i\bar{\theta }⟨OQ⟩+\mathcal{O}({\bar{\theta }}^2),$$

$$Q(t_f)=a^4\sum _{x}q(x,t_f),$$

$$q(x,t_f)=\frac{1}{64{\pi }^2}{ϵ}_{\mu \nu \rho \sigma }{G}_{\mu \nu }^a(x,t_f)G_{\rho \sigma }^a(x,t_f),$$

$$\chi (t_f)=\frac{a^8}{V}\sum _{x,y}⟨q(x,t_f)q(y,t_f)⟩$$

$$G_2({\mathbf{p}}^{\prime },t,\Pi )=a^3\sum _{\mathbf{x}}{e}^{-i{\mathbf{p}}^{\prime }\cdot \mathbf{x}}\mathrm{Tr}\left\{\Pi ⟨\mathcal{N}(\mathbf{x},t)\overline{\mathcal{N}}(0,0)⟩\right\},$$

$$G_2({\mathbf{p}}^{\prime },t,\Pi )=\frac{e^{-E_{{\beta }_0}t}}{2E_{{\beta }_0}}\mathrm{Tr}\left\{\Pi ⟨\Omega |\mathcal{N}|{\beta }_0⟩⟨{\beta }_0|\overline{\mathcal{N}}|\Omega ⟩\right\},$$

$$M_{\mathrm{eff}}({\mathbf{p}}^{\prime }=0,t,{\Pi }_{+})=\log \left[\frac{G_2({\mathbf{p}}^{\prime }=0,t,{\Pi }_{+})}{G_2({\mathbf{p}}^{\prime }=0,t+1,{\Pi }_{+})}\right]=m_{{\beta }_0^{+}},$$

$$G_2^{(Q)}({\mathbf{p}}^{\prime },t,\Pi ,t_f)=a^3\sum _{\mathbf{x}}{e}^{-i{\mathbf{p}}^{\prime }\cdot \mathbf{x}}\mathrm{Tr}\left\{\Pi ⟨\mathcal{N}(\mathbf{x},t)\overline{\mathcal{N}}(0,0)Q(t_f)⟩\right\}.$$

$${\alpha }_N=\frac{G_2^{(Q)}({\mathbf{p}}^{\prime }=0,t,{\gamma }_5{\Pi }_{+},t_f)}{G_2({\mathbf{p}}^{\prime }=0,t,{\Pi }_{+})},$$

$$\overline{Q}({\tau }_Q,t_f)=a^3\sum _{\mathbf{x}}q(\mathbf{x},{\tau }_Q;t_f),Q(t_f)=a\sum _{{\tau }_Q}\overline{Q}({\tau }_Q,t_f).$$

$${\Delta }_2^{(\overline{Q})}(t,\Pi ,t_f,{\tau }_Q)=a^3\sum _{\mathbf{x}}\mathrm{Tr}\left\{\Pi ⟨\mathcal{N}(\mathbf{x},t)\overline{Q}({\tau }_Q,t_f)\overline{\mathcal{N}}(0,0)⟩\right\}$$

$$G_2^{(Q)}(t,\Pi ,t_f)=a\sum _{\frac{{\tau }_Q}{a}=0}^{T/a}{\Delta }_2^{(\overline{Q})}(t,\Pi ,t_f,{\tau }_Q).$$

$$\overline{G}{}_2^{(\overline{Q})}(t,\Pi ,t_f,t_s)$$

$$G_2^{(Q)}(t,\Pi ,t_f)=\overline{G}{}_2^{(\overline{Q})}(t,\Pi ,t_f,t_s=T/2).$$

$$\frac{F_3^{p/n}(Q^2)}{2M_N}\stackrel{Q^2\ll m_{\pi }^2}{\to }{d}_{p/n}-S_{p/n}{Q}^2+O(Q^4),$$

$$G_3^{(Q)}({\mathbf{p}}^{\prime },t,\mathbf{q},\tau ,\Pi ,{\gamma }_{\mu },t_f)$$

$$G_3({\mathbf{p}}^{\prime },t,\mathbf{q},\tau ,\Pi ,{\gamma }_{\mu })$$

$${\mathcal{J}}_{\mu }(\mathbf{y},\tau )=\frac{4}{3}\overline{u}(\mathbf{y},\tau ){\gamma }_{\mu }u(\mathbf{y},\tau )-\frac{1}{3}\overline{d}(\mathbf{y},\tau ){\gamma }_{\mu }d(\mathbf{y},\tau ),$$

$$R({\mathbf{p}}^{\prime },t,\mathbf{q},\tau ,\Pi ,{\gamma }_{\mu })$$

$$R^{(Q)}({\mathbf{p}}^{\prime },t,\mathbf{q},\tau ,\Pi ,{\gamma }_{\mu },t_f)$$

$$K({\mathbf{p}}^{\prime },t,\mathbf{p},\tau )\equiv \sqrt{\frac{G_2({\mathbf{p}}^{\prime },\tau ,{\Pi }_{+})G_2({\mathbf{p}}^{\prime },t,{\Pi }_{+})G_2(\mathbf{p},t-\tau ,{\Pi }_{+})}{G_2(\mathbf{p},\tau ,{\Pi }_{+})G_2(\mathbf{p},t,{\Pi }_{+})G_2({\mathbf{p}}^{\prime },t-\tau ,{\Pi }_{+})}}.$$

$$R({\mathbf{p}}^{\prime },t,\mathbf{q},\tau ,\Pi ,{\gamma }_{\mu })=A(E_{{\mathbf{p}}^{\prime }},E_{\mathbf{p}})\text{Tr}\left\{\Pi (-i/p^{\prime }+m){\Gamma }_{\mu }(Q^2)(-i/p+m)\right\},$$

$${\Gamma }_{\mu }(Q^2)={\gamma }_{\mu }{F}_1(Q^2)+\frac{{\sigma }_{\mu \nu }{q}_{\nu }}{2m}{F}_2(Q^2).$$

$$A(E_{{\mathbf{p}}^{\prime }},E_{\mathbf{p}})=\frac{1}{4\sqrt{E_{{\mathbf{p}}^{\prime }}{E}_{\mathbf{p}}(E_{{\mathbf{p}}^{\prime }}+m)(E_{\mathbf{p}}+m)}}.$$

Full text

**Confirming the Existence of the strong CP Problem

in Lattice QCD with the Gradient Flow **

**Jack Dragosa, Thomas Luub, Andrea Shindlera,

Jordy de Vriesc,d, Ahmed Yousif${}^{\,a}$**

a Facility for Rare Isotope Beams, Physics Department, Michigan State University, East Lansing, Michigan.

b Institute for Advanced Simulation (IAS-4), Institut für Kernphysik (IKP-3) and JARA-HPC, FZJ, Jülich, Germany

c *Amherst Center for Fundamental Interactions, Department of Physics, University of Massachusetts Amherst, Amherst, MA 01003, USA *

d RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, New York 11973-5000, USA

We calculate the electric dipole moment of the nucleon induced by the QCD theta term. We use the gradient flow to define the topological charge and use $N_{f}=2+1$ flavors of dynamical quarks corresponding to pion masses of $700$, $570$, and $410$ MeV, and perform an extrapolation to the physical point based on chiral perturbation theory. We perform calculations at $3$ different lattice spacings in the range of $0.07~{}{\rm fm}<a<0.11$ fm at a single value of the pion mass, to enable control on discretization effects. We also investigate finite size effects using $2$ different volumes. A novel technique is applied to improve the signal-to-noise ratio in the form factor calculations. The very mild discretization effects observed suggest a continuum-like behavior of the nucleon EDM towards the chiral limit. Under this assumption our results read ${d_{n}=-0.00152(71)\ \bar{\theta}\ e\text{~{}fm}}$ and ${d_{p}=0.0011(10)\ \bar{\theta}\ e\text{~{}fm}}$. Assuming the theta term is the only source of CP violation, the experimental bound on the neutron electric dipole moment limits $\left|\bar{\theta}\right|<1.98\times 10^{-10}$ ($90\%$ CL). A first attempt at calculating the nucleon Schiff moment in the continuum resulted in ${S_{p}=0.50(59)\times 10^{-4}\ \bar{\theta}\ e}$ fm3 and ${S_{n}=-0.10(43)\times 10^{-4}\ \bar{\theta}\ e}$ fm3.

1 Introduction

A nonzero measurement of the electric dipole moment (EDM) of the nucleon in the foreseeable future would be a clear signal of new physics, since the known CP-violating phase of the CKM matrix leads to EDMs that lie orders of magnitude below current experimental limits. The source of a nonzero EDM could then either be the QCD $\bar{\theta}$-term or higher-dimension CP-violating quark-gluon operators that originate in beyond-the-Standard Model (BSM) physics, or a combination of these two. To interpret an EDM signal or lack thereof, and to possibly disentangle the source (e.g. $\bar{\theta}$-term or BSM), requires a non-perturbative calculation linking the CP-violating sources to the hadronic observables.

Lattice QCD can calculate the nucleon EDM directly in terms of CP-violating operators at the quark level. Various attempts have been made in this regard [1, 2, 3, 4, 5, 6, 7]. However, the renormalization of CP-violating operators within a lattice (discretised) formulation of QCD is very non-trivial, and for several operators presents large difficulties in interpreting lattice results. Further, the $\bar{\theta}$ term itself introduces a complex phase in the determinant of the quark matrix, which produces a sign problem and precludes the use of standard stochastic methods. Several techniques have been used to address the $\bar{\theta}$-term contribution to the EDM and attempts have been made to solve the complicated renormalization patterns of the CP-violating operators [8]. We refer to the recent review [9] for a summary.

We proposed to use the gradient flow to calculate all CP-violating source to the EDM in refs. [10, 11], and presented preliminary results [12, 13, 14, 15]. In this paper we consider the $\bar{\theta}$-term contribution to the EDM in a perturbative manner as discussed in [11]. This is well justified considering the stringent constraints on $\bar{\theta}$ set by EDM experiments. In this way we avoid the problem of a complex fermionic determinant. To define the QCD $\bar{\theta}$ term we use the gradient flow. The topological charge defined in this way has a finite and well defined continuum limit [16, 17, 18]. It is much faster to compute than using the Ginsparg-Wilson definition and it is theoretically more robust than definitions using cooling techniques. Another problem that hinders lattice calculations of the nucleon EDM is the very poor signal-to-noise ratio. In this respect we explore a novel technique to determine the space-time region where the signal in the relevant correlation functions is maximized. A first account of this technique has already been presented [12].

Additional insights into the EDM of the nucleon [19, 20, 21, 22, 23] and nuclei [24, 25] can be provided by chiral effective field theory. The CP-violating quark operators are translated to effective CP-violating hadronic operators and EDMs depend on the unknown low-energy constants (LECs) of the theory. The LECs can be estimated from dimensional analysis or, preferably, be determined from experiments and/or by lattice-QCD calculations. We use these insights from chiral calculations to understand the pion mass dependence of our results, and to connect our nucleon EDM calculations to nuclear EDMs.

The remainder of the paper is organized as follows: section 2 gives a cursory discussion of the phenomenology of the nucleon EDM, followed by an overview of the lattice details and parameters in section 3, where we discuss the general lattice strategy used and we define the basic observables, including the gradient flow. Sanity checks are given in section 4, where we compute and display the topological charge using the gradient flow. The nucleon two-point correlation function is explored in section 5, then the setup of the computation of the EDM is derived and results shown in section 6. A discussion follows the results section in section 7, where we discuss the ramifications of our results and compare our results to the literature. Finally, we conclude in section 8.

2 Phenomenology of the QCD theta term.

The discrete space-time symmetries parity $(P)$ and time-reversal $(T)$, and hence via the $CPT$ theorem also $CP$ symmetry, are broken in QCD by the $\bar{\theta}$ term. In the case of two quark flavors the QCD Lagrangian in Minkowski space is given by

[TABLE]

where $q=(u\,,d)^{T}$ denotes the quark doublet containing up and down quarks, $G^{a}_{\mu\nu}$ is the gluon field strength tensor, $\epsilon^{\mu\nu\alpha\beta}$ ($\epsilon^{0123}=+1$) is the completely antisymmetric tensor, $D_{\mu}$ the gauge-covariant derivative, $M$ the real $2\times 2$ quark-mass matrix, and $\bar{\theta}$ the coupling of the $CP$-odd interaction. In eq. (1) the complex phase of the quark-mass matrix has been absorbed in the physical parameter $\bar{\theta}=\theta+\mathrm{arg}\,\mathrm{det}(M)$. For the application of chiral perturbation theory ($\chi$PT) it is useful to perform an anomalous axial $U(1)$ transformation to replace the $CP$-odd gluonic term in favor of a complex mass term [19, 26]. Under the assumption that $\bar{\theta}\ll 1$ the QCD Lagrangian can then be written as

[TABLE]

where we have defined the average quark mass $\bar{m}=(m_{u}+m_{d})/2$, the quark-mass difference $\varepsilon=(m_{u}-m_{d})/(m_{u}+m_{d})$, and the reduced quark mass $m_{\ast}=m_{u}m_{d}/(m_{u}+m_{d})=\bar{m}(1-\varepsilon^{2})/2$.

The QCD $\bar{\theta}$ term induces an EDM in hadrons and nuclei. The first EDM search occured with the neutron in 1957 [27, 28]. Till this day, no signal has been found for the EDM, despite measurement sensivities having improved by six orders of magnitude. The current bound ${d_{n}<3.0\times 10^{-26}}$ e cm [29, 30] sets strong limits on the size of $\bar{\theta}$ and sources of $CP$ violation from physics beyond the SM [31].

In order to set a bound on the $\bar{\theta}$ term, it is necessary to calculate the dependence of the neutron EDM on $\bar{\theta}$ [19]. One way to do this is by using $\chi$PT. In the first step one derives interactions between the low-energy degrees of freedom, pion and nucleons (and heavier hadrons), that violate $CP$ and transform the same way under chiral symmetry as the complex mass term in eq. (2). In the second step, one combines the chiral $CP$-odd interactions with the standard $CP$-even chiral Lagrangian to calculate the nucleon EDM. This calculation has been done up to next-to-leading order (NLO) in both $SU(2)$ and $SU(3)$ $\chi$PT [20, 22] and gives in the two-flavored theory for the neutron ($d_{n}$) and proton ($d_{p}$) EDM:

[TABLE]

in terms of $g_{A}\simeq 1.27$ the strong pion-nucleon coupling constant, $F_{\pi}\simeq 92.4$ MeV the pion decay constant, $m_{\pi}$ and $m_{N}$ the pion and nucleon mass respectively, $e>0$ the proton charge, and three low-energy constants (LECs) of $CP$-odd chiral interactions $\bar{g}_{0}$ and $\bar{d}_{p/n}$. The first term in brackets in eq. (2) arises from the leading-order (LO) one-loop diagram involving the $CP$-odd vertex

[TABLE]

in terms of the nucleon doublet $N=(p\,n)^{T}$ and the pion triplet $\vec{\pi}$. The LO loop is divergent and the divergence and associated scale dependence have been absorbed into the counter terms $\bar{d}_{p/n}$ which signify contributions to the nucleon EDMs from short-range dynamics and appear at the same order as the LO loop diagrams. The second term in brackets in eq. (2) arises from finite next-to-leading-order (NLO) diagrams.

Because the $\bar{\theta}$ term breaks chiral symmetry as a complex quark mass, the LEC $\bar{g}_{0}$ can be related to known $CP$-even LECs using chiral symmetry arguments [19, 32, 33]

[TABLE]

where $(m_{n}-m_{p})^{\mathrm{strong}}$ is the quark-mass induced part of the proton-neutron mass splitting for which we used the recent lattice results [34, 35]. Inserting eq. (7.2) in eq. (2) we obtain

[TABLE]

Under the assumption that the terms $\bar{d}_{p/n}$ do not cancel against the calculable loop contributions, a comparison with the experimental bound gives the strong constraint $\bar{\theta}\leq 10^{-10}$. Clearly, a more reliable constraint on $\bar{\theta}$ requires a direct nonperturbative calculation of the full nucleon EDMs. This is the main goal of this work.

In the isoscalar combination $d_{n}+d_{p}$ the loop contribution cancels out to a large extent. For observables sensitive to this combination, such as the deuteron EDM whose measurement is the goal of the JEDI collaboration [36], a first-principle calculation of the total nucleon EDM is important. EDMs of light nuclei have been calculated as a function of $\bar{\theta}$ in ref. [25]. Nuclear EDMs get contributions from the single-nucleon EDMs and from the CP-violating nucleon-nucleon potential which is dominated by one-pion-exchange terms. The latter depend mainly on $\bar{g}_{0}$ and are therefore relatively well under control. The dominant remaining uncertainty is the size of the nucleon EDMs as a function of $\bar{\theta}$. With nonperturbative calculations of nucleon EDMs induced by the $\bar{\theta}$ term, we immediately obtain predictions for EDMs of light nuclei. With future improvements of nuclear theory even EDMs of diamagnetic atoms such as 199Hg and 225Ra could be directly given as a function of $\bar{\theta}$.

3 Lattice QCD action and numerical details

We discretize the QCD action on an hypercubic lattice with spacing $a$ and volume $L^{3}\times T$. The fermionic part of our QCD lattice action is the non-perturbatively O($a$)-improved Wilson action with $N_{f}=2+1$ dynamical quarks. The gauge part is the Iwasaki gauge action. For our calculation we have always used valence quarks with the same lattice action and the same bare parameters as the sea quark action, that is to say our framework is fully unitary.

We performed calculations using the publicly available PACS-CS gauge fields available through the ILDG [37]. We used $6$ different ensembles that allow us to study discretization effects, finite-size effects and pion mass dependence. We studied the pion-mass dependence with $3$ ensembles at $3$ different bare quark masses, at $L/a=32$ and $T/a=64$ and a lattice spacing $a=0.0907(13)$ fm. The lattice spacing and the physical point were determined with the experimental input of $m_{\pi},m_{K}$, and $m_{\Omega}$. More details on these ensembles are available in ref. [38] and are summarized in the first three M rows of tabs. 2, 2.

To study discretization effects we used $3$ ensembles with $3$ different lattice spacings but with the same volume, $L\simeq 1.8$ fm. The ratios of masses in the pseudoscalar and vector channels differ, between the $3$ ensembles, at most by $1\%$ in the light sector and at most of $3\%$ in the strange quark sector. These very small mismatches are irrelevant for all purposes for our scaling violation study. The lattice spacings and quark masses in these ensembles are also determined using $m_{\pi}$, $m_{K}$, and $m_{\phi}$. Details for these ensembles can be found in ref. [39] and summarized in the last three A rows of tabs. 2, 2.

Among the $6$ ensembles described above there are $2$ ensembles, M1 and A2, with the same bare parameters, $\beta=1.9$, $\kappa_{l}=0.13700$, $\kappa_{s}=0.1364$ and different lattice volumes with $L/a=20$, $L/a=32$ and $T=2L$. These $2$ ensembles allow us to investigate finite-size effects.

To improve the overlap with the ground state of the relevant matrix elements in the two- and three-point functions, we applied a Gaussian gauge-invariant smearing [40, 41] at the source and at the sink of our quark propagators. Using the notation of refs. [40, 41] we use $64$ iterations of the smearing algorithm with a smearing fraction of $\alpha=0.39$ using the definition in ref. [40]. These parameters corresponds to a spatial root-mean-square radius for the nucleon interpolating operator of around $0.4$ fm. The quality of our projection into the ground state can be evaluated from figs. 5(a), 5(b).

For the vector form factors studied here, we used the renormalization factor $Z_{V}$ determined using vector Ward identities in ref. [42] and summarized in tabs. 2, 2.

The strategy we use in this paper is a perturbative expansion in powers of $\bar{\theta}$ (the expansion is fully performed in Euclidean space). This is justified by the small value of $\bar{\theta}$ estimated from experimental constraints. With this strategy every correlation function $\braket{O}_{\bar{\theta}}$, evaluated in a $\bar{\theta}$ vacuum, is determined from a small-$\bar{\theta}$ expansion

[TABLE]

where $O$ is some multi-local operator, $\bar{\theta}$ is the coefficient for the CP-violating term, and $Q$ is the topological charge. The expectation values on the r.h.s of eq. (7) are computed on a standard QCD background. This allows us to use lattice QCD gauge configurations without generating new gauges for this specific calculation. We define the correlation functions used to determine the nucleon EDM induced by the $\bar{\theta}$ term in the next section.

The topological charge CP-violating operator that enters the correlation functions must in principle be normalized. We use the gradient flow [16] to define the topological charge which in this way has a finite continuum limit and does not need any additional normalization [16, 17, 18]. The reason is that the flowed fields are free from ultraviolet divergences [16, 17] for all positive flow times, $t_{f}>0$. Additionally it can be shown that the topological charge defined with the gradient flow is flow-time independent [18] for all positive flow times, $t_{f}>0$, in the continuum limit. Details on how we numerically perform the flowing of the fields have been described in ref. [11].

4 Topological charge and the gradient flow

We define the topological charge at finite lattice spacing as

[TABLE]

where the topological charge density reads

[TABLE]

and $G_{\mu\nu}^{a}(x,t_{f})$ is a lattice discretization of the continuum field tensor defined with flowed gauge fields. As a lattice definition for the field tensor, we use the discretization suggested in ref. [43]. This definition suffers from small discretization effects and, in fact, the corresponding topological susceptibility

[TABLE]

is flow-time independent starting from a flow-time radius, $\sqrt{8t_{f}}$, of about $1$ fm for all lattice spacings we have investigated. This can be seen in fig. 1, where we show the flow-time dependence of the topological susceptibility computed for all M- (left) and A-ensembles (right). As expected, the region where the susceptibility is independent of the flow time extends towards smaller flow-time values for smaller lattice spacings.

We used the topological charge to perform various checks on the quality of the ensembles. An important check for EDM calculations is to make sure that the ensembles sample the field space in such a way that no spurious CP-violation is induced. In other words we must check the expectation value $\langle Q(t_{f})\rangle=0$ within statistical errors.

In fig. 3 we show $\langle Q(t_{f})\rangle$ evaluated on all our ensembles for various pion masses and lattice spacings. To properly estimate the statistical uncertainties we evaluate the autocorrelation function and the corresponding integrated autocorrelation time, $\tau_{\rm int}$, as defined in ref. [44]. For all our ensembles the average topological charge vanishes within statistical errors. In fig. 3 we show the flow-time dependence of the integrated autocorrelation time for the topological charge. As expected the gradient flow, by smoothing out some of the short-distance fluctuations, allows a better determination of $\tau_{\rm int}$ that reaches a plateau for $\sqrt{8t_{f}}\simeq 0.2$ fm for all our ensembles [45, 46].

The integrated autocorrelation time $\tau_{\rm int}$ we obtain falls within the range $7<\tau_{\rm int}<35$ for the M1 and M2 ensembles and slightly smaller, $3<\tau_{\rm int}<10$, for our M3 ensemble. We attribute this behavior with the rather short Markov Chain for the M3 ensemble which most likely does not allow a more accurate determination of its $\tau_{\rm int}$. We also observe from fig. 3 that $\tau_{\rm int}$ increases as we decrease the lattice spacing. This is an expected result [47, 45] since the tunnelling between different topological sectors becomes increasingly difficult with decreasing lattice spacing, meaning that the sampling of different sectors, which would decrease $\tau_{\rm int}$, is lessened.

For completeness in fig. 4 we show the difference of error determination if we were to use a standard resampling technique, such as bootstrap, instead of the error determination using the autocorrelation function. In this case, $\langle Q(t_{f})\rangle\neq 0$ within uncertainties. This demonstrates how a robust uncertainty determination for the topological charge requires both the estimate of the autocorrelation function and its corresponding integrated autocorrelation time.

5 Two-point correlation functions and the nucleon mixing angle

In this section we analyze the nucleon two-point correlation function, to extract the effective mass, as well as the nucleon mixing angle. The standard two-point correlation function with sink momentum ${\bm{p}}^{\,\prime}$ has the form

[TABLE]

where $\Pi$ is some spin projector, and $\mathcal{N}$ is an interpolating field with the quantum numbers of a nucleon, inserted with a source-sink time separation of $t$. The spectral decomposition for this equation in the limit where $T\gg t\gg 0$, keeping implicit a sum over the polarizations, is

[TABLE]

where the lowest energy state $\beta_{0}$ for which $\bra{\Omega}\mathcal{N}\ket{\beta_{0}}\neq 0$ and $\bra{\beta_{0}}\mkern 1.5mu\overline{\mkern-1.5mu\mathcal{N}\mkern-1.5mu}\mkern 1.5mu\ket{\Omega}\neq 0$ arises from the approximation111It is clear from the context when we consider operators as in eq. (12) or interpolating fields as in eq. (11). $T\gg t\gg 0$.

The effective mass, which is shown in fig. 5, is given by the simple log ratio

[TABLE]

where $\Pi_{+}=(I+\gamma_{4})/2$ is the positive parity projector, $\beta^{+}_{0}$ is the lowest energy positive parity nucleon state, and again, $T\gg t\gg 0$. In fig. 5(a), we compare our effective mass determinations for the M-ensembles to those computed in ref. [38] and find agreement within statistical errors. As shown in fig. 5(b), we observe lattice-spacing dependence of the order of $10\%$ between the finest and coarsest lattices.

The nucleon mixing angle [1], $\alpha_{N}$, can be extracted by utilizing the two-point correlator from eq. (11), and the $\bar{\theta}$-modified two-point correlator

[TABLE]

The mixing angle is defined using the small-$\bar{\theta}$ expansion as

[TABLE]

in the region where $t\gg 0$ and $t_{f}\gg 0$.

In figs. 6(a) and 8(a), we show the dependence of the nuclear mixing angle on the source-sink separation $t$ (in fm) for the M- and A-ensembles, respectively. For the M-ensembles in figs. 6(a), there is little to no excited-state contamination effects for ${t>0.6}$ fm. The results suggest a non-trivial chiral behavior for $\alpha_{N}$. We will discuss this in detail in sec. 6 when we describe our EDM determination. For the lattice-spacing ensembles in fig. 8(a), we require a minimum source-sink separation of $t\approx a\,\{5,7,14\}$ fm (for $a=\{0.1095,0.0936,0.0684\}$ fm ensembles) to plateau and achieve ground state saturation. The final plateaued quantity for the different A-ensembles all lead to the same angle $\alpha_{N}$ and all results are consistent within statistical uncertainties. Similarly, fig. 8(b) demonstrates that computing $\alpha_{N}$ at $m_{\pi}=700$ MeV with different box sizes lead to consistent results. In tab. 3 we summarize the fit ranges and resulting values for $\alpha_{N}$.

The value of the flow-time radius, $\sqrt{8t_{f}}$, for all these analyses is fixed around $0.5-0.6$ fm where the topological charge is least affected by lattice artifacts. Figs. 6(c), 8(c) show the nucleon mixing angle plotted against the flow time $\sqrt{8t_{f}}$ at a fixed source-sink separation $t$ for the M- and A-ensembles. For the M-ensembles, we see no flow-time dependence after $\sqrt{8t_{f}}>0.2-0.3$ fm, confirming that the results obtained in this region are free from gradient-flow discretization effects. A similar conclusion is reached for the A-ensembles.

In figs. 6(b) we show the integrated autocorrelation time of $\alpha_{N}$ for the M-ensemble results shown in fig. 6(a). For the M1 and M2 ensembles a factor of $2-4$ increase in autocorrelation as the source-sink separation approaches 0. Fortunately, a minimum source-sink separation of ${t\approx 1}$ fm greatly decreases the autocorrelation correction that we apply in the determination of the nucleon mixing angle $\alpha_{N}$. Most importantly, in comparison to $\braket{Q}$ from fig. 3 (i.e. not in the presence of a nucleon), the autocorrelation effect is dramatically decreased by a factor of at least $\simeq 4$. We attribute this effect to the presence of a fermionic part, $\mathcal{N}\mkern 1.5mu\overline{\mkern-1.5mu\mathcal{N}\mkern-1.5mu}\mkern 1.5mu$, in the correlation function. Numerical evidence suggests that the observables considered in this work containing fermion lines, such as $\alpha_{N}$ and the EDM are less coupled to the slow modes contributing to the spectral decomposition of the autocorrelation function [45]. As this effect will be greater when analyzing the EDM (from three-point correlation functions), we resort to our standard bootstrap error propagation technique for the final EDM computation. We checked explicitly that error estimates from a bootstrap and an autocorrelation analysis give consistent results.

5.1 Improving the Nucleon Mixing Angle

In this section we describe a method previously explored in [48], that aims to reduce the statistical uncertainty of the determination of the nucleon mixing angle $\alpha_{N}$. The strategy can be described as an attempt to understand the space-time region where the overlap between the topological charge density and the fermionic part of the correlation function is maximal. To perform this investigation we define a spatially-summed topological charge density

[TABLE]

We then numerically study the dependence on $\tau_{Q}$ of $\alpha_{N}$ and corroborate our numerical findings with a spectral decomposition of the relevant correlators.

The ratio $\alpha_{N}$ and the modified two-point correlator $G_{2}^{(Q)}$ have the same $\tau_{Q}$ dependence and we therefore focus on the latter. Setting $\bm{p}^{\,\prime}=\bm{0}$ and omitting it in our expressions, we define

[TABLE]

where the correlator in eq. (14) can be obtained by summing $\tau_{Q}$ from [math] to the time extent of the lattice $T$

[TABLE]

To focus on the region where the signal resides, we sum the spatially-summed topological charge density, $\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu(\tau_{Q},t_{f})$, symmetrically starting from the source location. That is we sum $\tau_{Q}$ starting from [math] (and $T$) up to a value $t_{s}$ (and $T-t_{s}$). The goal is to find a summation window $t_{s}$ small enough such that we capture all the signal and avoid the summation of unnecessary “noise”. We define the partial summed correlator

[TABLE]

from which, using the periodicity of our lattice, the original correlator in eq. (14) is obtained as

[TABLE]

Although there are other choices for the starting point of our summation in $\tau_{Q}$, we only consider starting from $\tau_{Q}=0$. In app. A we derive a spectral decomposition for the correlator in eq. (17). We argue that in the limit $t_{s}\gg t\gg 0$, the partially-summed correlator $\mkern 1.5mu\overline{\mkern-1.5muG\mkern-1.5mu}\mkern 1.5mu^{(\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu)}_{2}(t,\Pi,t_{f},t_{s})$ is independent of $t_{s}$ and $t$, up to exponentially suppressed corrections. These corrections seem to be rather small, and in fact our numerical experiments indicate that we can safely stop the summation over $\tau_{Q}$ at $t_{s}\simeq t$. In this way we avoid to sum in the region between $t$ and $T/2$ where numerically the correlators seem to vanish up to statistical fluctuations.

We first fix the source-sink separation $t$ to a large enough value such that effects from excited states are suppressed. We then study the dependence of $\alpha_{N}$ on the summation window $t_{s}$. In fig. 10 we show the $t_{s}$ dependence of $\alpha_{N}$, for the M-ensembles (left), and A-ensembles (middle) and two different physical volumes corresponding to M1 and A2 ensembles (right). In all ensembles we observe that $\alpha_{N}$ reaches a plateau when $t_{s}\simeq t$, consistent with the expectation that contributions for $t_{s}>t$ are exponentially suppressed and below our statistical accuracy. We do observe a very small drift of $\alpha_{N}$ for larger values of $t_{s}$ for the ensembles M1, and a smaller drift for the ensemble M3, for small values of $t$. We attribute this to statistical fluctuations that could arise from small local parity-violating effects induced by non-vanishing matrix elements of $\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu(\tau_{Q},t_{f})$ between two states of the same parity, $\bra{\beta}\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu(\tau_{Q},t_{f})\ket{\beta}\neq 0$. These local fluctuations are averaged out when the charge density is summed over the whole space-time volume as shown in fig. 3. Nevertheless, all values of $\alpha_{N}$ determined with the improved method are statistically compatible with the results obtained with the standard analysis.

To compare the improved extraction of $\alpha_{N}$ to the standard determination described in sec. 5, we show in figs. 12, 12 the standard and improved determination of $\alpha_{N}$ as a function of the Euclidean source-sink separation $t$. The values of $t_{s}$ considered are summarized in tab. 4. We observe a signal-to-noise improvement in all our ensembles, up to a factor $2$, with the most significant observed in the ensembles M1, M2 and A3. We observe the largest discrepancy between the improved and unimproved methods, of the order of $2.3~{}\sigma$, in the M1 and M3 ensembles. We attribute this discrepancy to standard statistical fluctuations of the gauge fields. A summary of the fit ranges and results for the improved nucleon mixing angle is given in tab. 4, where for companions we added the values of $\alpha_{N}$ determined in the standard way.

6 Electric Dipole Moment Results

The neutron (n) and proton (p) EDMs, $d_{p/n}$, can be extracted from the CP-odd electric dipole form factor222The general form not requiring $Q^{2}\ll m_{\pi}^{2}$, is given in eq. (53) and discussed in detail in secs. 6.2 and 6.3. We performed the same analysis with the fit function in eq. (53) and found insignificant changes to the EDM results (see secs. 6.2 and 6.3).

[TABLE]

which requires a lattice QCD computation of $F_{3}(Q^{2})$. The variable $Q^{2}$ in this case refers to the momentum transfer and should not be confused with the topological charge. The small $\bar{\theta}$ expansion provides us a way of accessing $F_{3}$ from three-point correlation functions without the need for generating new gauge configurations at finite $\bar{\theta}$ and without relying on a problematic analytical continuation to imaginary $\bar{\theta}$. To access $F_{3}(Q^{2})$, we calculate the following three-point correlation functions with and without the insertion of the topological charge, respectively,

[TABLE]

where the electromagnetic current in terms of the quark currents is given by

[TABLE]

and $\mathcal{N}$ denotes standard proton or neutron interpolating fields.

Once the three-point correlation functions are computed, we remove the leading Euclidean time dependence and nucleon-to-vacuum amplitude contributions via the ratios:

[TABLE]

where we have implicitly defined ratios for the proton and the neutron. We define the square-root factor as

[TABLE]

The spectral decomposition of the ratio function $R$ in eq. (6), in the limit $T\gg t\gg 0$, reads

[TABLE]

where the vector form factor contains all terms allowed by the symmetries of the theory

[TABLE]

For completeness the expression of $A(E_{\bm{p}^{\,\prime}},E_{\bm{p}})$ reads

[TABLE]

The data that we computed coming from the fixed-sink method is such that $\bm{p}^{\,\prime}=\bm{0}$ (which implies $\bm{q}=-\bm{p}$), simplifying eq. (26) to

[TABLE]

The analogous modified ratio function with the insertion of the topological charge $R^{(Q)}$ in eq. (6) has, retaining only the ground state contribution, the following spectral decomposition at leading order in $\bar{\theta}$

[TABLE]

where

[TABLE]

Due to subtleties between lattice quantities and physical, the form factor decomposition in presence of a CP-violating operator insertion, is written in terms of modified form factors, $\widetilde{F}_{2}(Q^{2})$ and $\widetilde{F}_{3}(Q^{2})$, related to the physical form factors by [6]

[TABLE]

The rotated form factor $F_{3}(Q^{2})$ corresponds to the actual electric dipole form factor as measured in experiments. From now on, we will focus on this quantity.

The ratio functions $R$ and $R^{(Q)}$ become constant, as long as the large-time approximation $T\gg t\gg\tau\gg 0$ is satisfied to ensure ground-state dominance. As the fixed-sink method is employed to compute the three-point correlation functions, a region in which this large time approximation is satisfied for $\tau$ can be found and we denote the results of the fits as

[TABLE]

The technique for fitting these ratio functions over $\tau$ is described in app. B. With this construction, a system of equations can be solved for form factors $F_{i}(Q^{2})$ , $i=1,2,3$ of the form:

[TABLE]

where the collective index $A$ denotes any combination of the indices $A=\{j,k,l\}$. In other words, we run over all possible combinations of projectors $\Pi$, all current momentum $\bm{q}$ within a given $Q^{2}$, and operator gamma matrix $\gamma_{\mu}$. The index $A$ of the matrix $\mathcal{A}_{Ai}(Q^{2})$ corresponds to the coefficients for each form factor $F_{i}$ for the corresponding ratio function $R$ or $R^{(Q)}$. These coefficients are found by analyzing the spectral decomposition of $R$ or $R^{(Q)}$, which needs to be done for every evaluated index $A$.

Using eq. (21), we extrapolate to $F^{p/n}_{3}(Q^{2}\rightarrow 0)/(2M_{N})=d_{p/n}$. We use a linear plus constant fit function, giving the extrapolated value $d_{p/n}$ at $Q^{2}\rightarrow 0$ (as well as slope in $Q^{2}$ providing $S_{p/n}$).

The final extraction of the neutron (left) and proton (right) CP-odd form factor $\frac{F_{3}(Q^{2})}{2M_{N}}$ is shown for the M-ensembles in fig. 14 and for the A-ensembles in fig. 14. Fig. 14 shows that all M-ensembles are statistically consistent evaluated, and with zero. Fig. 14 shows that there are no major discretization effects, as all the extrapolated $Q^{2}\rightarrow 0$ results are consistent.

The following figs. 15, 17, 17 are all displayed to understand the systematic effects resulting from varying the flow time $t_{f}$, and different methods of determining the nucleon mixing angle $\alpha_{N}$ used in the form factor decomposition ($\mathcal{A}_{Aj}(Q^{2})$ in eq. (35)). In fig. 15, for example, we show how the form factors $F_{3}$, determined at different flow-time radii ${\sqrt{8t_{f}}=0.60,0.65,0.70}$ fm (green, red and blue), are statistically consistent for all three M-ensembles (left to right). From both fig. 17, where the improved method (see sec. 5.1) of determining the nucleon mixing angle $\alpha_{N}$ (in red) is compared to the standard method for $\alpha_{N}$ (in blue), and fig. 17, where we vary the fit range for extracting $\alpha_{N}$, it is clear that a more precise determination of $\alpha_{N}$ has a negligible impact on improving the precision of the results for the CP-odd form factor $F_{3}$. A summary of the $Q^{2}\rightarrow 0$ extrapolations for different ensembles is given in tabs. 5, 6.

6.1 Improving the Modified Three-Point Correlation Function

In this section, we utilize a similar improvement technique used for $\alpha_{N}$, but now applied to the modified three-point correlation function $G^{(Q)}_{3}$. The improvement starts by analyzing the time dependence of the spatially integrated topological charge density

[TABLE]

where $\tau_{Q}$ signifies the temporal location of the topological charge $\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu$ defined in eq. (16). The spectral decomposition for this correlator has the form:

[TABLE]

where $\alpha,\beta,\gamma$ and $\delta$ are labels for the states propagating, and the 0 subscript indicates the lowest energy state propagating with the appropriate quantum numbers. We stress that $t_{f}\neq 0$ implies the absence of any contact terms. From fig. 18, a clear signal is observed at $\tau_{Q}=0$ on all ensembles. This motivates summing $\tau_{Q}$ symmetrically around $\tau_{Q}=0$ to obtain the summed three-point correlator:

[TABLE]

The resulting fit function to the sum range $t_{s}$, for the range $t_{s}>t$, is:

[TABLE]

Where $\gamma_{\pm}$ and $\beta_{\pm}$ represent the positive and negative ($\pm$) parity nucleon states. $A_{0}$, $A_{\gamma_{\pm}},A_{\beta_{\pm},\gamma_{\mp}}$, and $A_{\beta_{\pm}}$ are combinations of nucleon matrix elements, $E_{\gamma_{\pm}}$ is the energy of the propagating state $\gamma_{\pm}$, and $E_{0_{\pm}}$ is the lowest energy of the positive and negative parity nucleon states $0_{\pm}$. We construct the improved ratio function $\mkern 1.5mu\overline{\mkern-1.5muR\mkern-1.5mu}\mkern 1.5mu^{(\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu)}$ in the same way as in eq. (6), but with the replacement of $G^{(Q)}_{3}\rightarrow\mkern 1.5mu\overline{\mkern-1.5muG\mkern-1.5mu}\mkern 1.5mu^{(\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu)}_{3}$. The value of the correlation function $G^{(Q)}_{3}$, used to extract the CP-odd form factor, is obtained in the limit $t_{s}\rightarrow T$. If the summation over $\tau_{Q}$ is performed up to a value $t_{s}<T$ the neglected terms will be exponentially small as one can deduce from eq. (39). Our numerical results seem to indicate that indeed the neglected contributions for intermediate values of $t_{s}$ are well below the statistical accuracy of our calculation.

In fig. 19, the results for the symmetrically summed topological charge ratio function $\mkern 1.5mu\overline{\mkern-1.5muR\mkern-1.5mu}\mkern 1.5mu^{(\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu)}$ are shown as a function of the sum range $t_{s}$. In all cases, a plateau can be observed at $t_{s}=\tau$. This indicates that all the exponential terms in eq. (39) are suppressed for $t_{s}>\tau$. Coupled with the large statistical noise inherent in the data, we fit the result with a constant value once the plateau has formed. These fit ranges are displayed in tables 7 and 8, and are used for the form factor analysis in sec. 6.2.

Finally, fig. 20 displays a standard modified ratio function $R_{3}^{(Q)}$ plot over current insertion time $\tau$, where the improved ratio function (blue) is compared with the standard method (red). The improved ratio function uses the “min” time from tables 7 and 8.

In fig. 21, a comparison between the improved ratio functions (blue) and the standard integrated topological charge (red) used in the extraction of the neutron CP-violating form factor $\frac{F_{3}(Q^{2})}{2M_{N}}$ is shown. In all cases, a two-to-three times increase in the signal-to-noise is observed and all results are statistically consistent333We have at most $1.5~{}\sigma$ disagreement between the two methods at $Q^{2}\rightarrow 0$ for the $a=0.0684$ fm ensemble..

6.2 Continuum Extrapolated Results with Improved Ratio Functions

Armed with the improved results for the nucleon EDMs, the next step entails the extrapolation to the physical pion mass and the continuum limit. From $\chi$PT we learn that the leading dependence of the nucleon EDMs on the pion mass is given by [49]

[TABLE]

where $C_{1}$ and $C_{2}$ are fit constants. To account for the finite lattice spacing, we include an additional fit parameter, $C_{3}$,

[TABLE]

The additional term ensures that the EDM only vanishes in the chiral limit after taking the continuum limit. We have performed a global fit with eq. (41) taking into account our $6$ data points from ensembles A1-A3 and M1-M3. In the four plots in figs. 22, 23, we show the EDM results for the proton and neutron separately as function of the pion mass and lattice spacing.

Specifically, in Fig. 22 we show the extraction of the neutron (left) and proton (right) EDM plotted against their $m_{\pi}^{2}$ values (in MeV). The blue band shows the extrapolation using the fit function in eq. (41), evaluated at $d_{p/n}(a=0,m_{\pi})$. This function evaluated at the physical pion mass is what we are interested in. In red we show the same extrapolation, where the fit is evaluated instead at $d_{p/n}(a=0.09\text{ fm},m_{\pi})$, to study the role of discretization errors. In particular, we observe an uncertainty of the EDMs at the physical pion mass that is roughly twice larger at $a=0.09$ fm. It is perhaps surprising that the uncertainty at the physical point reduces in the continuum limit. But the reason is clear. By fitting the nucleon EDMs to the fit function in eq. (41), the uncertainty on the fit parameters $C_{1}$ and $C_{2}$ is increased by the presence of the $C_{3}$ term. Now that the $a^{2}$ dependence is taken into account, we can perform an interpolation between the EDM in the chiral limit and the pion masses of our ensembles. In the continuum limit, $a=0$, the resulting nucleon EDM at the physical pion mass has now less uncertainty because $d_{n,p}(a=0,m_{\pi}=0)$ while $d_{n,p}(a>0,m_{\pi}=0)\neq 0$ and unconstrained.

The final continuum extrapolation values for the neutron and proton EDM are

[TABLE]

and we include the determination for the fit parameters $C_{1}$, $C_{2}$, $C_{3}$ of eq. (41), as well as the chi-squared per degree of freedom parameter, $\chi^{2}_{PDF}$, in tab. 9. The error on $\chi^{2}_{PDF}$ is determined from the bootstrap samples distribution. Since the correlators for the proton and the neutron EDM are different, it is possible to obtain different relative uncertainties in the two cases. It is not clear to us though, why we observe a relative larger uncertainty for the proton than for the neutron.

In fig. 23 we show the dependence of our EDM results on the lattice spacing $a$ for the neutron (left) and proton (right) EDM. Overlaid on top, we have the evaluation of the fit function eq. (41) at two different values of $m_{\pi}$: $m_{\pi}=700$ MeV (purple band) and $m_{\pi}=m_{\pi}^{phys}$ (green band). The ensembles analyized in this work do not allow us to study mass-dependent discretization effects, but we can still observe the impact of the chiral interpolation on the continuum limit. The continuum extrapolation has less uncertainty, thanks to the constraint that the EDM vanishes in the chiral limit. Adding more ensembles to study mass dependence cutoff effects is certainly desirable, but it does not change the main conclusion of this analysis.

We can extract a value of the CP-odd pion-nucleon LEC, $\bar{g}_{0}$, which plays an important role in the EDMs of nuclei and diamagnetic atoms, by identifying our result for fit parameter $C_{2}$ with the coefficient of the log term in eq. (2) for the neutron EDM. This gives the relation

[TABLE]

leading to the extraction

[TABLE]

at the physical pion mass. This result is in good agreement with the chiral perturbation theory prediction in eq. (5) and confirms the applicability of the fit function in eq. (41). A consistent result, with larger uncertainties is obtained for the proton EDM (see tab. 9).

6.3 Schiff Moment of the Proton and Neutron

Apart from the EDMs of the neutron and the proton, the nucleon electric dipole form factor (EDFF) contains additional information. The EDFF can be decomposed as

[TABLE]

where $d_{p/n}$ denotes the proton or neutron EDM, $S_{p/n}$ denotes the proton or neutron Schiff moments defined by ${S_{p/n}=(2M_{N})^{-1}(dF^{p/n}_{3}/dQ^{2})|_{Q^{2}=0}}$, and $H_{p/n}$ are functions that capture the remaining $Q^{2}$ dependence. Chiral perturbation theory allows for a calculation of the Schiff Moments and the $H_{p/n}$ functions from the analogous isovector and isoscalar quantities [49]. At leading order in the chiral expansion the nucleon EDMs are given in eq. (2). The leading-order Schiff moments are isovector and given by

[TABLE]

where we have used eq. (5). NLO corrections have been calculated in ref. [22] and reduce the leading-order result by roughly $50\%$ and provide a tiny contribution to $S_{n}+S_{p}=\mathcal{O}(10^{-5}\,\bar{\theta}\,e\,\mathrm{fm}^{3})$. $\chi$PT thus predicts that the neutron and proton Schiff moments are equal in magnitude but with opposite sign. The leading-order $H_{p/n}$ are also isovector and given by

[TABLE]

where

[TABLE]

In the limit $Q^{2}\ll m_{\pi}^{2}$ the nucleon EDFFs become

[TABLE]

such that the Schiff moments provide the dominant $Q^{2}$ dependence of the EDFFs. The nucleon EDMs and the LEC $\bar{g}_{0}$ are induced by the $\bar{\theta}$ term and scale as $d_{p/n}\sim\bar{g}_{0}\sim\bar{m}_{*}\bar{\theta}\sim m_{\pi}^{2}\bar{\theta}$. As such, the Schiff moments scale as $S_{p/n}\sim\bar{g}_{0}/m_{\pi}^{2}$ which is pion mass independent. This statement is potentially confusing as we infer from eq. (2) that the $\bar{\theta}$ term decouples in the chiral limit and the whole nucleon EDFF should vanish. eq. (50), however, requires $Q^{2}\ll m_{\pi}^{2}$. In the opposite limit, we obtain

[TABLE]

and the EDFFs vanish in the chiral limit as expected.

The goal is to extract $S_{p/n}$ from our lattice data as this allows for a direct comparison to the $\chi$PT prediction in eq. (46) and the extraction in the previous section based on the pion mass dependence of the nucleon EDMs. To extract $S_{p/n}$, we first extrapolate our results to small $Q^{2}$ by fitting the EDFF to the function

[TABLE]

The effects of the $h_{1}^{(0)}$ function turns out to have minimal impact on the extraction of $d_{p/n}(m_{\pi}^{2},a^{2})$ and $S_{p/n}(m_{\pi}^{2},a^{2})$, and we obtain similar results if we use the fit function

[TABLE]

This shows that our results are not precise enough to isolate the more subtle $Q^{2}$ behavior.

Once we have obtained $S_{p/n}(m_{\pi}^{2},a^{2})$ we can extrapolate to the continuum limit and the physical pion mass. LO $\chi$PT predicts no dependence on the pion mass, and, having an $\mathcal{O}(a)$ improved lattice action, we add a quadratic dependence on the lattice spacing $a$

[TABLE]

with $C_{4}$ and $C_{5}$ fit constants. The results for the Schiff moments along with the continuum extrapolation are shown in figs. 24, 25. In fig. 25 we show the fit results with the $a^{2}$ dependence. We observe minimal discretization effects over the range $a=\{0\rightarrow 0.12\}$ fm. In fig. 24 we show the fit results as a function of the pion mass $m_{\pi}$. At $a=0$ we perform a constant fit in the pion mass, to obtain the proton and neutron Schiff moments at the physical point. We do not extrapolate to the chiral limit because the $\chi$PT prediction that $S_{p/n}$ are pion-mass independent will break down at some point as inferred from eq. (51). We obtain for the Schiff moments at the physical point

[TABLE]

as well as the fit parameters $C_{4}=S_{p/n}$ and $C_{5}$ from performing this fit in tab. 10. The uncertainties are significant and the magnitudes are somewhat below the LO $\chi$PT predictions in eq. (46), but in better agreement once $\chi$PT NLO corrections are included. There is some evidence for a dominantly isovector Schiff moment as predicted from $\chi$PT, but the uncertainties are too large to make strong statements. We perform a sanity check of our result by comparing the ChPT predictions for the fit coefficient $C_{2}$ and $C_{4}$. From eqs. (43) and (46), we infer the LO ChPT prediction

[TABLE]

Our fit values for this ratio are given in tab. 10, and agree with this prediction within (large) statistical errors.

7 Discussion

In this section, we discuss the EDM and Schiff moment results for the neutron and proton. The succeeding Section 7.1 compares our determination of the EDM to previous lattice QCD EDM computations. Then following in Section 7.2, the phenomenological ramifications of our results for the EDM and Schiff moments are discussed.

7.1 Comparison with other works

In this section, we compare the results obtained for the neutron EDM $d_{n}$ with few lattice QCD results from the literature. As noted in [6], it is sometimes problematic to compare different EDM calculations, as most results preceding this paper do not consider the rotation of the CP-odd form factor $F_{3}$ with $F_{2}$ and $\alpha_{N}$ computed on the lattice. 444We note that no general consensus has been reached about the need to perform this rotation of the form factors. The EDM $d_{n}$ (rotated to $F_{3}$) is shown in Fig. 26. Good agreement is seen from our results to the others [1, 3] at $m_{\pi}\approx 475$ MeV, but we see a slight tension between the results of [4, 3] and our results at $m_{\pi}\approx 350$ MeV. It must be stressed that the rotation requires knowledge of the phase $\alpha_{N}$ and the unrotated form factors $F_{2}$ and $F_{3}$ which are not always easy to extract. To rotate the “C. Alexandrou et al, 2016” [4] results, an estimation of $F_{2}$ was determined from [50]. To rotate the “F.-K. Guo et al, 2015” [3] results, $F_{2}$ was determined from [51] (at $\bar{\theta}=0$) and $\alpha_{N}$ and $F_{3}$ estimated via a linear+cubic fit in $\bar{\theta}$ performed by [6].

In particular, the lattice results for $F_{3}$ not obtained in this work do not take into account correlations between $F_{2}$, $F_{3}$ and $\alpha_{N}$. As such, Fig. 26 is mainly shown for illustrative purposes and the error estimates for results not obtained in this work should be taken with a grain of salt.

7.2 Impact on EDMs of light nuclei

Armed with a non-perturbative determination of the nucleon EDMs as a function of $\bar{\theta}$ we can revisit EDMs of systems with more than a single nucleon. EDM experiments so far have mainly focused on neutral systems, but EDMs of charged particles can be probed if the particles are trapped in electromagnetic storage rings [52]. This technique has lead to a direct limit on the EDM of the muon [53], and to plans to pursue EDM measurements of protons and light nuclei in dedicated storage rings. Such measurements are still far away but impressive progress has been reported in refs. [54, 55]. EDMs of light nuclei have been calculated using chiral EFT [24, 25]

[TABLE]

in terms of the EDMs of nucleons and the CPV pion-nucleon coupling constants $\bar{g}_{0}$ and $\bar{g}_{1}$ associated to the interactions

[TABLE]

Values of $\bar{g}_{0}$ and $\bar{g}_{1}$ can be obtained from chiral-symmetry arguments [56, 33] that link these LECs to the hadron spectrum

[TABLE]

and the smallness of $\bar{g}_{1}/\bar{g}_{0}$ is due to approximate isospin symmetry.

In absence of direct lattice calculations of $d_{n}$ and $d_{p}$ we could only predict values for the combinations

[TABLE]

But with our lattice determination of $d_{n}(\bar{\theta})$ and $d_{p}(\bar{\theta})$, we can now estimate the EDMs of light ions directly in terms of $\bar{\theta}$

[TABLE]

Due the dependence on the isoscalar nucleon EDM, $d_{n}+d_{p}$, the deuteron EDM is still very uncertain. The tri-nucleon EDMs however are predicted more than three standard deviations from zero and with a fixed sign. The total uncertainty arises in roughly equal amounts from uncertainties in eq. (7.2) and in the determination of the nucleon EDMs in eq. (6.2). If nonzero EDMs are measured these relations can be used to differentiate between the SM theta term and BSM sources of CP violation. They can also provide indirect evidence for the existence of a Peccei-Quinn mechanism, by finding EDM patterns in disagreement with eq. (7.2) [57].

8 Conclusion

In this paper we computed the proton and neutron EDM from dynamical lattice QCD using various pion masses at different lattice spacings and volumes, as enumerated in tabs. 2, 2. We found our results have rather small (within our statistical uncertainties) discretization effects, which greatly simplified our continuum limit extrapolations. We found satisfactory agreement with existing results, as discussed in section 7. With our measurements at multiple pion masses, we performed a chiral interpolation to obtain, at the physical pion mass and in the continuum limit ${d_{n}=-0.00152(71)\ \bar{\theta}\ e\text{~{}fm}}$ and ${d_{p}=0.0011(10)\ \bar{\theta}\ e\text{~{}fm}}$. The nonzero result for the neutron EDM confirms the existence of the strong CP problem at two standard deviations and limits $\left|\bar{\theta}\right|<1.98\times 10^{-10}$ at $90\%$ C.L. The dependence of the nucleon EDMs on the pion mass allowed us to extract of the CP-odd pion-nucleon coupling, $\bar{g}_{0}$. The resulting value is in good agreement with chiral perturbation theory.

Important to our analysis was the implementation of the gradient flow on the topological charge. In fact we can perform the continuum limit at fixed flow time with no need to calculate the normalization of the topological charge. As we have discussed and documented in section 5 and ensuing subsections, the gradient flow also allows a more robust determination of the integrated autocorrelation time that must be taken into account when estimating statistical uncertainties.

Be that as it may, the extraction of a non-zero EDM is notoriously difficult due to its poor signal-to-noise ratio. To address this issue we have employed a novel technique, first documented in [12], for reducing the noise in our measured observables. Instead of summing all space-time points in the calculation of ratios relevant for the extraction of our 3-point function and $\alpha_{N}$ term, this method focuses on the space-time region where the signal is strongest. We argued that the neglected space-time region gives exponentially suppressed contributions to the correlation functions and this expectation has been confirmed by our numerical results. On some ensembles, this method enabled us to increase the signal to noise by a factor of $\approx 2$. This method was described in detail in subsection 6.1 and Appendix A.

We have also analysed the $Q^{2}$ dependence of our form factors and performed an extraction of the Schiff moment with dynamical fermions. Our results of $S_{p}=0.50(59)\times 10^{-4}\ \bar{\theta}\ e$ fm3 for the proton and $S_{n}=-0.10(43)\times 10^{-4}\ \bar{\theta}\ e$ fm3 for the neutron that are in reasonable agreement with chiral perturbation theory predictions. Our estimates for this value can be improved upon with more statistics and calculations on larger lattices (and thus lower $Q^{2}$ points), which would allow for a more robust extraction.

Our calculation is a big step towards a precise determination of the nucleon EDM and Schiff moment. Improvement of these results will most definitely come from increased statistics, and more calculations at different pion masses at several lattice spacings. We comment here on the necessity to perform calculations at the physical pion mass. In the chiral limit the EDM induced by the $\theta$ term vanishes (i.e. $d_{p/n}=0$ at $m_{\pi}=0$). In our view, given the small value of the EDM induced by the $\theta$ term and the additional standard reduction in signal-to-noise as the pion mass is lowered for nucleon correlators, calculations of these quantities at the physical pion mass have possibly less to gain than those at higher pion masses. Because of this constraint, it could be advantageous to have results at slightly heavier-than-physical pion masses and then robustly interpolate to the physical pion mass using $\chi$PT. The subsequent errors of the interpolation are stable and easily quantified precisely because one is doing an interpolation and not an extrapolation. We remark though that, for a chirally–breaking lattice action, such as the non-perturbatively O($a$) improved Wilson-clover fermion action we have adopted, the nucleon EDMs vanish in the chiral limit only after performing the continuum limit. This emphasize the importance of the continuum limit when using a chirally-breaking action. In this respect the gradient flow allows us to perform a safe study of discretization effects.

To summarize, the ideal scenario of a direct determination at the physical point with statistical uncertainties under control, can be circumvented by simply investing more time in lattice QCD calculations at slightly heavier pion masses (where the signal-to-noise is not as prohibitive). It goes without saying that calculations of the EDM at heavier-than-physical pion masses can potentially be more cost effective than the physical-pion-mass calculations only if one has a robust description of the lattice data with $\chi$PT.

Acknowledgments

We would like to thank Mattia Bruno, Taku Izubuchi and Sergey Syritsyn for valuable discussions. This work was supported in part by Michigan State University through computational resources, provided by the Institute for Cyber-Enabled Research, and in part through financial support from the Deutsche Forschungsgemeinschaft (Sino-German CRC 110). The authors gratefully acknowledge the computing time granted through JARA-HPC on the supercomputer JURECA [58] at Forschungszentrum Jülich.

Appendix A Alpha Improvement Derivation

Staring with the general three-point correlation function

[TABLE]

we handle the time ordering in the next two sections, by performing the spectral decomposition for $t>\tau$ and $\tau<t$. This is the general expression for an arbitrary operator $O$ and spin projector $\Pi$. For the computation of the nucleon mixing angle $\alpha_{N}$, we have $O=\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu$.

A.1 Case $t>\tau$

Starting with the specific time ordering $t>\tau$ in eq. (66), we perform the standard spectral decomposition to produce the correlation function

[TABLE]

where the sum over states $\alpha,\beta,\gamma$ have been reduced to states that only contain momenta $\mathbf{p}_{\gamma}=\mathbf{q}$, $\mathbf{p}_{\beta}=\mathbf{p}-\mathbf{q}$ and $\mathbf{p}_{\alpha}=\mathbf{p}$. The two approximation one can apply to this equation are $T\gg t$ and $t\gg 0$, which are related to the source-sink separation of the two-point correlation function

[TABLE]

where $\alpha_{0}$ is the lowest lying energy state that gives a non-zero contribution to $\Delta_{2}^{(O)}$

[TABLE]

A.2 Case $t<\tau$

This has the same form, with replacing $\mathcal{N}\leftrightarrow O$

[TABLE]

The two approximations $T\gg t$ and $t\gg 0$ are again applied

[TABLE]

where this time, $\gamma_{0}$ is the lowest lying state that gives a non-zero contribution to $\Delta_{2}^{(O)}$

[TABLE]

A.3 Total Form

Over the total range $\tau\in[0,T]$, the expression is

[TABLE]

A.4 Symmetric Partially Summed Current $t_{s}>t$

This region is where the fit will take place, so we omit the derivation for $t_{s}<t$.

[TABLE]

noting the second sum is shifted to $\tau\in[t+a,t_{s}]$ using the lattice spacing increment $a$.

One thing to note here, is the terms $\tau=0$ and $\tau=t$ are contact terms, which need to be handled properly (operator product expansion, or gradient flow). Next we group $\tau$ terms in preparation for the $t_{s}$ sum

[TABLE]

The sums can be computed by using:

[TABLE]

substituting in this expression, and including the terms where $\alpha=\beta$ and $\gamma=\beta$ separately

[TABLE]

As the final 2 terms are only exponentially dependent on $t_{s}$, we write these terms as exponentials of single energy indices

[TABLE]

and clumping like terms

[TABLE]

A.5 Explicit form for $O=\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu$

As $\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu$ is a parity violating operator, the nucleon states that propagate before and after this operator must be opposite in parity. This removes the first and second terms as $\bra{\beta}\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu\ket{\beta}=0$. As well as this, the terms with sums over two terms either require ($\alpha,\beta=\alpha_{+},\beta_{-}$) or ($\alpha,\beta=\alpha_{-},\beta_{+}$) where the subscript $\pm$ refers to the state having positive or negative parity. The projector is selected to be $\Pi=\gamma_{5}\Pi_{+}=\gamma_{5}\frac{I+\gamma_{4}}{2}$, which results in only the trace term with ($\alpha,\beta=\alpha_{+},\beta_{-}$) being non-zero

[TABLE]

The terms $e^{-E_{\alpha_{+}}(T+t)}$ and $e^{-E_{\alpha_{+}}(T-t_{s})}$ in the final sum are exponentially suppressed as ${T\gg T/2\geq t_{s}}$ and ${T\gg t}$

[TABLE]

From this complicated expression, the $t_{s}$ dependence only appears exponentially in the final term. Therefore, we can fit the two-point correlation function with

[TABLE]

Due to the statistical noise of the data and high correlation in the data with respect to $t_{s}$, we elected to neglect the excited state term by fitting a constant in the region where $Be^{-Et_{s}}\ll A$.

Appendix B Ratio Function Fit Range Selection

In this appendix, we present the technique used for extracting the CP-odd form factor $F_{3}(Q^{2})$ from the ratio function in eq. (6). Since only constant (”one-state”) fits are implemented for the ratio functions, careful consideration to excited state effects is needed.

The method employed to account for fit range dependence in our error estimates, is to include multiple fit ranges that satisfy some $\chi^{2}$ per degree of freedom ($\chi^{2}_{PDF}$) criterium. For this study, we only select fits that satisfy $\chi^{2}_{PDF}\in[0.5,1]$. Using the multiple fit range deturminations of $R$ and $R^{Q}$, we extend eq. (35) to include different fit ranges:

[TABLE]

where the extra index $f(A)$ refers to which fit range is used, which depends on the collective index $A=\{j,k,l\},\,A\in[1,\ldots,N_{A}]$.

Since each ratio function selected by index $A$ has $f(A)$ different ways to extract the quantity, the system is solved for every combination of $f(A)\,\forall\,A\in[1,\ldots,N_{A}]$. This results in $\prod_{A}F(A)$ independent system of equations to solve, where $F(A)$ is the number of different fits accepted (using the $\chi^{2}_{PDF}$ criterium) for index $f(A)=1,2,\ldots,F(A)$.

Once the form factors have been solved over different fit range combinations, the result we obtain is $F_{i,f}(Q^{2})$, where the ($A$ missing) index $f$ refers to which combined set of fit ranges were used. Since the extrapolation to $Q^{2}\rightarrow 0$ must be performed to compute the nucleon EDM, this must be performed for every $f(Q^{2})$ combination (analogous to $f(A)$ above). So in addition to above, we increase the number of fits to $\prod_{Q^{2}}F(Q^{2})$, where $F(Q^{2})$ is the number of fits computed for index $f(Q^{2})=1,2,\ldots,F(Q^{2})$.

Combining both these studies together, the resulting nucleon EDM has been computed using different fit ranges, indexed by $\frac{F^{p/n}_{3,f}(Q^{2}\rightarrow 0)}{2M_{N}}=d_{p/n,f}$. So to obtain a final result where the statistical uncertainty from the gauge fields and the systematic errors arising from the fit ranges can be combine into a single uncertainty, we extend the bootstrap samples which are already used to compute the statistical uncertainty $d_{p/n,B}$ where $B$ runs over $[1,N_{b}]\otimes[1,N_{F}]$ where $N_{f}$ is the number of fits which each have $N_{b}$ bootstrap samples.

B.1 Computational Viability

As one may notice, the above formulation is of order $O(A!)$, assuming a fixed number of fit ranges selected. A stochastic estimation of the fit range variation is highly recommended, which can be employed when solving the form factor eq. (83), as well as when taking the form factor $Q^{2}\rightarrow 0$ extrapolation.

At the form factor solving stage, this is employed by randomly selecting $N_{\chi}$ different fit range that satisfy the $\chi^{2}_{PDF}$ criterium. The resulting number of systems of equations to be solved are $N_{\chi}^{N_{A}}$.

For the form factor extrapolation in $Q^{2}\rightarrow 0$, a random selection of $N_{F}$ results of index $f(Q^{2})$ in $F_{N,f(Q^{2})}(Q^{2})$. The resulting number of fits to be performed using this estimation is $N_{F}^{N_{Q}}$ for $N_{Q}$ number of transfer momentum $Q^{2}$.

B.2 Results computed in this paper

The results computed in this paper use the fit criterium $\chi^{2}_{PDF}\in[0.5,1]$, and excluded fits of length 2. The cutoff for the number of fit ranges per ratio function is $N_{\chi}=4$ and the cutoff for the form factor extrapolation is $N_{Q}=4$ as well.

As for the the number of equations to solve, results at lattice $(\bm{q}/a)^{2}=1,4$ has $1024$ equations, $(\bm{q}/a)^{2}=2$ has $4096$ equations and $(\bm{q}/a)^{2}=3$ has $16384$ equations. Multiplying these numbers by $200$ bootstrap samples, will give the individual number of system of equations solved. Once this is complete, we avoid computing $\sim 70$ trillion equations by performing the stochastic estimate which only requires $256$ equations to solve for. Although it may seem the values for $N_{\chi}$ and $N_{F}$ are insufficient in size, the results shown in fig. 27 demonstrates minimal variation when analyzing each individual $\frac{F^{p/n}_{3,f}(Q^{2})}{2M_{N}}$ over different fit ranges $f$.