Contribution of the QCD Θ-term to nucleon electric dipole moment
Tanmoy Bhattacharya, Vincenzo Cirigliano, Rajan Gupta, Emanuele Mereghetti, Boram Yoon
LLA-UR-20-30515
Contribution of the QCD Θ -term to nucleon electric dipole moment Tanmoy Bhattacharya, ∗ Vincenzo Cirigliano, † Rajan Gupta, ‡ Emanuele Mereghetti, § and Boram Yoon ¶ Los Alamos National Laboratory, Theoretical Division T-2, Los Alamos, NM 87545 Los Alamos National Laboratory, Computer Computational and Statistical Sciences, CCS-7, Los Alamos, NM 87545 (Dated: January 19, 2021)We present a calculation of the contribution of the Θ-term to the neutron and proton electric dipolemoments using seven 2+1+1-flavor HISQ ensembles. We also estimate the topological susceptibilityfor the 2+1+1 theory to be χ Q = (66(9)(4) MeV) in the continuum limit at M π = 135 MeV. Thecalculation of the nucleon three-point function is done using Wilson-clover valence quarks. The (cid:8)(cid:8) CPform factor F is calculated by expanding in small Θ. We show that lattice artifacts introduce aterm proportional to a that does not vanish in the chiral limit, and we include this in our chiral-continuum fits. A chiral perturbation theory analysis shows that the N ( ) π ( ) state should providethe leading excited state contribution, and we study the effect of such a state. Detailed analysisof the contributions to the neutron and proton electric dipole moment using two strategies forremoving excited state contamination are presented. Using the excited state spectrum from fits tothe two-point function, we find d Θ n is small, | d Θ n | (cid:46) .
01 Θ e · fm, whereas for the proton we get | d Θ p | ∼ .
02 Θ e · fm. On the other hand, if the dominant excited-state contribution is from the Nπ state, then | d Θ n | could be as large as 0 .
05 Θ e · fm and | d Θ p | ∼ .
07 Θ e · fm. Our overall conclusion isthat present lattice QCD calculations do not provide a reliable estimate of the contribution of theΘ-term to the nucleon electric dipole moments, and a factor of ten higher statistics data are neededto get better control over the systematics and possibly a 3 σ result. PACS numbers: 11.15.Ha, 12.38.GcKeywords: neutron electric dipole moment, Θ-term, CP violation, lattice QCD, form factors
I. INTRODUCTION
The permanent electric dipole moments (EDMs) ofnondegenerate states of elementary particles, atoms andmolecules are very sensitive probes of CP violation ( (cid:8)(cid:8)
CP).Since the EDMs are necessarily proportional to theirspin, and under time-reversal the direction of spin re-verses but the electric dipole moment does not, a nonzeromeasurement confirms CP violation assuming CPT isconserved. Of the elementary particles, atoms and nucleithat are being investigated, the electric dipole momentsof the neutron (nEDM) and the proton (pEDM) are thesimplest quantities for which lattice QCD can provide thetheoretical part of the calculation needed to connect theexperimental bound or value to the strength of (cid:8)(cid:8)
CP in agiven theory [1, 2].EDMs can shed light on one of the deepest mysteries ofthe observed universe, the origin of the baryon asymme-try: the universe has 6 . +0 . − . × − baryons for everyblack body photon [3], whereas in a baryon symmetricuniverse, we expect no more than about 10 − baryonsand anti-baryons for every photon [4]. It is difficult toinclude such a large excess of baryons as an initial condi-tion in an inflationary cosmological scenario [5]. The wayout of the impasse lies in generating the baryon excess ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] ¶ Electronic address: [email protected] dynamically during the evolution of the universe. But,if the matter-antimatter symmetry was broken post in-flation and reheating, then one is faced with Sakharov’sthree necessary conditions [6] on the dynamics: the pro-cess has to violate baryon number, evolution has to oc-cur out of equilibrium, and charge-conjugation and CPinvariance have to be violated.CP violation exists in the electroweak sector of thestandard model (SM) of particle interactions due to aphase in the Cabibbo-Kobayashi-Maskawa (CKM) quarkmixing matrix [7], and possibly due to a similar phase inthe Pontecorvo–Maki–Nakagawa–Sakata (PNMS) matrixin the leptonic sector [8, 9]. The effect of these on nEDMand pEDM is, however, small: that arising from the CKMmatrix is about O (10 − ) e cm [10–12], much smallerthan the current 90% confidence level (CL) experimentalbound d n < . × − e cm [13], and than the reachof ongoing experiments, d n < . × − e cm at 90%confidence [15].In principle, the SM has an additional source of CPviolation arising from the effect of QCD instantons.The presence of these localized finite action nonper-turbative configurations in a non-Abelian theory leadsto inequivalent quantum theories defined over various‘Θ’-vacua [16, 17]. Because of asymptotic freedom,all nonperturbative configurations including instantons The slightly stronger 95% CL bounds d n < . × − e cm and d p < . × − e cm can be obtained from the experimentallimit on the Hg [14] EDM, assuming that nucleon EDMs arethe dominant contributions to the nuclear EDM. a r X i v : . [ h e p - l a t ] J a n are strongly suppressed at high temperatures [18, 19]where baryon number violating processes occur. Becauseof this, CP violation due to such vacuum effects doesnot lead to appreciable baryon number production [20].Nonetheless, understanding the contribution of such aterm to the nucleon EDM is very important for two rea-sons. First, the Θ term constitutes a ‘background’ contri-bution to all hadronic EDMs that needs to be understoodbefore one can claim discovery of new sources of CP vio-lation through nucleon or hadronic EDM measurements;and second, besides generating higher-dimensional CP-odd operators, new sources of CP-violation beyond theStandard Model (BSM) also generate a so-called ‘inducedΘ term’ [1, 21, 22] if one assumes that the Peccei-Quinnmechanism is at work [23]. Therefore, in the large classof viable models of CP violation that incorporate thePeccei-Quinn mechanism, quantifying the contributionof the induced Θ to the nucleon EDM (operationally, thecalculation is the same as in the first case) is essential tobound or establish such sources of CP violation.Until recently, the calculation of hadronic matrix ele-ments needed to connect nucleon EDMs to SM and BSMsources of CP violation relied on chiral symmetry sup-plemented by dimensional analysis [24–32] or QCD sumrules [1, 22, 33–36], both entailing large theoretical errors.Large-scale simulations of lattice QCD provide a first-principles method for calculating these matrix elementswith controlled uncertainties. Several groups have re-ported results of lattice QCD calculations of the neutronEDM induced by the QCD Θ term [37–44] and by higher-dimensional operators, such as the quark EDM [45, 46]and at a more exploratory level the quark chromo-EDM [47–49]. In this paper, we present a new calcu-lation of the contribution of the Θ-term to the nEDMand pEDM and show that the statistical and systematicuncertainties are still too large to extract reliable esti-mates.This paper is organized as follows: In Section II, we de-scribe our notation by introducing the Lagrangian with (cid:8)(cid:8) CP and the needed matrix elements. In Section III, wedescribe the decomposition of the matrix elements intothe electromagnetic form factors. Section IV providesthe lattice parameters used in the calculations. In Sec-tion V, we present the implementation of the gradientflow scheme, and in Sec.VI the calculation of the topo-logical susceptibility. Section VII describes the method-ology for extracting the (cid:8)(cid:8)
CP phase α for the ground statecreated by the nucleon interpolating operator used, fromthe two-point function. This phase controls the CP trans-formation of the asymptotic nucleon state. Section VIIIdescribes the calculation strategy for obtaining the formfactors when this phase α is nonzero and gives the formu-lae used to extract the (cid:8)(cid:8) CP form factor F from the ma-trix elements. In Section IX, we discuss the extraction of F ( q ) and the removal of the excited states contamina-tion. The extrapolation of F ( q ) to q = 0 is presentedin Sec. X. Section XI discusses the lattice-spacing arti-facts. Our results with the excited state spectrum taken from the two-point function are presented in Sec. XII andthose with an N π excited state in Sec. XIII. These resultsare compared to previous calculations in Section XIV.Conclusions are presented in Section XV. Further detailson the connection between Minkowski and Euclidean no-tation, the extraction of the form factors, the chiral ex-trapolation, excited-state contamination, and the O ( a )corrections in the Wilson-clover theory are presented infive appendices. II. THE QCD Θ -TERM QCD allows for the existence of a P and T (and (cid:8)(cid:8)
CPif CPT is conserved) violating dimension-four operator,i.e., the Θ-term. In its presence, the QCD Lagrangiandensity in Euclidean notation becomes L QCD −→ L (cid:8)(cid:8)
CPQCD = L QCD + i Θ G aµν (cid:101) G aµν π (1)where G aµν is the chromo-field strength tensor, (cid:101) G aµν = (cid:15) µνλδ G aλδ is its dual, and Θ is the coupling. G µν (cid:101) G µν is a total derivative of a gauge-variant current and itsspace-time integral gives the topological charge Q = (cid:90) d x G aµν (cid:101) G aµν π . (2)Non-zero values of Q are tied to the topological structureof QCD and the U (1) axial anomaly. In addition, higherdimension operators that arise due to novel (cid:8)(cid:8) CP couplingsat the TeV scale generate this term under renormaliza-tion in a hard cutoff scheme like lattice regularization orgradient flow [45]. Also, BSM models in which the Peccei-Quinn mechanism is operative induce such a term [1].Under a chiral transformation, one can rotate Θ intoa complex phase of the quark matrix and vice versa. Itis, therefore, necessary to work with the convention in-dependent Θ = Θ + Arg Det M q , which includes both,Θ from all sources and the overall phase of the quarkmatrix M q . Since, the argument of the determinant isill-defined when it is zero, all physical effects of Θ vanishin the presence of even a single massless quark flavor.If the overall Θ is nonzero, then this operator wouldinduce an nEDM d n of size d n = Θ X (3) X ≡ lim q → F ( q )2 M N Θ . (4) Throughout the paper, we work in Euclidean space, using q forthe Euclidean 4-momentum and Q for the topological charge.The gauge field includes a factor of the strong coupling, g , sothat the kinetic term is G aµν G aµν / g . Also, our conventions forconnecting the Euclidean and Minkowski metrics are given inAppendix A. Here X is obtained from the (cid:8)(cid:8) CP part of the matrix el-ement of the electromagnetic vector current within theneutron state in the presence of the Θ-term and F isthe (cid:8)(cid:8) CP violating form factor defined in Eq. (6). This isobtained, at the leading order, from the (cid:8)(cid:8)
CP part of thematrix element (cid:104) N | J EM µ | N (cid:105) (cid:12)(cid:12) Θ ≈ (cid:104) N | J EM µ | N (cid:105) (cid:12)(cid:12) Θ=0 − i Θ (cid:42) N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) J EM µ (cid:90) d x G aµν (cid:101) G aµν π (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N (cid:43) , (5)where we have assumed that the Θ-term is the onlysource of (cid:8)(cid:8) CP. In other words, X provides the connec-tion between the (cid:8)(cid:8) CP coupling (Θ) and the nEDM ( d n ).At present, the upper bound on the nEDM, | d n | < . × − e cm (90% CL) [13], is used along with an es-timate X ∼ (2 . ± . × − e cm [1] to set a limit onthe size of Θ (cid:46) − . This is an unnaturally small num-ber! One solution to this unnaturalness is the dynamicaltuning of Θ = 0 using the Peccei-Quinn mechanism [23].Our goal is to calculate X using lattice QCD, whichmultiplied by the cumulative value, Θ, from all sources(SM or BSM), gives the full contribution to nEDM fromthe dimension-4 G ˜ G operator in Eq. (1). Knowing X willallow current and future bounds on (or measured valueof) d n to more stringently constrain or pin down Θ.In the rest of the paper, all the analyses are carriedout assuming that the only (cid:8)(cid:8) CP coupling arises fromthe Θ-term, whose strength is Θ. Results are presentedfor Θ = 0 .
2, which we have checked is small enough sothat O (Θ ) corrections are negligible for all quantitiesof interest ( α and F defined later).The lattice calculation consists of the evaluation of theconnected and disconnected diagrams shown in Fig. 1.The disconnected diagram gets contributions from all γ µ Q γ µ Q FIG. 1: The connected (left) and disconnected (right) dia-grams with the insertion of the bilinear vector current (redfilled circle) in the nucleon two-point function. The signal isgiven by the correlation between this 3-point function and thetopological charge shown by the filled yellow circle. The Peccei-Quinn mechanism relaxes Θ dynamically to Θ ind , thepoint where the effective potential achieves its minimum. In theabsence of other sources of CP violation in the theory, Θ ind = 0. quark flavors in the loop—but their contributions tothe CP-conserving form-factors of the vector current aresmall [50]. In this work, we assume the same holds forthe CP-violating ones and neglect the contribution to theelectric dipole moment coming from these diagrams.
III. FORM FACTOR OF THEELECTROMAGNETIC CURRENT
The parameterization of the matrix element of the elec-tromagnetic current, J EM µ ( q ), defined in Eq. (5), withinthe nucleon state in terms of the most general set of formfactors consistent with the symmetries of the theory is (cid:104) N ( p (cid:48) , s (cid:48) ) | J EM µ | N ( p, s ) (cid:105) Θ (cid:54) CP = u N ( p (cid:48) , s (cid:48) ) (cid:34) γ µ F ( q )+ 12 M N σ µν q ν (cid:16) F ( q ) − iF ( q ) γ (cid:17) + F A ( q ) M N ( /qq µ − q γ µ ) γ (cid:35) u N ( p, s ) , (6)where M N is the nucleon mass, q = p (cid:48) − p is theEuclidean 4-momentum transferred by the electromag-netic current, σ µν = ( i/ γ µ , γ ν ], and u N ( p, s ) repre-sents the free neutron spinor of momentum p and spin s obeying ( i/p + M N ) u N ( p, s ) = 0, with γ implementingthe asymptotic (i.e., free) parity operation. Through-out, we work in Euclidean space and refer the readerto Appendix A for details on our conventions. F and F are the Dirac and Pauli form factors, in terms ofwhich the Sachs electric and magnetic form factors are G E = F − ( q / M N ) F and G M = F + F , respec-tively. The anapole form factor F A and the electricdipole form factor F violate parity P; and F violatesCP as well. The zero momentum limit of these form fac-tors gives the charges and dipole moments: the electriccharge is G E (0) = F (0), the magnetic dipole moment is G M (0) / M N = ( F (0) + F (0)) / M N , and the EDM isdefined in Eq. (4).In all the discussions in this paper, the current J EM µ used is the renormalized local vector current Z V (cid:80) i e i ψ i γ µ ψ i , where e i is the electric charge of a quarkwith flavor i . The renormalization is carried out by tak-ing ratios of all three-point fermion correlators with thelattice estimate of the vector charge, g V ≡ /Z V , whichis given by the forward matrix element of ψ i γ ψ i . Theseratios are constructed with identical source, sink, andcurrent insertion positions and within the single jack-knife loop used for the statistical analysis of the data to We emphasize that we use q for the Euclidean four-momentum-squared that is denoted by Q in our previous work and through-out the literature. As noted in the Appendix A, it is the negativeof the Minkowski four-momentum-squared. Ensemble a M valπ L × T M val π L τ /a aM N N conf Confs. N HP N LP χ / Q ID [fm] [MeV] Per Bin [MeV] a m
310 0 . . .
8) 24 ×
64 4 . { , , } . ,
052 64 ,
832 145.9(2.7) a m
220 0 . . .
9) 32 ×
64 4 . { , , } . ,
000 64 ,
000 145.3(2.4) a m L . . .
7) 40 ×
64 5 . { , , } . ,
000 128 ,
000 141.3(2.5) a m
310 0 . . .
8) 32 ×
96 4 . { , , } . ,
784 140 ,
544 129.5(2.3) a m
220 0 . . .
8) 48 ×
96 4 . { , , } . ,
844 123 ,
008 115.0(2.2) a m
130 0 . . .
0) 64 ×
96 3 . { , , } . ,
156 164 ,
992 106.8(1.7) a m
310 0 . . ×
144 4 . a m
220 0 . . ×
144 4 . a m
135 0 . . .
4) 96 ×
192 3 . { , , , } . ,
812 28 ,
992 89.3(2.8)TABLE I: Lattice parameters, nucleon mass M N , number of configurations analyzed, and the total number of high precision(HP) and low precision (LP) measurements made. We also give the bin size (Confs. per bin) used in the statistical analysis oftwo- and three-point functions. The last column gives the topological susceptibility χ Q calculated at flow time τ gf = 0 .
68 fmand with a bin size of 20 configurations. The ensembles a m
310 and a m
220 have been used only for the calculation of χ Q ,and 861 configurations were used to calculate χ Q on the a m
135 ensemble. take advantage of error reduction due to correlated fluc-tuations. IV. LATTICE PARAMETERS
We present results on seven ensembles, whose param-eters are defined in Table I. These were generated by theMILC collaboration [51] using 2+1+1-flavors of highlyimproved staggered quarks (HISQ) action. For the con-struction of the nucleon correlation functions we use theclover-on-HISQ formulation that has been used exten-sively by us in the calculation of the nucleon charges andform factors as described in Refs. [52, 53]. These ensem-bles cover three values of the lattice spacing, a ≈ . .
09 and 0 .
06 fm and three values of the pion mass M π ≈ ,
220 and 130 MeV. Further details of the lattice pa-rameters and methodology, statistics, and the interpolat-ing operator used to construct the nucleon 2- and 3-pointcorrelation functions can be found in Refs. [52, 53].
V. TOPOLOGICAL CHARGE UNDERGRADIENT FLOW
We calculate the topological charge using the gradi-ent flow scheme to implement operator renormalizationand to reduce lattice discretization effects [41, 54]. Theprimary advantage of the scheme is that at finite flow This forward matrix element has very small excited state con-tamination and, therefore, does not affect our excited state fitsat this level of precision. times , i.e., for τ gf >
0, the flow time provides an ul-traviolet cutoff, and the continuum limit, a →
0, of alloperators built solely from gauge fields is finite. More-over, since topological sectors arise dynamically as wetake the continuum limit, the gradient flowed topologicalcharge takes on integer values, and no renormalization isneeded to convert it to a scheme that preserves this prop-erty; in particular, correlators of the topological chargeare flow-time independent [54].These statements are, however, not true at finite latticespacing and volume. At small τ gf , we get O ( a /τ ) arte-facts. In Fig. 2, we show the distribution of the topologi-cal charge Q as a function of the flow time τ gf in physicalunits. Its distribution has stabilized by τ gf = 0 .
24 fm forthe a = 0 .
12 fm ensembles, and by τ gf = 0 .
17 fm for the a = 0 .
09 and 0 .
06 fm ensembles. The large values of Q that form the long tail of the distribution at τ gf = 0 aresmoothed out, indicating that they are lattice artifacts.In Fig. 3, we show the distribution of the differencefrom the nearest integer. This distribution stabilizesmore slowly and it is only by τ gf = 1 .
31 fm ( τ gf = 0 .
76 fm)on the a ≈ .
12 fm ( a ≈ .
09 and 0 .
06 fm) ensembles thatthe charges are close to integers. The relevant distribu-tion important for the calculation of the nucleon corre-lation functions is, however, likely to be the distributionof Q shown in Fig. 2. To explore this, we show in Fig. 4the value of F as a function of τ gf for the a ≈ . .
09 fm ensembles, and find that indeed the corre-lation functions, and thus F , do stabilize early but the We use the notation τ gf ≡ √ t for the flow time, where t isthe parameter in the flow equations in Ref. [54]. We used theRunge-Kutta integrator given in that reference for integratingthe flow equations, with a step size of 0 .
01. Changing the stepsize to 0 .
002 changed the results on topological susceptibility byless than 0.2%.
75 50 25 0 25 50 75Q0.000.010.020.030.040.050.06 P r o b a b ili t y a12m310 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.68 fm gf =0.86 fm gf =1.04 fm gf =1.40 fm
40 20 0 20 40Q0.000.010.020.030.040.050.060.07 P r o b a b ili t y a09m310 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.51 fm gf =0.68 fm gf =0.76 fm
75 50 25 0 25 50 75Q0.000.010.020.030.04 P r o b a b ili t y a12m220 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.68 fm gf =0.86 fm gf =1.04 fm gf =1.40 fm
40 20 0 20 40Q0.000.010.020.030.040.050.06 P r o b a b ili t y a09m220 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.51 fm gf =0.68 fm gf =0.76 fm
75 50 25 0 25 50 75Q0.0000.0050.0100.0150.0200.0250.030 P r o b a b ili t y a12m220L gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.69 fm gf =0.86 fm gf =1.04 fm gf =1.41 fm
40 20 0 20 40Q0.000.010.020.030.04 P r o b a b ili t y a09m130 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.51 fm gf =0.68 fm gf =0.76 fm
40 20 0 20 40Q0.000.010.020.030.040.050.06 P r o b a b ili t y a06m135 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.51 fm gf =0.68 fm gf =0.76 fm FIG. 2: The distribution of the topological charge Q as afunction of the flow time τ gf . The panels on the left (right)show data for the a = 0 .
12 fm ( a = 0 .
09 fm) ensembles P r o b a b ili t y a12m310 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.68 fm gf =0.86 fm gf =1.04 fm gf =1.40 fm P r o b a b ili t y a09m310 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.51 fm gf =0.68 fm gf =0.76 fm P r o b a b ili t y a12m220 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.68 fm gf =0.86 fm gf =1.04 fm gf =1.40 fm P r o b a b ili t y a09m220 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.51 fm gf =0.68 fm gf =0.76 fm P r o b a b ili t y a12m220L gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.69 fm gf =0.86 fm gf =1.04 fm gf =1.41 fm P r o b a b ili t y a09m130 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.51 fm gf =0.68 fm gf =0.76 fm P r o b a b ili t y a06m135 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.51 fm gf =0.68 fm gf =0.76 fm FIG. 3: The panels show the distribution of the difference, Q − Q int , of the measured Q from the nearest integer Q int . -0.6-0.4-0.20.00.20.40.60.8 0 0.2 0.4 0.6 0.8 1 F ~ , n / Θ – τ gf (fm) a12m310, q = 0.18 GeV a12m220L, q = 0.07 GeV -0.6-0.4-0.20.00.20.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 F ~ , n / Θ – τ gf (fm) a09m310, q =0.18 GeV a09m220, q =0.09 GeV a09m130, q =0.05 GeV FIG. 4: Data for (cid:101) F ,n / Θ, defined in Eq. (28), at the smallest value of q , respectively, on the a
12 (left panel) and the a
09 (rightpanel) ensembles. The estimates show no significant change after τ gf ≈ . a
09 ensembles and τ gf ≈ . a
12 ensembles. τ gf required for the coarser lattices is longer. Thus, tobe conservative, the results presented below are obtainedwith flow times τ gf ( a
06) = 0 .
68 fm, τ gf ( a
09) = 0 .
68 fmand τ gf ( a
12) = 0 .
86 fm respectively.In Fig. 5, we show the distribution of the nearestinteger, Q int , to the topological charge at τ gf ≈ . τ gf = 0 .
76 fm) on the a ≈ .
12 fm ( a ≈ .
09 and 0 .
06 fm)ensembles, by which time the Q int identified with agiven configuration has stabilized. This distribution isapproximately symmetric about zero as expected since (cid:104) Q (cid:105) = 0, and no gaps are visible in the distribution. InFig. 6, we show the autocorrelation function of Q versus
30 20 10 0 10 20 30NearestInteger(Q)0.000.010.020.030.040.050.06 P r o b a b ili t y a12m310 gf =1.40 fm
20 10 0 10 20NearestInteger(Q)0.000.010.020.030.040.050.060.07 P r o b a b ili t y a09m310 gf =0.76 fm
40 20 0 20 40NearestInteger(Q)0.000.010.020.030.04 P r o b a b ili t y a12m220 gf =1.40 fm
30 20 10 0 10 20 30NearestInteger(Q)0.000.010.020.030.040.050.06 P r o b a b ili t y a09m220 gf =0.76 fm
40 20 0 20 40NearestInteger(Q)0.0000.0050.0100.0150.0200.0250.0300.035 P r o b a b ili t y a12m220L gf =1.41 fm
40 20 0 20 40NearestInteger(Q)0.000.010.020.030.04 P r o b a b ili t y a09m130 gf =0.76 fm
20 10 0 10 20NearestInteger(Q)0.000.010.020.030.040.05 P r o b a b ili t y a06m135 gf =0.76 fm FIG. 5: The distribution of the nearest integer charge, Q int , associated with a given configuration at τ gf ≈ . a
12 ensembles) and 0 .
76 fm ( a
09 and a
06 ensembles), bywhich time the Q int identified with a given configurationhas stabilized. t acf A u t o c o rr e l a t i o n a12m310 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.68 fm gf =0.86 fm gf =1.04 fm gf =1.40 fm 0 5 10 15 20 t acf A u t o c o rr e l a t i o n a09m310 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.51 fm gf =0.68 fm gf =0.76 fm0 5 10 15 20 t acf A u t o c o rr e l a t i o n a12m220 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.68 fm gf =0.86 fm gf =1.04 fm gf =1.40 fm 0 5 10 15 20 t acf A u t o c o rr e l a t i o n a09m220 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.51 fm gf =0.68 fm gf =0.76 fm0 5 10 15 20 t acf A u t o c o rr e l a t i o n a12m220L gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.69 fm gf =0.86 fm gf =1.04 fm gf =1.41 fm 0 5 10 15 20 t acf A u t o c o rr e l a t i o n a09m130 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.51 fm gf =0.68 fm gf =0.76 fm0 5 10 15 20 t acf A u t o c o rr e l a t i o n a06m135 gf =0.00 fm gf =0.17 fm gf =0.34 fm gf =0.51 fm gf =0.68 fm gf =0.76 fm FIG. 6: The autocorrelation function for different valuesof the flow time. The data show that the long τ gf behaviorstabilizes by τ gf = 0 .
34 fm in all cases. Q a09m130 gf =0.68 fm0 200 400 600 800Configurations20020 Q a06m135 gf =0.68 fm FIG. 7: The time history of Q on the a m
130 (upper) and a m
135 (lower) ensembles at τ gf = 0 .
68 fm. No long timefreezing of the topological charge is observed. the flow time. The data show no significant change after τ gf (cid:38) . Q in any of the ensembles analyzed as illustratedusing the a m
130 and a m
135 ensembles at flow time τ gf = 0 .
68 fm in Fig. 7. The autocorrelation is lessthan about 10 configurations for all but the a m VI. TOPOLOGICAL SUSCEPTIBILITY
The topological susceptibility χ Q is defined as χ Q = (cid:90) d x (cid:104) Q ( x ) Q (0) (cid:105) . (7)Its value in the pure gauge theory, χ quenched Q , is relatedto the mass of the η (cid:48) meson in a theory with N f lightflavors in the chiral limit via the axial anomaly, viz., the gf (fm) / [ M e V ] a12m220L 0.2 0.3 0.4 0.5 0.6 0.7 gf (fm) / [ M e V ] a09m220 FIG. 8: Illustration of the flow-time dependence of the topological susceptibility at small flow times showing that it is almostindependent of the flow time when the flow time is much larger than the lattice spacing. -2-1 0 1 2 3 4 5 0 0.2 0.4 0.6 0.8 1 χ / ( τ g f ) – χ / ( L / √ ) [ M e V ] τ gf2 / L a12m220 a12m220L FIG. 9: Comparison of the flow-time dependence of the topo-logical susceptibility at large flow times on two ensembles dif-fering only in lattice volume, showing that the dependence isa finite size effect.
Witten-Veneziano relation [55, 56] M η (cid:48) ≈ N f F π χ quenched Q , (8)where F π is the pion decay constant in the conventionwhere its physical value is about 93 MeV. FollowingRef. [57], we can include the effects of the quark masses.Including SU (3) breaking at leading order in χ PT butneglecting the heavier quarks gives χ quenched Q ≈ F π ( M η (cid:48) − M η )6 (cid:32)(cid:114) − δ Kπ δ Kπ (cid:33) ,χ quenched Q ≈ F π ( M η (cid:48) − M η )6 (cid:32) M η − M K M η (cid:48) − M η (cid:33) , (9)where δ Kπ ≡ ( M K − M π ) / ( M η (cid:48) − M η ) is an SU(3)breaking ratio. The two expressions, which can be de-rived independently, give χ quenched Q ≈ (172 MeV) and (179 MeV) respectively, thus quantifying the accuracyof the expansion.With dynamical fermions, however, the susceptibilityshould vanish in the chiral limit. For SU ( N f ) flavorgroup with finite but degenerate quark masses, it shouldbehave as [58–60]:1 χ Q ≈ χ quenched Q + 2 N f M π F π . (10)For N f = 2 light flavors and the strange quark, butneglecting the heavier quarks that give negligible cor-rections, leading order chiral perturbation theory ( χ PT)modifies this to1 χ Q ≈ χ quenched Q + 4 M π F π (cid:18) − M π M η (cid:19) − . (11)We calculate χ Q on the 2+1+1 flavor HISQ ensem-bles, which are O ( a ) improved. The results are given inTable I. In addition to the seven ensembles used to calcu-late F , we include data from the a m
310 and a m F .As discussed in Section V, the topological suscepti-bility at finite flow time needs no renormalization, andshould be independent of flow time up to O ( a /τ ) ef-fects. As shown in Fig. 8, this is true up to a small,almost linear, downward drift with increasing flow time.In Fig. 9, we compare the results on a m
220 and a m L ensembles, and show that this is a τ /L ef-fect, where L is the lattice size. At the flow times and For asymmetric lattices like ours, we expect the smaller spatialextent to dominate the finite volume effect. volumes we use in the calculation, this is a small effectand therefore neglected.To obtain χ Q at M π = 135 MeV and a = 0, we use thefit ansatz χ Q ( a, M π ) = c a + c M π + c a M π , (12)which assumes χ Q is zero in the chiral-continuum limit.We do not find a viable χ /dof on including all ninedata points. Reasonable fits are found on neglecting(i) all three a ≈ .
12 fm points and (ii) all three a ≈ .
12 fm and the a m
310 point. These two fitsgive χ Q = [70(6) MeV] and χ Q = [63(9) MeV] , re-spectively, at M π = 135 MeV. We take the average χ Q = [66(9)(4) MeV] as our best estimate, the larger ofthe two errors and an additional systematic uncertainty,which is half the difference. These results are in goodagreement with the expected value, (79 MeV) , obtainedusing the physical meson masses and decay constants inEqs. (9) and (11). The data and the fit case (i) are shownin Fig. 10. VII. CALCULATION OF THE (cid:8)(cid:8)
CP PHASE α In a field theory in which parity is not conserved, thedefinition of parity of a composite state, e.g., the neutronstate, needs care [22, 33, 37]. To explain this, we startwith the most general spectral decomposition of the time-ordered 2-point nucleon correlator (cid:104) Ω |T N ( p , τ ) N ( p , | Ω (cid:105) = (cid:88) i, s e − E i τ A ∗ i A i M s i , (13)where A i is the amplitude for creating state i , E i is itsenergy, the Euclidean time τ is the separation betweenthe source and the sink, and, for notational convenience,we are assuming a discrete spectrum. A common choiceon the lattice of the neutron interpolating operator N is N ≡ (cid:15) abc [ d a T Cγ γ u b ] d c , (14)where C = γ γ (the sign is conventional and does notaffect the nucleon correlators we study; see Appendix Afor details of our convention) is the charge conjugationmatrix, a , b , c are the color indices and u , d are the quarkflavors. The 4 × M s i in Eq. (13) dependson the state and the momentum p . Its most general formconsistent with Lorentz covariance is (cid:88) s M s i = e iα i γ ( − i/p i + M i )2 E pi e iα ∗ i γ (15) ≡ e iα i γ (cid:88) s u iN ( p , s ) u iN ( p , s ) e iα ∗ i γ , (16) Up to a possible extra factor of γ , which, however, is prohibitedby PT symmetry in our calculations. where p i ≡ iE i . It is clear that because of the presenceof the phases α i , the parity operator that transformsthe spinor associated with the i th asymptotic stateis P α i ≡ e iα i γ P e − iα i γ , where P ≡ ηγ is the usualparity operator for a particle with intrinsic parity η .The phases α i depend on the realization of discretesymmetries: If the interpolating field is chosen such that P implements parity in the free theory, Im α i = 0 for aPT symmetric theory, Re α i = 0 for the CP symmetrictheory, α i = 0 for a P symmetric theory. For our case ofonly (cid:8)(cid:8) CP, all α i are, therefore, real, which will be implicitexcept in Appendix B. It is important to note that thevalue of α i depends on the interpolating operator N ,the state, and the source of (cid:8)(cid:8) CP. Its value for the groundstate can be extracted from the large τ behavior of theimaginary part of the nucleon 2-point function. Consider r α ( τ ) ≡ Im C P ( τ )Re C ( τ ) (17) ≡ Im Tr (cid:2) γ (1 + γ ) (cid:104) N ( τ ) N (0) (cid:105) (cid:3) Re Tr (cid:2) (1 + γ ) (cid:104) N ( τ ) N (0) (cid:105) (cid:3) (18)= (cid:80) i M i sin(2 α i ) |A i | / (2 E i ) e − E i τ (cid:80) i ( E i + M i cos(2 α i )) |A i | / (2 E i ) e − E i τ . (19)Keeping only the first two states one gets r α ( τ ) ≈ M sin(2 α ) E + M cos(2 α ) (20) × M E M E sin(2 α )sin(2 α ) | ˜ A | e − ( E − E ) τ ( E + M cos(2 α )) E ( E + M cos(2 α )) E | ˜ A | e − ( E − E ) τ , where ˜ A i = A i / A . At zero three-momentum ( E i = M i )the above expression simplifies to r α ( τ ) ≈ tan α × sin(2 α )sin(2 α ) | ˜ A | e − ( M − M ) τ cos ( α )cos ( α ) | ˜ A | e − ( M − M ) τ . (21)The data for r α versus τ are shown in Fig. 11 for allseven ensembles. The α for the ground state obtainedfrom the two-state fit agrees with the plateau at large τ ,where the lowest state dominates, and is independent ofmomentum. VIII. THREE-POINT FUNCTIONS IN THEPRESENCE OF THE PHASE α In the presence of the phase α N ≡ α corresponding tothe ground-state nucleon [47], the most straightforwardway to extract the matrix element of the electromagneticcurrent J EM µ within the neutron ground state in the pres-ence of (cid:8)(cid:8) CP is to calculate the correlation function e − iα N γ (cid:104) Ω | N ( p (cid:48) , τ ) J EM µ ( q , t ) N ( p , | Ω (cid:105) (cid:12)(cid:12) (cid:54) CP e − iα N γ -1012345678 0 0.005 0.01 0.015 χ Q × - ( M e V ) a (fm ) a12m310a12m220a12m220La09m310 a09m220a09m130a06m310a06m220 a06m135Extrap χ /dof = 1.2 -1012345678 0 0.03 0.06 0.09 0.12 χ Q × - ( M e V ) M π (GeV ) a12m310a12m220a12m220La09m310 a09m220a09m130a06m310a06m220 a06m135Extrap χ /dof = 1.2 FIG. 10: Fits to the data for the topological susceptibility, χ Q , using the ansatz given in Eq. (12). a12m310 χ /dof = 0.64 r α / Θ – τ p = 0.00 GeV p = 0.34 GeV a09m310 χ /dof = 0.53 r α / Θ – τ p = 0.00 GeV p = 0.36 GeV a12m220 χ /dof = 1.61 r α / Θ – τ p = 0.00 GeV p = 0.21 GeV a12m220L χ /dof = 0.47 r α / Θ – τ p = 0.00 GeV p = 0.13 GeV a09m220 χ /dof = 0.15 r α / Θ – τ p = 0.00 GeV p = 0.17 GeV a09m130 χ /dof = 0.97 r α / Θ – τ p = 0.00 GeV p = 0.10 GeV a06m135 χ /dof = 0.05 r α / Θ – τ p = 0.00 GeV p = 0.10 GeV FIG. 11: The extraction of the phase α N / Θ with Θ = 0 . r α defined in Eq. (21). It is a Lorentz scalar and independent of the momentum as confirmed by the latticedata. The χ /dof values presented are from fully correlated fits, except for the case of the a m
220 and a m
135 ensembleson which we use uncorrelated fits to avoid instabilities. ∝ ( − i/p (cid:48) + M N ) O µ ( − i/p + M N ) , (22)where p (cid:48) ≡ p + q and O µ ≡ γ µ F + 12 M N σ µν q ν ( F − iF γ )+ F A M N ( /qq µ − q γ µ ) γ . (23)Here, the current J EM µ is inserted at times t between theneutron source and sink operators located at time 0 and τ , and a sum over the spin labels is implicit. We also as-sume that t and τ are large enough that only the groundstate dominates the correlation function. This form re-sults from the realization that γ remains the parity op-erator for the ground state nucleon when working withthe interpolating field defined to be e − iα N γ N instead of N in all correlation functions.This approach, however, requires, evaluating the full4 × P pt ≡
12 (1 + γ )(1 + iγ γ ) , (24)so the contribution of a nonzero α N has to be incorpo-rated at the time of the decomposition of the matrix ele-ment into the form factors. As discussed in Appendix B,by taking a suitable ratio of 3- and 2-point functions,one can isolate the four-vector V µ encoding the nucleonground state contribution to the matrix element of theelectromagnetic current, V µ ≡
14 Tr (cid:2) e iα N γ P pt e iα N γ ( − i/p (cid:48) + M N ) O µ ( − i/p + M N ) (cid:3) , (25)where O µ is given in Eq. (23). The full expressions for V , , , , along with a general strategy for extracting F ,from the four coupled complex equations is given in Ap-pendix B.To extract F , the (cid:8)(cid:8) CP part of the three-point func-tions, a very significant simplification of the analysis andimprovement in the signal is achieved by subtracting theΘ = 0 contribution from each component of the currentin Eq. (25) before making the excited state fits and de-composing the resulting ground state matrix element interms of form factors. This is implemented by analyzingthe ground state contribution in terms of the combina-tion ¯ V µ = V µ (Θ) − V µ (0). Working to first order in Θ,and recalling that s α N ≡ sin α N cos α N ∼ α N ∼ O (Θ),and F ∼ O (Θ), the expressions for the ground state con-tributions of the three-point functions ¯ V , , , in terms ofform factors simplify to¯ V = − q q G , (26a)¯ V = − q q G , (26b)¯ V = 12 (cid:16) M N ( E N − M N ) s α N G − q G (cid:17) , (26c) ¯ V = i (cid:16) q ( E N + M N ) G − q M N s α N G (cid:17) = iq M N (cid:16) ( E N + M N )2 M N F − s α N G E (cid:17) , (26d)where G = F + F and G = F + s α N F . We solve theabove system for G and G . At q = 0 there is a furthersimplification because G (0) = Q N + F (0) where Q N isthe nucleon charge. With this, we get F (0) = G (0) − s α N ( G (0) − Q N ) . (27)Though the nucleon anomalous magnetic moment G (0) − Q N = F (0) ≡ κ N has been measured very pre-cisely, the largest contribution to G comes from s α N F ,and the statistical error is much smaller when extrapo-lating G ( q ) − s α N ( G ( q ) − Q N ), rather than extrap-olating only G ( q ) and then combining it with s α N µ N to get the right hand side of Eq. (27). Also, note that G ( q ) can be obtained uniquely from ¯ V and ¯ V for anumber of values of q , which provides a useful check.One can extend Eq. (27) to define˜ F ( q ) ≡ G ( q ) − s α N (cid:0) G ( q ) − Q N (cid:1) . (28)To get F (0) = ˜ F (0), we find better control by extrapo-lating ˜ F ( q ) to q → V and V ) significantly reduces the statistical errorsand improves the analysis of excited state contamination(ESC) discussed next. IX. REMOVING ESC IN F In order to extract the ground state contribution ¯ V µ from lattice data on the ratio R µ ( τ, t, q ) of three- andtwo-point functions defined in Eq. (B7), we need to re-move all excited states that make a significant contribu-tion.We have analyzed data on R µ ( τ, t, q ) in terms of a two-state fit, following two strategies. In the first, we havetaken the first excited-state energies from a three-state fitto the two-point function. In the second strategy, we haveset the first excited-state energy to the non-interactingenergy of the N π state, motivated by the χ PT expecta-tion that the leading excited state is the
N π state, withamplitude of the same size as the ground state contri-bution (see Appendix D for more details). In Fig. 13we compare the two strategies for Im( R ( τ, t, q )). The χ /dof of the fits are similar for the two cases on allthree ensembles, but the ground state estimate is vastlydifferent and thus the contribution to the nEDM. With1 τ =8 a , t=4 a , q =0.18 GeV F , n / Θ – Θ – OriginalVarianceReduction -0.004-0.002 0 0.002 0.004 0.006 0.008 -6 -4 -2 0 2 4 6 a09m310; q =(1,1,1)2 π / La With Θ –=0 subtractionWith momentum average R e ( R ( τ , t , q ) ) t τ =08 τ =10 τ =12 τ =14 -0.004-0.002 0 0.002 0.004 0.006 0.008 -6 -4 -2 0 2 4 6 a09m310; q =(1,1,1)2 π / La With Θ –=0 subtractionWithout momentum average R e ( R ( τ , t , q ) ) t τ =08 τ =10 τ =12 τ =14 FIG. 12: The improvement in signal under subtraction of the Θ = 0 contribution and averaging over equivalent momenta. Thepanel on the left shows, using data from the a m
310 ensemble, (i) the improvement in F ,n as Θ → R ( τ, t, q ) with both the Θ = 0 and momentum averagingon the a m
310 ensemble with Θ = 0 . q = (1 , , π/La , while that on the right is without averaging over equivalentmomenta. the current data, picking between them is the key un-resolved challenge for this calculation. The very largeextrapolation for τ → ∞ in the N π case, however, leadsus to question whether a two-state fit is sufficient if the
N π state is included and whether a similar effect mightcontaminate our extraction of α N . We therefore firstperform the analysis taking the excited state energy, E ,from a three-state fit to the two-point function and returnto an analysis including a N π state in Sec. XIII.A second issue arising from the small signal in F isthat two-state fits to many of the correlation functionswith the full covariance matrix are unstable with respectto variations in the values of τ and t skip , the number ofpoints skipped in the fits adjacent to the source and sinkfor each τ . Examples of this are shown in Fig. 14 forRe( R ( τ, t, q )). This has two consequences for the anal-ysis. First, we have carried out the final analysis usingonly the diagonal elements of the covariance matrix. Wehave, however, checked that in cases where fully covariantfits are possible, the two results are consistent. Since weuse uncorrelated fits for removing excited-state contam-ination, we do not quote a χ /dof for these fits. Second,the system of four equations, Eqs. (26), over determines G and G . While we solve the full set of equations asexplained in appendix B, the data from Re( R , ( τ, t, q )),which have poor signal, do not make a significant con-tribution. We have checked this by removing them fromthe analysis and the results are essentially unchanged,i.e., the results are dominated by Re( R ( τ, t, q )) andIm( R ( τ, t, q )). X. EXTRAPOLATION OF F ( q ) TO q → The ansatz used to extrapolate F ( q ) to q → F ( q ), de-fined in Eq. (28), instead of F ( q ) as they are consistentto leading order and the extraction of ˜ F ( q ) is bettercontrolled. We examine three fits based on Eq. (C1): • Linear: the quantities d i and S (cid:48) i are free parameters and H i is set to zero. • χ PT: Only d i is a free parameter, S (cid:48) i are given inEq. (C11), g in Eq. (C7), and the H i in Eq. (C13). • χ PTg0: Same as χ PT except g is left as a freeparameter.The data and fits for the neutron and proton are pre-sented in Figs. 15 and 16. The data are, within errors,flat in all cases and the extrapolated values from the threetypes of fits are consistent. Since in most cases, we havereliable data at only three values of q , we take the finalresult from the χ PT fit. At the end, we will take thedifference between the Linear and χ PT fits to estimatethe associated systematic uncertainty.
XI. ADDITIONAL O ( a ) ARTIFACTS
Before performing a chiral-continuum extrapolation ofthe results, in this section we justify our continuumextrapolation formula for d n (Θ) that includes an M π -independent term that does not vanish in the chiral limit,i.e., a term proportional to am q .There are multiple sources of O ( a ) corrections that weneed to consider. First, since our clover coefficient c SW is set to its tadpole-improved tree-level value, the ac-tion, and hence all matrix elements, have residual O ( α s a )corrections. Because of the use of smeared gauge fields,however, the tadpole-improved tree-level approximationis extremely good, and these are expected to be tiny ef-fects. Second, the vector current we insert is not im-proved [62], and, hence, we expect its renormalizationcoefficient to have O ( am q ) corrections. Such multiplica-tive terms, however, are unimportant near the chiral-continuum limit, where the (cid:8)(cid:8) CP form factors vanish. Athird source of O ( a ) effects is the required improvementof the vector current by an O ( am q ) mixing with thederivative of the tensor current, which can give rise toa nonzero F , but only in the presence of CP violation2 a09m310; q =(0,0,1)2 π / La χ /dof = 1.22 I m ( R ( τ , t , q ) ) t τ =08 τ =10 τ =12 τ =14 -0.006-0.004-0.002 0 0.002 0.004 0.006 -6 -4 -2 0 2 4 6 a09m310; q =(0,0,1)2 π / La χ /dof = 1.28 I m ( R ( τ , t , q ) ) t τ =08 τ =10 τ =12 τ =14 -0.0005 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 -6 -4 -2 0 2 4 6 a09m220; q =(0,0,1)2 π / La χ /dof = 1.14 I m ( R ( τ , t , q ) ) t τ =08 τ =10 τ =12 τ =14 -0.008-0.006-0.004-0.002 0 0.002 0.004 -6 -4 -2 0 2 4 6 a09m220; q =(0,0,1)2 π / La χ /dof = 1.29 I m ( R ( τ , t , q ) ) t τ =08 τ =10 τ =12 τ =14 -0.001 0 0.001 0.002 0.003 0.004 0.005 0.006 -6 -4 -2 0 2 4 6 a09m130; q =(0,0,1)2 π / La χ /dof = 1.04 I m ( R ( τ , t , q ) ) t τ =08 τ =10 τ =12 τ =14 -0.03-0.025-0.02-0.015-0.01-0.005 0 0.005 0.01 -6 -4 -2 0 2 4 6 a09m130; q =(0,0,1)2 π / La χ /dof = 0.73 I m ( R ( τ , t , q ) ) t τ =08 τ =10 τ =12 τ =14 FIG. 13: Comparison of the two-state fit to the ratio Im( R ( τ, t, q )) defined in Eqs. (B7) with the first excited-state energiestaken from a three-state fit to the two-point function (left panels) and set equal the non-interacting energy of the Nπ state(right panels). The data for the three ensembles with a ≈ .
09 fm are shown in the three rows. The χ /dof of the two setsof fits are comparable, but the extrapolated ground sate value (solid black line) is vastly different. The data are shown for q = (0 , , π/La and the four largest values of τ . All data are with Θ = 0 . in the theory. Since the topological charge does not in-troduce CP violation in the chiral limit, we would expectthe behavior of d n to be dominantly O ( a ) in the chirallimit, if these were the only O ( a ) effects.In Appendix E, we analyze the Wilson-clover theorybased on the framework of a continuum EFT for thelattice action and the axial Ward Identities. Following Refs. [63–65], we show that the topological charge gives O ( a ) (cid:8)(cid:8) CP corrections, and identify this as effectively dueto the insertion of the isoscalar quark chromo-EDM op-erator, which the topological term can mix with. Sincethis term is expected to survive in the chiral limit, weinclude an O ( am q ) term in our chiral continuum fits.3 -0.002-0.0015-0.001-0.0005 0 0.0005 0.001 -6 -4 -2 0 2 4 6 a09m310; q =(0,1,1)2 π / La χ /dof = 0.94 R e ( R ( τ , t , q ) ) t τ =08 τ =10 τ =12 τ =14 -0.0015-0.001-0.0005 0 0.0005 0.001 0.0015 -6 -4 -2 0 2 4 6 a09m220; q =(0,1,1)2 π / La χ /dof = 1.42 R e ( R ( τ , t , q ) ) t τ =08 τ =10 τ =12 τ =14 -0.003-0.002-0.001 0 0.001 0.002 0.003 0.004 -6 -4 -2 0 2 4 6 a09m130; q =(0,1,1)2 π / La χ /dof = 1.43 R e ( R ( τ , t , q ) ) t τ =08 τ =10 τ =12 τ =14 FIG. 14: Examples of unstable two-state fits to the ratio Re( R ( τ, t, q )) defined in Eqs. (B7) with the first excited-stateenergies taken from a three-state fit to the two-point function. The data are for the three ensembles with a ≈ .
09 fm, for q = (0 , , π/La and the values of τ are specified in the labels. All data are with Θ = 0 . -0.015-0.01-0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0 0.1 0.2 0.3 0.4 0.5 F ~ , n / M N Θ – (f m ) q (GeV ) Linear [0.7] χ PT [0.4] χ PTg0 [0.8]a12m310 -0.002 0 0.002 0.004 0.006 0.008 0.01 0.012 0 0.1 0.2 0.3 0.4 0.5 F ~ , n / M N Θ – (f m ) q (GeV ) Linear [0.1] χ PT [0.1] χ PTg0 [0.1]a09m310 -0.09-0.08-0.07-0.06-0.05-0.04-0.03-0.02-0.01 0 0.01 0 0.05 0.1 0.15 0.2 0.25 0.3 F ~ , n / M N Θ – (f m ) q (GeV ) Linear [0.1] χ PT [0.1] χ PTg0 [0.0]a12m220 F ~ , n / M N Θ – (f m ) q (GeV ) Linear [0.3] χ PT [0.4] χ PTg0 [0.4]a12m220L -0.02-0.015-0.01-0.005 0 0.005 0.01 0.015 0 0.05 0.1 0.15 0.2 0.25 0.3 F ~ , n / M N Θ – (f m ) q (GeV ) Linear [0.9] χ PT [0.6] χ PTg0 [0.9]a09m220 -0.04-0.03-0.02-0.01 0 0.01 0.02 0.03 0 0.05 0.1 0.15 F ~ , n / M N Θ – (f m ) q (GeV ) Linear [1.5] χ PT [0.9] χ PTg0 [1.7]a09m130 -0.08-0.06-0.04-0.02 0 0.02 0.04 0.06 0 0.05 0.1 0.15 F ~ , n / M N Θ – (f m ) q (GeV ) Linear [0.1] χ PT [0.3] χ PTg0 [0.0]a06m135
FIG. 15: The extrapolation of ˜ F ( q ) to q → χ PT and χ PTg0,are defined in the text. The χ /dof of the fits are given within square parentheses. All data are with Θ = 0 . XII. CHIRAL-CONTINUUM EXTRAPOLATIONAND RESULTS
In this section, we present the chiral-continuum (CC)extrapolation of data for d n (and, similarly, d p ) obtainedon the seven ensembles. For each, we examine four cases. These consist of two CC fits, Linear and χ PT, using theleading order terms d n ( a, M π ) = c M π + c aM π + c a (29) d n ( a, M π ) = c M π + c L M π ln (cid:18) M π M N (cid:19) + c a , (30)4 -0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0 0.1 0.2 0.3 0.4 0.5 F ~ , p / M N Θ – (f m ) q (GeV ) Linear [0.3] χ PT [0.3] χ PTg0 [0.3]a12m310 F ~ , p / M N Θ – (f m ) q (GeV ) Linear [0.3] χ PT [0.1] χ PTg0 [0.3]a09m310 -0.04-0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.05 0.1 0.15 0.2 0.25 0.3 F ~ , p / M N Θ – (f m ) q (GeV ) Linear [0.2] χ PT [0.3] χ PTg0 [0.0]a12m220 -0.04-0.02 0 0.02 0.04 0.06 0.08 0.1 0 0.05 0.1 0.15 0.2 0.25 0.3 F ~ , p / M N Θ – (f m ) q (GeV ) Linear [1.8] χ PT [1.5] χ PTg0 [2.1]a12m220L F ~ , p / M N Θ – (f m ) q (GeV ) Linear [0.0] χ PT [0.0] χ PTg0 [0.0]a09m220 -0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.05 0.1 0.15 F ~ , p / M N Θ – (f m ) q (GeV ) Linear [2.0] χ PT [1.0] χ PTg0 [1.6]a09m130 -0.15-0.1-0.05 0 0.05 0.1 0.15 0 0.05 0.1 0.15 F ~ , p / M N Θ – (f m ) q (GeV ) Linear [0.1] χ PT [0.1] χ PTg0 [0.2]a06m135
FIG. 16: The extrapolation of ˜ F ( q ) to q → where the term c a is the O ( a ) effect discussed in Sec. XI,because of which d n,p do not vanish in the chiral limitat finite a . The ansatz are distinguished by the termsproportional to c (Linear) and c L ( χ PT). In these fits, M N is set to its physical value 940 MeV. We make thesetwo fits to the data for d n,p obtained using (i) the linearand (ii) χ PT extrapolation in q , which leads to fourestimates. These four CC fits for the neutron and theproton are shown in Figs. 17 and 18. The results and thefit coefficients c i are given in Table II.As discussed in Appendix C, at NLO in χ PT the co-efficient of the chiral logarithm c L is fixed in terms ofthe isovector scalar charge, the quark condensate and thepion decay constant, leading to ( c L ) n = − ( c L ) p = 0 . . Although the central values of the fits are ap-proximately one order of magnitude larger, our resultsare compatible with this estimate at the 1 σ –2 σ level.For the central value we take the χ PT( q ) | χ PT(CC)result and the full spread between the four for the error.The final results, using the definition in Eq. (4), are d n = − . e · fm (31) d p = 0 . e · fm (32) where the second systematic error is the spread in thefour estimates given in Table II. XIII. ANALYSIS INCLUDING THE Nπ EXCITED STATE
In this section, we describe how all ground state quan-tities change when the
N π excited state is included. Thisanalysis should be considered exploratory because (i) theextrapolations in the fits to remove ESC (see Fig. 13), (ii)the errors, and (iii) the cancellations when combining dif-ferent terms to get F using Eqs. (26) are all large.In Fig. 19, we show the increase in the value of α forthe two physical mass ensembles as compared to the datapresented in Fig. 11. The q behavior is similar to thatshown in Figs. 15 and 16, and the final results for the fourstrategies are given in Table II. The CC fits for the neu-tron and the proton using the χ PT( q ) | χ PT(CC) strat-egy are shown in Fig. 20For the central value we again take the χ PT( q ) | χ PT(CC) result and the full spread for the error. This5 -0.10-0.050.000.050.10 0 0.03 0.06 0.09 0.12 0.15 χ /dof = 0.80 F , n ( q = ) / M N Θ – (f m ) a [fm] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.10 0 0.03 0.06 0.09 0.12 χ /dof = 0.80 F , n ( q = ) / M N Θ – (f m ) M π [GeV ] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.10 0 0.03 0.06 0.09 0.12 0.15 χ /dof = 1.21 F , n ( q = ) / M N Θ – (f m ) a [fm] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.10 0 0.03 0.06 0.09 0.12 χ /dof = 1.21 F , n ( q = ) / M N Θ – (f m ) M π [GeV ] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.10 0 0.03 0.06 0.09 0.12 0.15 χ /dof = 0.78 F , n ( q = ) / M N Θ – (f m ) a [fm] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.10 0 0.03 0.06 0.09 0.12 χ /dof = 0.78 F , n ( q = ) / M N Θ – (f m ) M π [GeV ] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.10 0 0.03 0.06 0.09 0.12 0.15 χ /dof = 1.21 F , n ( q = ) / M N Θ – (f m ) a [fm] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.10 0 0.03 0.06 0.09 0.12 χ /dof = 1.21 F , n ( q = ) / M N Θ – (f m ) M π [GeV ] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap FIG. 17: The chiral-continuum extrapolation of d n using the ansatz given in Eq. (30). The four rows show (i) a linear CC fitto the data obtained using a linear extrapolation in q discussed in Sec. X; (ii) a linear CC fit to the data obtained using the χ PT extrapolation in q ; (iii) a χ PT CC fit to the data obtained using a linear extrapolation in q ; and (iv) a χ PT CC fit tothe data obtained using the χ PT extrapolation in q . All data are with Θ = 0 . -0.10-0.050.000.050.100.15 0 0.03 0.06 0.09 0.12 0.15 χ /dof = 0.46 F , p ( q = ) / M N Θ – (f m ) a [fm] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.100.15 0 0.03 0.06 0.09 0.12 χ /dof = 0.46 F , p ( q = ) / M N Θ – (f m ) M π [GeV ] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.100.15 0 0.03 0.06 0.09 0.12 0.15 χ /dof = 0.58 F , p ( q = ) / M N Θ – (f m ) a [fm] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.100.15 0 0.03 0.06 0.09 0.12 χ /dof = 0.58 F , p ( q = ) / M N Θ – (f m ) M π [GeV ] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.100.15 0 0.03 0.06 0.09 0.12 0.15 χ /dof = 0.60 F , p ( q = ) / M N Θ – (f m ) a [fm] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.100.15 0 0.03 0.06 0.09 0.12 χ /dof = 0.60 F , p ( q = ) / M N Θ – (f m ) M π [GeV ] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.100.15 0 0.03 0.06 0.09 0.12 0.15 χ /dof = 0.69 F , p ( q = ) / M N Θ – (f m ) a [fm] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.10-0.050.000.050.100.15 0 0.03 0.06 0.09 0.12 χ /dof = 0.69 F , p ( q = ) / M N Θ – (f m ) M π [GeV ] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap FIG. 18: The chiral-continuum extrapolation of d p using the ansatz given in Eq. (30). The rest is the same as in Fig. 17. NeutronFit types F / M N χ /dof c c c L c (fm) fm-GeV GeV fm-GeV Linear( q ) | Linear(CC) − − q ) | χ PT(CC) − χ PT( q ) | Linear(CC) 0.0005(17) 1.213 0.028(92) 0.8(1.2) − χ PT( q ) | χ PT(CC) − q ) | Linear(CC) 0.0076(46) 0.455 0.42(25) − q ) | χ PT(CC) 0.037(18) 0.597 − − − χ PT( q ) | Linear(CC) 0.0027(23) 0.578 0.15(13) − χ PT( q ) | χ PT(CC) 0.0238(98) 0.687 − − − Nπ excited state)Linear( q ) | Linear(CC) − − − q ) | χ PT(CC) − χ PT( q ) | Linear(CC) 0.0039(42) 2.246 0.22(23) 8.4(3.8) − χ PT( q ) | χ PT(CC) − − Nπ excited state)Linear( q ) | Linear(CC) 0.019(12) 0.347 1.04(66) − q ) | χ PT(CC) 0.140(54) 0.358 − − − χ PT( q ) | Linear(CC) 0.0040(50) 0.398 0.22(27) − χ PT( q ) | χ PT(CC) 0.068(25) 0.522 − − − d n and d p for the four fit strategies defined in the text. Also given arethe fit parameters c i defined in Eqs. (29)-(30) and the χ /dof of the fit. Results are given for two choices of the first excitedstate energy: (top) from a three-state fit to the two-point function, and (bottom) the noninteracting Nπ state. gives d n | Nπ = − . e · fm (33) d p | Nπ = 0 . e · fm (34)where the second systematic error is the spread in thefour estimates given in Table II. XIV. COMPARISON TO PREVIOUS WORK
There are two estimates [44, 66] of the contribution ofthe Θ-term to the nEDM since the clarification of theimpact of the phase α that arises in the nucleon spinorin a theory with (cid:8)(cid:8) CP in Ref. [47]. That work also containsa review of previous results, which after correction wereconsistent with zero. No estimate is given in Ref. [47], butthere is a preliminary value in a subsequent conferenceproceedings, Ref. [67]. All three of these calculationsuse the small Θ expansion and gradient flow method fortopological charge renormalization as in this work. Allresults are summarized in Table III.The ETM collaboration [66] has performed the calcula-tion on one 2+1+1-flavor twisted mass clover-improvedensemble with a = 0 . M π = 139(1) MeV, M π L = 3 .
62. Data are presented for a single value of τ = 12 so there is no information on excited state effects,continuum extrapolation, chiral behavior, or finite-sizeeffects. They also implicitly implement the Θ = 0 sub-traction (see Eqs. (26)) that we find reduces the statis- tical noise by using the spin projector (1 + γ ) iγ γ k / F (0) by making a constant fit to thelowest three q points. Their final result is taken usingthe spectral projectors method, which they find reducesthe errors by a factor of about two compared to the field-theoretic definition of the topological charge used in thiswork. They do not, however, assess a systematic errorassociated with excited-state effects, extrapolation in q ,or the chiral-continuum fit.The calculation presented in Ref. [44] uses six 2+1-flavor Wilson-Clover ensembles but only one below M π =567 MeV, with M π = 410 MeV. The values of latticespacings range between 0 . < a < .
11 fm. A linearfit in q is made to obtain F (0). Also, artifacts due toESC are not analyzed and, in any case, data with theheavy pion masses studied, M π >
410 MeV, would notprovide sensitivity to analyses with or without includinga
N π state. This is the only other calculation that haspresented a chiral extrapolation using the χ PT ansatz(Eq. (30) but with a O ( a ) discretization correction in-stead of our c a term). As shown in the bottom rightpanels in Figs. 17 and 18, such chiral fits have an inflec-tion point close to the smallest M π data point in orderto satisfy the constraint F = 0 at M π = 0. In the caseof Ref. [44], this occurs around M π = 400 MeV, raisingquestions on the reliability of the extrapolation.8 a09m130 χ /dof = 0.85 r α / Θ – τ p = 0.00 GeV p = 0.10 GeV a06m135 χ /dof = 0.11 r α / Θ – τ p = 0.00 GeV p = 0.10 GeV FIG. 19: The increase in the value of α when fits to the two-point functions are made including a Nπ excited state as comparedto data in Fig. 11. -0.4-0.3-0.2-0.10.00.10.20.30.4 0 0.03 0.06 0.09 0.12 0.15 χ /dof = 2.43 F , n ( q = ) / M N Θ – (f m ) a [fm] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.4-0.3-0.2-0.10.00.10.20.30.4 0 0.03 0.06 0.09 0.12 χ /dof = 2.43 F , n ( q = ) / M N Θ – (f m ) M π [GeV ] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.4-0.20.00.20.40.60.8 0 0.03 0.06 0.09 0.12 0.15 χ /dof = 0.52 F , p ( q = ) / M N Θ – (f m ) a [fm] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap -0.4-0.20.00.20.40.60.8 0 0.03 0.06 0.09 0.12 χ /dof = 0.52 F , p ( q = ) / M N Θ – (f m ) M π [GeV ] a12m310a12m220a12m220L a09m310a09m220a09m130 a06m135Extrap FIG. 20: The chiral-continuum extrapolation of d n (top) and d p (bottom) using the ansatz given in Eq. (30), using Nπ as theexcited state fits, and with the χ PT( q ) | χ PT(CC) strategy. All data are with Θ = 0 . · fm Θ e · fmThis Work d n = − . d p = 0 . Nπ d n = − . d p = 0 . | d n | = 0 . d n = − . d p = 0 . d n ≈ .
001 –TABLE III: Summary of lattice results for the contribution of the Θ-term to the neutron and proton electric dipole moment. XV. CONCLUSIONS
This paper presents a calculation of the contribution ofthe Θ-term to the nucleon electric dipole moment using2+1+1-flavor HISQ ensembles and Wilson-clover valencequarks. Two of the seven ensembles analyzed are at thephysical pion mass, which anchor our chiral fits. Thecalculation has been done using the small Θ expansionmethod. Significant effort has been devoted to getting areliable signal in the (cid:8)(cid:8)
CP violating form factor F . Thegradient flow scheme has been used to renormalize the Θ-term and the results are shown to be independent of theflow time. Our estimate of the topological susceptibilityfor the 2+1+1 theory is χ Q = (66(9)(4) MeV) in thecontinuum limit at M π = 135 MeV.We also present two technical issues. First, in Ap-pendix D, we show that, in chiral perturbation theory,the N π excited state should provide the dominant con-tamination. We have, therefore, used two strategies forremoving excited state contamination. In the first, themass gaps are taken from fits to the spectral decomposi-tion of the nucleon two-point function, and in the secondwe assume they are given by the non-interacting energyof the N ( ) π ( ) state. We find a very significant differ-ence between the two as shown in Secs. IX and XIII, andby the results summarized in Tables II and III.The second technical issue discussed in Sec. XI andappendix E is that lattice artifacts introduce a term pro-portional to am q , because of which d n does not vanishin the chiral limit at finite a . Our chiral-continuum fitshave been made including this term.The analysis of the q dependence of F Θ3 has been car-ried out using both a linear and the leading order χ PTexpression as described in Sec. X. The current data donot distinguish between the two. Similarly, the chiralfit is also carried out using a linear and the leading order χ PT expression as described in Sec. XII. The results fromthese four sets of fits and the two strategies to removeexcited-state contributions are summarized in Table II.Our preferred values are obtained using the leadingorder χ PT expressions. The analysis using excited statesfrom fits to the two-point function indicate that d Θ n issmall, | d Θ n | (cid:46) .
01 Θ e · fm, whereas for the proton we get | d Θ p | ∼ .
02 Θ e · fm. On the other hand, if the dominantexcited-state contribution is from the N π state, then | d Θ n | could be as large as 0 .
05 Θ e · fm and | d Θ p | ∼ .
07 Θ e · fm.Lastly, we find the sign of d Θ p to be opposite to that of d Θ n .From the final summary of results presented in Ta-ble III, which also includes estimates from previousworks, it is clear that, at present, lattice calculations donot yet provide a reliable estimate. To improve the cur-rent 100% uncertainty to a 3 σ result will require a factorof at least ten improvement in statistics. Acknowledgments
The calculations used the Chroma software suite [68].This research used resources at (i) the National EnergyResearch Scientific Computing Center, a DOE Office ofScience User Facility supported by the Office of Sci-ence of the U.S. Department of Energy under ContractNo. DE-AC02-05CH11231; (ii) the Oak Ridge Leader-ship Computing Facility, which is a DOE Office of Sci-ence User Facility supported under Contract DE-AC05-00OR22725; (iii) the USQCD Collaboration, which arefunded by the Office of Science of the U.S. Department ofEnergy, and (iv) Institutional Computing at Los AlamosNational Laboratory. T. Bhattacharya and R. Guptawere partly supported by the U.S. Department of En-ergy, Office of Science, Office of High Energy Physicsunder Contract No. DE-AC52-06NA25396. We acknowl-edge support from the U.S. Department of Energy, Officeof Science, Office of Advanced Scientific Computing Re-search and Office of Nuclear Physics, Scientific Discoverythrough Advanced Computing (SciDAC) program, andof the U.S. Department of Energy Exascale ComputingProject. T. Bhattacharya, V. Cirigliano, R. Gupta, E.Mereghetti and B.Yoon were partly supported by theLANL LDRD program.
Appendix A: Connection between Minkowski andEuclidean notations
To make our conventions explicit, we present the con-nection between Minkowski and Euclidean variables inTable IV.To connect the Lagrangian density for the Θ term inMinkowski and Euclidean spaces, we take the Minkowskiaction associated with the QCD Θ term to be S M Θ = − Θ32 π (cid:90) d x ( G M ) a µν ( x ) ( ˜ G M ) aµν ( x ) (A1)where ( ˜ G M ) aµν = (1 / (cid:15) µναβM ( G M ) aαβ and ( (cid:15) M ) =+1 = − (cid:15) M . Upon rotating to the Euclidean space onegets d x M = − id x E and (cid:15) µναβM ( G M ) aµν ( G M ) aαβ = i ( (cid:15) E ) µναβ ( G E ) aµν ( G E ) aαβ . (A2)The factor of + i arises from the transformation of thefield strength and because each term in the sum has onefactor of G i (or G i ) and one factor of G jk . Moreover,we used (cid:15) ijkM ≡ ( (cid:15) E ) ijk = ( (cid:15) E ) ijk , (A3)which implies ( (cid:15) E ) ijk = − ( (cid:15) M ) ijk (A4)and hence ( (cid:15) E ) = +1.0 Quantity Minkowsky ↔ Euclidean Remarks4-vector v µ v M = v M = − iv E = − iv E Ensures v M · v (cid:48) M = − v E · v (cid:48) E ; v iM = − v Mi = v iE = iv Ei In particular, v M = − v E . t ≡ x M = − ix E ≡ − iτp M = − ip E = E Derivatives ∂ M = i∂ E ∂ µ = ∂/∂x µ and ∂ µ = ∂/∂x µ in both E and M ∂ Mi = − ∂ Mi = − ∂ Ei = − ∂ Ei Gauge Fields A M = iA E D µ = ∂ µ − A µ transforms homogeneously A Mi = − A Mi = − A Ei = − A Ei ( G M ) i = − i ( G E ) i , ( G M ) i = i ( G E ) i ( G M ) ij = ( G E ) ij , ( G M ) ij = ( G E ) ij γ matrices γ E = γ M , γ iE = − iγ iM We adopt the DeGrand-Rossi basis [69]. These γ E = γ E γ E γ E γ E = − iγ M γ M γ M γ M = − γ M = − γ c Euclidean gamma matrices are Hermitean. γ µM ≡ γ c γ c γ µc γ c γ c Minkowski gamma matrices are unitarily transformedfrom the standard chiral basis, γ µc [70] /p M = − i/p E /D M = i /D E ψ M = γ c γ c ψ c and ¯ ψ M = ¯ ψ c γ c γ c .Charge C M = iγ M γ M Conjugation C c = iγ c γ c Matrix C E = γ E γ E TABLE IV: Connection between Euclidean and Minkowsky variables.
Putting together the change in the measure and thechange in the Lagrangian density we have S M Θ = − Θ64 π (cid:15) µναβE (cid:90) d x E ( G E ) aµν ( x ) ( G E ) aαβ ( x ) , (A5)and hence ( iS M = − S E ) S E Θ = + i Θ64 π (cid:15) µναβE (cid:90) d x E ( G E ) aµν ( x ) ( G E ) aαβ ( x ) , (A6)consistently with Eq. (1). Appendix B: Extraction of F The Euclidean four-vector V µ ( q ) defined in Eq. (25)can be determined from lattice data by taking appropri-ate ratios of 3-pt function and 2-pt functions. This isachieved by defining the projected 2- and 3-point func-tions as follows, C pt ( t, p ) = Tr (cid:2) P pt (cid:104) Ω | N ( p , t ) ¯ N ( p , | Ω (cid:105) (cid:3) (B1) C µ pt ( τ, t, q ) = Tr (cid:2) P pt (cid:104) Ω | N ( p (cid:48) , τ ) J EM µ ( t ) ¯ N ( p , | Ω (cid:105) (cid:3) , (B2)with q = p (cid:48) − p , p (cid:48) = 0, P pt given in Eq. (24), P pt = 12 (1 + γ ) , (B3)and, neglecting the contributions of heavier quarks, J EM µ = e (cid:16) (2 / uγ µ u − (1 /
3) ¯ dγ µ d − (1 / sγ µ s (cid:17) . (B4) The ratio˜ R µ ≡ C µ pt ( τ, t, q ) C pt ( τ, p (cid:48) ) × (cid:18) C pt ( t, p (cid:48) ) C pt ( τ, p (cid:48) ) C pt ( τ − t, p ) C pt ( t, p ) C pt ( τ, p ) C pt ( τ − t, p (cid:48) ) (cid:19) / (B5)becomes independent of t and τ if t, τ are sufficiently largethat excited state effects can be neglected, and takes theform V µ ( q ) (cid:112) E p E p (cid:48) ( E p + M N cos(2 α N ))( E p (cid:48) + M N cos(2 α N )) . (B6)In our plots to demonstrate the signal and excited states,we, therefore, choose to show the quantity R µ ( τ, t, q ) ≡ ˜ R µ g V × (cid:113) E p E p (cid:48) ( E p + M N cos(2 α N ))( E p (cid:48) + M N cos(2 α N )) , (B7)where g V ≡ C µ pt ( τ, t, ) / C pt ( t, ), and α N is calculatedfrom fits to the 2-pt functions with momentum p or p (cid:48) asdiscussed in Section VII.The components of V µ are expressed in terms of formfactors F , , ,A ( q ) defined in Eq. (6) as follows:1 V = ic α N M N ( q + iq ) F ( q )+ (cid:26) − c α N M N q −
12 [ s α N q q + ic α N q ( E N − M N )] (cid:27) F ( q ) − i [ s α N q ( E N − m N ) − c α N q q ] F A ( q ) −
12 [ c α N q q − is α N q ( E N − M N )] F ( q ) , (B8a) V = c α N M N ( q + iq ) F ( q )+ (cid:26) c α N M N q −
12 [ s α N q q + ic α N q ( E N − M N )] (cid:27) F ( q )+ 2 i [ s α N q ( E N − m N ) + c α N q q ] F A ( q ) −
12 [ c α N q q − is α N q ( E N − M N )] F ( q ) , (B8b) V = M N [ ic α N q + s α N ( E N − M N )] F ( q )+ 12 (cid:8) − ic α N ( E N − M N ) q − s α N q + 2 s α N M N ( E N − M N )] (cid:9) F ( q ) − ic α N (cid:2) q + q (cid:3) F A ( q ) − (cid:2) c α N q − is α N q ( E N − M N ) (cid:3) F ( q ) , (B8c) V = M N [ c α N ( E N + M N ) − is α N q ] F ( q ) − (cid:8) c α N ( E N − M N ) − is α N q ( E N − M N )] (cid:9) F ( q )+ 12 (cid:2) ic α N q ( E N + M N ) + s α N ( E N − M N ) (cid:3) F ( q ) , (B8d)where c α ≡ (cos 2 Re α + cosh 2 Im α ) / s α ≡ (sin 2 Re α + i sinh 2 Im α ) /
2. For PT symmetric theories,where α is real, these expressions simplify to c α = cos α and s α = cos α sin α .From the above expressions we want to extract F ( q ),that gives the neutron EDM. It turns out that the RHS ofEqs. (B8) is most naturally expressed in terms of G , , given by G = F + F (B9a) G = F − q E m F + sc q E m F , (B9b) G = F + sc F (B9c)where q E = q + q and s ≡ sin α cos α , c ≡ cos α .For a given momentum transfer q = ( q , q , q ),Eqs. (B8) thus represents eight equations for G , , . Theycan be written in a compact form as follows: K ( q ) G G G − V ( q ) = 0 , (B10)where K ( q ) is an 8 × K ( q ) = (cid:18) X ( q ) 0 X ( q )0 Y ( q ) 0 (cid:19) (B11a) X ( q ) = m − cq cq s ( E − m ) − isq (B11b) X ( q ) = − c q q q q − i ( E + m ) (B11c) Y ( q ) = m c q q q − i ( E + m ) , (B11d)and V ( q ) is an eight-dimensional array given by V ( q ) = (cid:18) V R ( q ) V I ( q ) (cid:19) , (B12a) V R ( q ) = (cid:18) Re (cid:126) V ( q ) i Im V ( q ) (cid:19) , (B12b) V I ( q ) = (cid:18) Im (cid:126) V ( q ) − i Re V ( q ) (cid:19) . (B12c)To solve for G , , ( q ), for a given three-momentumtransfer q = ( q , q , q ) we can use a least squares esti-2mator. Namely, we minimize the function F ( G , , ) = (cid:88) (cid:126)q ∈ P ( (cid:126)q ) 8 (cid:88) i,j =1 w ij ( q ) E i ( q ) E j ( q ) (B13)where E i ( q ) = (cid:88) β =1 K iβ ( q ) G β − V i ( q ) (B14) w ij ( q ) = (cid:2) C − V ( q ) (cid:3) ij (B15)where the weights matrix is the inverse of the covariancematrix of lattice “measurements” V i ( q ):[ C V ( q )] ij = Cov ( V i ( q ) , V j ( q )) . (B16)For independent variables V i ( q ) the covariance matrix C V and it inverse are positive definite. This guarantees that F ( G , , ) is minimized if and only if E i ( q ) = 0 for all i .The sum over momenta runs over the six permutations( q , q , q ), ( q , q , q ), ( q , q , q ), ( q , q , q ), ( q , q , q ),( q , q , q ).The function F ( G , , ) is stationary for ∂F∂G α = 0 α = 1 , , . (B17)Explicitly, since ∂E j /∂G α = K jα , one finds2 (cid:88) (cid:126)q ∈ P ( (cid:126)q ) 8 (cid:88) i,j =1 w ji ( q ) E i ( q ) K jα ( q ) = 0 , α = 1 , , . (B18)or even more explicitly (cid:88) (cid:126)q ∈ P ( (cid:126)q ) 8 (cid:88) i,j =1 w ji ( q ) (cid:88) β K iβ ( q ) G β − V i ( q ) K jα ( q )= 0 , α = 1 , , , (B19)which is a system of three equations for G , , . The ex-tremum condition for F ( G , , ) implies the following lin-ear equation for G , , ( q ): A αβ G β = B α (B20)where the 3 × A and the three dimensional array B are given by A αβ = (cid:88) (cid:126)q ∈ P ( (cid:126)q ) 8 (cid:88) i,j =1 K jα ( q ) w ji ( q ) K iβ ( q ) (B21a) For ease of notation, we are ignoring current conservation, whichrelates the various components V i ( q ). Strictly speaking, we needto eliminate the dependent components of V i ( q ) when using aconserved current to get an invertible covariance matrix. B α = (cid:88) (cid:126)q ∈ P ( (cid:126)q ) 8 (cid:88) i,j =1 K jα ( q ) w ji ( q ) V i ( q ) . (B21b)So from the lattice data on V i ( q ), their covariance ma-trix, and the explicit form of the matrix K iα ( q ) given inEq. (B11) one can construct A αβ and B α and solve for G , , . Error on G , , can be assigned with the boot-strap method. Appendix C: Chiral extrapolation formulae
We can express the electric dipole form factor as F i ( q )2 M N = d i − S (cid:48) i q + H i ( q ) , (C1)where d i is the EDM, S (cid:48) i the Schiff moment (with someabuse of notation), and H i ( q ) account for the higherorder dependence on q . Here, i is an isospin label, andthe results are more conveniently expressed in terms ofan isoscalar ( i = 0) and isovector ( i = 1) component.The neutron and proton form factors are F , p ( q ) = F ( q ) + F ( q ) ,F , n ( q ) = F ( q ) − F ( q ) . (C2)At NLO in χ PT, the EDMs are given by [24, 28, 30, 32], d = e ¯ d + eg A ¯ g (4 πF π ) (cid:20) πM π M N (cid:21) , (C3) d = e ¯ d ( µ ) + eg A ¯ g (4 πF π ) (cid:20) − ln M π µ + 5 π M π M N (cid:21) , (C4)where the renormalization scale dependence of the LEC¯ d cancels the µ in the logarithm. Here g A = 1 . F π =92 . g is a (cid:8)(cid:8) CP pion-nucleon coupling, defined as L = − ¯ g F π ¯ N π · τ N, (C5)which is related by chiral symmetry to the neutron-proton mass splitting [24]¯ g = (cid:18) M n − M p ¯ mε + O (cid:18) M π Λ χ (cid:19)(cid:19) m ∗ ¯Θ = g S ¯ m ¯Θ , (C6)where m − ∗ = m − u + m − d , 2 ¯ m = m u + m d , and Λ χ ∼ χ PT expansion breaksdown. g S is the isovector scalar charge, and the lastequality holds in the isospin limit. At the physical pionmass, one obtains [71]¯ g F π = (15 . ± . · − ¯Θ , (C7)but the last term in Eq. (C6) allows to extend the relationto arbitrary masses in the regime of validity of χ PT. Inparticular, in the χ PT fits to F ( q ) we use¯ g = g S B M π ¯Θ , (C8)3with g S = 1 . B = 2 . d , are two low-energyconstants, which, by naive-dimensional-analysis, scale as¯ d , = O (cid:18) M π (4 πF π ) (cid:19) (C9)The first derivative of the form factor is [28, 30, 32] S (cid:48) = 0 , (C10) S (cid:48) = eg A ¯ g πF π ) M π (cid:20) − π M π M N (cid:21) . (C11)At N LO there are additional long- and short-distancecontributions to both isoscalar and isovector components.The remaining momentum dependence of the EDFF isgiven by the functions H i ( q ) introduced in Eq. (C1), H ( q ) = 0 , (C12) H ( q ) = 4 eg A ¯ g πF π ) (cid:20) h a ( x ) − π M π M N h b ( x ) (cid:21) , (C13)with x ≡ q / M π . h a appears at leading order, h a ( x ) = − (cid:34)(cid:114) x ln (cid:32) (cid:112) /x + 1 (cid:112) /x − (cid:33) − (cid:16) x (cid:17)(cid:105) , (C14)while h b is generated at NLO h b ( x ) = − (cid:34) x ) (cid:18) (cid:18) √ x arctan √ x − x (cid:19)(cid:19) − x (cid:35) . (C15)Since these behave as h ( n ) i ( x ) = x + O ( x ) for x (cid:28) O ( q ), dependence of H i is consistent withthe definition in Eq. (C1). Appendix D: Excited state contamination in chiralperturbation theory
In this appendix, we show that, in χ PT, the gap be-tween the ground state and excited state contributions tothe CP-odd components of the three-point function C µ pt is expected to be of order of the pion mass M π . This canbe intuitively understood from the fact that the nucleonEDM induced by the QCD ¯Θ term receives a LO contri-bution from a long-range pion loop [24], shown in Figure21. In Minkowski space, this diagram has a branch cutwhen the intermediate pions and nucleon go on-shell. InEuclidean space, this translates into a N π excited state,whose amplitude is of the same size as the ground statecontribution. For simplicity, we focus only on the dia-gram shown in Figure 21, and assume that the nucleoninterpolating field does not couple to nucleon plus pions. p -k ,-k-q p -k ,-kp ,-q k ,k p ,0 FIG. 21: Leading order diagram for the excited states contri-bution to the three-point function C µ pt in chiral perturbationtheory. A black square denotes an insertion of the CP-oddpion-nucleon couplings ¯ g . Filled circles denote CP-even pion-nucleon and pion-photon couplings. We start from the 4 th component of the three-pointfunction. Carrying out the Dirac traces in Eq. (25), inthe limit M N (cid:29) q , we find C pt = q τ ¯ g g A (4 πF π ) e − M N t B − E N t (cid:40) f ( M π , q, L )+ (4 π ) L M π E π (cid:32) e − M π t + e − M π t B + E π M π (cid:0) e − E π t + e − E π t B (cid:1) − M π + E π M π (cid:0) e − E π t − M π t B + se − M π t − E π t B (cid:1) + ( E π − m π ) M π ( E π + M π ) × (cid:16) e − ( M π + E π ) t + e − ( M π + E π ) t B (cid:17) (cid:33) + . . . (cid:41) , (D1)where t B = τ − t , E N = (cid:112) M N + q ∼ M N , E π = (cid:112) M π + q and . . . denotes terms with a gap with twoor more units of momentum. f ( M π , q, L ) denotes theground state loop function, which we write as an infinitevolume term f ∞ and a correction ∆ f ( M π , q, L ) = f ∞ ( M π , q ) + ∆( M π , q, L ) . (D2)In the non-relativistic limit, f ∞ ( M π , q ) is given by f ∞ ( M π , q ) = (4 π ) (cid:32) (cid:90) d k (2 π ) k + (cid:126)k + M π k + ( (cid:126)k + (cid:126)q ) + M π (cid:33) , (D3)4and is ultraviolet divergent. In dimensional regulariza-tion and in the MS scheme f ∞ ( M π , q ) = log µ M π + 2 − (cid:114) x ln (cid:113) x + 1 (cid:113) x − , (D4)with x = q / (4 M π ), which is of course the same functionas in Section C. The finite volume correction is given by∆( M π , q, L ) = (4 π ) (cid:90) dk π L (cid:88) (cid:126)k − (cid:90) d k (2 π ) k + (cid:126)k + M π k + ( (cid:126)k + (cid:126)q ) + M π , (D5)which can be written in terms of Bessel functions as [72]∆( M π , q, L )=2 (cid:88) (cid:126)n (cid:54) =0 (cid:90) dxK (cid:16) L (cid:112) M π + q x (1 − x ) | (cid:126)n | (cid:17) . (D6)At q = 0, for M π L ∼
4, ∆ amounts to a 0 .
1% correc-tion. Eqs. (D1) and (D4) thus show that the excitedstates have a gap of O ( M π ). The ratio of the ground andexcited state contributions is determined by the quan-tity (4 π ) / ( LM π ) , which is a number of order 1 for LM π = 4. We thus do not expect a significant sup-pression of the excited states. A similar calculation canbe performed for the spatial components C i pt , yielding aresult similar to Eq. (D1), but with a sinh rather thancosh behavior. Appendix E: O ( a ) corrections in the Wilson-Clovertheory In this appendix, we analyze CP violation due to thetopological charge in the Wilson-Clover theory at O ( a ).We will denote by O ( d ) n , ˜ O ( d ) n , O ( d ) , ren n , the set of bare,subtracted, and renormalized operators of dimension d ,respectively. Subtracted operators, i.e., operators free ofpower divergences, are defined by˜ O ( d ) n (cid:48) = O ( d ) n (cid:48) − (cid:88) d (cid:48) 13 Tr [ ˆ m ] ∂ µ (cid:0) ¯ ψγ µ γ ψ (cid:1) (E20) O (5)8 = ¯ ψiγ ˆ m ψ (E21) O (5)9 = Tr (cid:2) ˆ m (cid:3) ¯ ψiγ ψ (E22) O (5)10 = Tr [ ˆ m ] ¯ ψiγ ˆ mψ (E23) O (5)11 ≡ P EE = i ¯ ψ E γ ψ E (E24) O (5)12 ≡ ∂ · A E = ∂ µ [ ¯ ψ E γ µ γ ψ + ¯ ψγ µ γ ψ E ] (E25) O (5)13 ≡ A ∂ = ¯ ψγ /∂ψ E − ¯ ψ E ←− /∂ γ ψ (E26) O (5)14 ≡ A A ( γ ) = ie (cid:16) ¯ ψQ /A ( γ ) γ ψ E − ¯ ψ E Q /A ( γ ) γ ψ (cid:17) , (E27)where ˜ σ µν ≡ ( σ µν γ + γ σ µν ) and ψ E = ( /D + ˆ m ) ψ .Keeping in mind that O ( x , ..., x n ) has the structure N ( x ) J EM µ ( x ) ¯ N ( x ), in terms of the neutron source andsink operator and the electromagnetic current, the vari-ous O (5) n contribute to Eq. (E12) as follows: • O (5)1 is the isoscalar chromo-EDM operator andcontributes an O ( a ) term to the LHS of Eq. (E12).In fact, as shown below, this is the leading O ( a )contribution, thus proving a linear relation betweenisovector insertions of the pseudoscalar density andthe chromo-EDM. • O (5)2 , , are total derivatives and their insertion inEq. (E12) vanish upon integration over (cid:82) d x . • O (5)3 , involve one and two powers of the electromag-netic field strength. In order to eliminate the pho-ton field in the correlation functions in Eq. (E12),one needs electromagnetic loops, making the con-tribution of O (5)3 , to Eq. (E12) of O ( a α EM /π ), andthus negligible to the order we are working. • O (5)5 provides a correction of O ( am ) proportionalto ( G ˜ G ) in the LHS of Eq. (E12). • O (5)8 , , become ˆ m ¯ ψiγ ψ when ˆ m ∝ I . Therefore,their contributions have the same form of the pseu-doscalar insertion in Eq. (E12), but suppressed by O ( am ). • The operators O (5)11 , , , vanish by using the quarkequations of motion and can contribute contactterms to the LHS of Eq. (E12). However, it turnsout that none of them actually contributes at thisorder. O (5)11 contains two equation of motion op-erators. Therefore, when inserted in Eq. (E12), itwill always involve a contraction with a quark field6in the neutron source or sink operator, and thusit will not contribute to the residue of the neu-tron pole. O (5)12 is a total derivative and drops outof Eq. (E12). O (5)13 is gauge-variant operator anddrops out of Eq. (E12) as long as O ( x , ..., x n ) is agauge singlet, which is the case for O ( x , x , x ) ∝ N ( x ) J EM µ ( x ) ¯ N ( x ). O (5)14 involves the photonfield and therefore can contribute to Eq. (E12) onlyto O ( aα EM /π ).So in summary, for ˆ m ∝ I , Eq. (E12) becomes (cid:90) d x (cid:28) O ( x , ..., x n ) (cid:16) − m ¯ ψ ( x ) γ ψ ( x ) (cid:16) O ( am ) (cid:17) − N F π ( G ˜ G ) ren (cid:16) O ( am ) (cid:17) − aK X ˜ O (5)1 (cid:17)(cid:29) = − (cid:90) d x (cid:28) δO ( x , ..., x n ) δ ( iα ( x )) (cid:29) . (E28)If ˆ m (cid:54) = I , the singlet AWI, Eq. (E12), involves ¯ ψ ˆ mγ ψ .All the arguments above go through, except for the effectof O (5)8 , , . O (5)10 gives a correction of O ( am ) proportionalto ¯ ψ ˆ mγ ψ , while O (5)8 , contribute nonmultiplicative terms involving the nonsinglet pseudoscalar densities of O ( a ˆ m )in Eq. (E28). The presence of these additional terms doesnot affect our conclusion about the existence of O ( am q )corrections. [1] M. Pospelov and A. Ritz, Annals Phys. , 119 (2005),arXiv:hep-ph/0504231.[2] J. Engel, M. J. Ramsey-Musolf, and U. van Kolck, Prog.Part. Nucl. Phys. , 21 (2013), arXiv:1303.2371 [nucl-th].[3] C. Bennett et al. (WMAP Collaboration), Astrophys. J.Suppl. , 1 (2003), arXiv:astro-ph/0302207.[4] E. W. Kolb and M. S. Turner, Front. Phys. , 1 (1990).[5] P. Coppi, eConf C040802 , L017 (2004).[6] A. Sakharov, Pisma Zh. Eksp. Teor. Fiz. , 32 (1967).[7] M. Kobayashi and T. Maskawa, Prog. Theor. Phys. ,652 (1973).[8] Z. Maki, M. Nakagawa, and S. Sakata, Prog. Theor. Phys. , 870 (1962).[9] H. Nunokawa, S. J. Parke, and J. W. Valle, Prog. Part.Nucl. Phys. , 338 (2008), arXiv:0710.0554 [hep-ph].[10] I. Khriplovich and A. Zhitnitsky, Phys. Lett. B , 490(1982).[11] A. Czarnecki and B. Krause, Phys. Rev. Lett. , 4339(1997), arXiv:hep-ph/9704355.[12] C.-Y. Seng, Phys. Rev. C , 025502 (2015),arXiv:1411.1476 [hep-ph].[13] C. Abel et al. (nEDM), Phys. Rev. Lett. , 081803(2020), arXiv:2001.11966 [hep-ex].[14] B. Graner, Y. Chen, E. Lindahl, and B. Heckel,Phys. Rev. Lett. , 161601 (2016), [Erratum: Phys.Rev. Lett. 119, 119901 (2017)], arXiv:1601.04339[physics.atom-ph].[15] T. Ito, SNS nEDM progress and plans, https://indico.frib.msu.edu/event/13/contributions/191/attachments/79/384/SNSnEDM_Ito.pdf (2019).[16] R. Jackiw and C. Rebbi, Phys. Rev. Lett. , 172 (1976).[17] J. Callan, Curtis G., R. Dashen, and D. J. Gross, Phys.Lett. B , 334 (1976).[18] D. J. Gross, R. D. Pisarski, and L. G. Yaffe, Rev. Mod.Phys. , 43 (1981).[19] A. Dolgov, Phys. Rept. , 309 (1992). [20] V. Kuzmin, M. Shaposhnikov, and I. Tkachev, Phys. Rev.D , 466 (1992).[21] I. I. Bigi and N. Uraltsev, Nucl. Phys. B , 321 (1991).[22] M. Pospelov and A. Ritz, Phys. Rev. Lett. , 2526(1999), arXiv:hep-ph/9904483.[23] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. , 1440(1977).[24] R. Crewther, P. Di Vecchia, G. Veneziano, and E. Witten,Phys. Lett. B , 123 (1979), [Erratum: Phys.Lett.B 91,487 (1980)].[25] A. Pich and E. de Rafael, Nucl. Phys. B , 313 (1991).[26] P. L. Cho, Phys. Rev. D , 3304 (1993), arXiv:hep-ph/9212274.[27] B. Borasoy, Phys. Rev. D , 114017 (2000), arXiv:hep-ph/0004011.[28] W. Hockings and U. van Kolck, Phys. Lett. B , 273(2005), arXiv:nucl-th/0508012.[29] S. Narison, Phys. Lett. B , 455 (2008),arXiv:0806.2618 [hep-ph].[30] K. Ottnad, B. Kubis, U.-G. Meissner, and F.-K. Guo,Phys. Lett. B , 42 (2010), arXiv:0911.3981 [hep-ph].[31] J. de Vries, R. Timmermans, E. Mereghetti, and U. vanKolck, Phys. Lett. B , 268 (2011), arXiv:1006.2304[hep-ph].[32] E. Mereghetti, J. de Vries, W. Hockings, C. Maekawa,and U. van Kolck, Phys. Lett. B , 97 (2011),arXiv:1010.4078 [hep-ph].[33] M. Pospelov and A. Ritz, Phys. Rev. D , 073015(2001), arXiv:hep-ph/0010037.[34] O. Lebedev, K. A. Olive, M. Pospelov, and A. Ritz, Phys.Rev. D , 016003 (2004), arXiv:hep-ph/0402023.[35] K. Fuyuto, J. Hisano, and N. Nagata, Phys. Rev. D ,054018 (2013), arXiv:1211.5228 [hep-ph].[36] U. Haisch and A. Hala, JHEP (11), 154,arXiv:1909.08955 [hep-ph].[37] E. Shintani, S. Aoki, N. Ishizuka, K. Kanaya,Y. Kikukawa, Y. Kuramashi, M. Okawa, Y. Tanigchi, A. Ukawa, and T. Yoshie, Phys. Rev. D , 014504(2005), arXiv:hep-lat/0505022.[38] F. Berruto, T. Blum, K. Orginos, and A. Soni, Phys. Rev.D , 054509 (2006), arXiv:hep-lat/0512004.[39] A. Shindler, J. de Vries, and T. Luu, PoS LAT-TICE2014 , 251 (2014), arXiv:1409.2735 [hep-lat].[40] F.-K. Guo, R. Horsley, U.-G. Meissner, Y. Naka-mura, H. Perlt, P. Rakow, G. Schierholz, A. Schiller,and J. Zanotti, Phys. Rev. Lett. , 062001 (2015),arXiv:1502.02295 [hep-lat].[41] A. Shindler, T. Luu, and J. de Vries, Phys. Rev. D ,094518 (2015), arXiv:1507.02343 [hep-lat].[42] C. Alexandrou, A. Athenodorou, M. Constantinou,K. Hadjiyiannakou, K. Jansen, G. Koutsou, K. Ottnad,and M. Petschlies, Phys. Rev. D , 074503 (2016),arXiv:1510.05823 [hep-lat].[43] E. Shintani, T. Blum, T. Izubuchi, and A. Soni, Phys.Rev. D , 094503 (2016), arXiv:1512.00566 [hep-lat].[44] J. Dragos, T. Luu, A. Shindler, J. de Vries, andA. Yousif, Confirming the Existence of the strong CPProblem in Lattice QCD with the Gradient Flow (2019),arXiv:1902.03254 [hep-lat].[45] T. Bhattacharya, V. Cirigliano, R. Gupta, H.-W. Lin,and B. Yoon, Phys. Rev. Lett. , 212002 (2015),arXiv:1506.04196 [hep-lat].[46] R. Gupta, B. Yoon, T. Bhattacharya, V. Cirigliano, Y.-C. Jang, and H.-W. Lin, Phys. Rev. D98 , 091501 (2018),arXiv:1808.07597 [hep-lat].[47] M. Abramczyk, S. Aoki, T. Blum, T. Izubuchi, H. Ohki,and S. Syritsyn, Phys. Rev. D96 , 014501 (2017),arXiv:1701.07792 [hep-lat].[48] T. Bhattacharya, B. Yoon, R. Gupta, and V. Cirigliano,in Proceedings of LATTICE2018: The 36th Annual In-ternational Symposium on Lattice Field Theory (2018)arXiv:1812.06233 [hep-lat].[49] J. Kim, J. Dragos, A. Shindler, T. Luu, and J. de Vries,PoS LATTICE2018 , 260 (2019), arXiv:1810.10301[hep-lat].[50] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finken-rath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, andA. Vaquero Aviles-Casco, Phys. Rev. D , 014509(2019), arXiv:1812.10311 [hep-lat].[51] A. Bazavov et al. (MILC Collaboration), Phys. Rev. D87 , 054505 (2013), arXiv:1212.4768 [hep-lat].[52] R. Gupta, Y.-C. Jang, B. Yoon, H.-W. Lin, V. Cirigliano,and T. Bhattacharya, Phys. Rev.