Contribution of the QCD $\Theta$-term to nucleon electric dipole moment

We present a calculation of the contribution of the $\Theta$-term to the neutron and proton electric dipole moments using seven 2+1+1-flavor HISQ ensembles. We also estimate the topological susceptibility for the 2+1+1 theory to be $\chi_Q = (66(9)(4) \rm MeV)^4$ in the continuum limit at $M_\pi = 135$ MeV. The calculation of the nucleon three-point function is done using Wilson-clover valence quarks. The CP-violating form factor $F_3$ is calculated by expanding in small $\Theta$. We show that lattice artifacts introduce a term 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(0) \pi(0)$ state should provide the leading excited state contribution, and we study the effect of such a state. Detailed analysis of the contributions to the neutron and proton electric dipole moment using two strategies for removing excited state contamination are presented. Using the excited state spectrum from fits to the two-point function, we find $d_n^\Theta$ is small, $|d_n^\Theta| \lesssim 0.01 \overline \Theta e$ fm, whereas for the proton we get $|d_p^\Theta| \sim 0.02 \overline \Theta e$ fm. On the other hand, if the dominant excited-state contribution is from the $N \pi$ state, then $|d_n^\Theta|$ could be as large as $0.05 \overline \Theta e$ fm and $|d_p^\Theta| \sim 0.07 \overline \Theta e$ fm. Our overall conclusion is that present lattice QCD calculations do not provide a reliable estimate of the contribution of the $\Theta$-term to the nucleon electric dipole moments, and a factor of ten higher statistics data are needed to get better control over the systematics and possibly a $3\sigma$ result.


I. INTRODUCTION
The permanent electric dipole moments (EDMs) of nondegenerate states of elementary particles, atoms and molecules are very sensitive probes of CP violation (CP). Since the EDMs are necessarily proportional to their spin, and under time-reversal the direction of spin reverses but the electric dipole moment does not, a nonzero measurement confirms CP violation assuming CPT is conserved. Of the elementary particles, atoms and nuclei that are being investigated, the electric dipole moments of the neutron (nEDM) and the proton (pEDM) are the simplest quantities for which lattice QCD can provide the theoretical part of the calculation needed to connect the experimental bound or value to the strength ofCP in a given theory [1,2].
EDMs can shed light on one of the deepest mysteries of the observed universe, the origin of the baryon asymmetry: the universe has 6.1 +0. 3 −0.2 × 10 −10 baryons for every black body photon [3], whereas in a baryon symmetric universe, we expect no more than about 10 −20 baryons and anti-baryons for every photon [4]. It is difficult to include such a large excess of baryons as an initial condition in an inflationary cosmological scenario [5]. The way out of the impasse lies in generating the baryon excess * Electronic address: tanmoy@lanl.gov † Electronic address: cirigliano@lanl.gov ‡ Electronic address: rajan@lanl.gov § Electronic address: emereghetti@lanl.gov ¶ Electronic address: boram@lanl.gov dynamically during the evolution of the universe. But, if the matter-antimatter symmetry was broken post inflation and reheating, then one is faced with Sakharov's three necessary conditions [6] on the dynamics: the process has to violate baryon number, evolution has to occur out of equilibrium, and charge-conjugation and CP invariance have to be violated. CP violation exists in the electroweak sector of the standard model (SM) of particle interactions due to a phase in the Cabibbo-Kobayashi-Maskawa (CKM) quark mixing matrix [7], and possibly due to a similar phase in the Pontecorvo-Maki-Nakagawa-Sakata (PNMS) matrix in the leptonic sector [8,9]. The effect of these on nEDM and pEDM is, however, small: that arising from the CKM matrix is about O(10 −32 ) e cm [10][11][12], much smaller than the current 90% confidence level (CL) experimental bound d n < 1.8 × 10 −26 e cm [13], 1 and than the reach of ongoing experiments, d n < 3.4 × 10 −28 e cm at 90% confidence [15].
In principle, the SM has an additional source of CP violation arising from the effect of QCD instantons. The presence of these localized finite action nonperturbative configurations in a non-Abelian theory leads to inequivalent quantum theories defined over various 'Θ'-vacua [16,17]. Because of asymptotic freedom, all nonperturbative configurations including instantons are strongly suppressed at high temperatures [18,19] where baryon number violating processes occur. Because of this, CP violation due to such vacuum effects does not lead to appreciable baryon number production [20]. Nonetheless, understanding the contribution of such a term to the nucleon EDM is very important for two reasons. First, the Θ term constitutes a 'background' contribution to all hadronic EDMs that needs to be understood before one can claim discovery of new sources of CP violation through nucleon or hadronic EDM measurements; and second, besides generating higher-dimensional CPodd operators, new sources of CP-violation beyond the Standard Model (BSM) also generate a so-called 'induced Θ term' [1,21,22] if one assumes that the Peccei-Quinn mechanism is at work [23]. Therefore, in the large class of viable models of CP violation that incorporate the Peccei-Quinn mechanism, quantifying the contribution of the induced Θ to the nucleon EDM (operationally, the calculation is the same as in the first case) is essential to bound or establish such sources of CP violation.
Until recently, the calculation of hadronic matrix elements needed to connect nucleon EDMs to SM and BSM sources of CP violation relied on chiral symmetry supplemented by dimensional analysis [24][25][26][27][28][29][30][31][32] or QCD sum rules [1,22,[33][34][35][36], both entailing large theoretical errors. Large-scale simulations of lattice QCD provide a firstprinciples method for calculating these matrix elements with controlled uncertainties. Several groups have reported results of lattice QCD calculations of the neutron EDM induced by the QCD Θ term [37][38][39][40][41][42][43][44] and by higherdimensional operators, such as the quark EDM [45,46] and at a more exploratory level the quark chromo-EDM [47][48][49]. In this paper, we present a new calculation of the contribution of the Θ-term to the nEDM and pEDM and show that the statistical and systematic uncertainties are still too large to extract reliable estimates.
This paper is organized as follows: In Section II, we describe our notation by introducing the Lagrangian witḧ CP and the needed matrix elements. In Section III, we describe the decomposition of the matrix elements into the electromagnetic form factors. Section IV provides the lattice parameters used in the calculations. In Section V, we present the implementation of the gradient flow scheme, and in Sec.VI the calculation of the topological susceptibility. Section VII describes the methodology for extracting theCP phase α for the ground state created by the nucleon interpolating operator used, from the two-point function. This phase controls the CP transformation of the asymptotic nucleon state. Section VIII describes the calculation strategy for obtaining the form factors when this phase α is nonzero and gives the formulae used to extract theCP form factor F 3 from the matrix elements. In Section IX, we discuss the extraction of F 3 (q 2 ) and the removal of the excited states contamination. The extrapolation of F 3 (q 2 ) to q 2 = 0 is presented in Sec. X. Section XI discusses the lattice-spacing artifacts. Our results with the excited state spectrum taken from the two-point function are presented in Sec. XII and those with an N π excited state in Sec. XIII. These results are compared to previous calculations in Section XIV. Conclusions are presented in Section XV. Further details on the connection between Minkowski and Euclidean notation, the extraction of the form factors, the chiral extrapolation, excited-state contamination, and the O(a) corrections in the Wilson-clover theory are presented in five appendices.

II. THE QCD Θ-TERM
QCD allows for the existence of a P and T (andCP if CPT is conserved) violating dimension-four operator, i.e., the Θ-term. In its presence, the QCD Lagrangian density in Euclidean notation becomes where G a µν is the chromo-field strength tensor, G a µν = 1 2 µνλδ G aλδ is its dual, and Θ is the coupling. 2 G µν G µν is a total derivative of a gauge-variant current and its space-time integral gives the topological charge Non-zero values of Q are tied to the topological structure of QCD and the U (1) axial anomaly. In addition, higher dimension operators that arise due to novelCP couplings at the TeV scale generate this term under renormalization in a hard cutoff scheme like lattice regularization or gradient 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 Θ into a complex phase of the quark matrix and vice versa. It is, therefore, necessary to work with the convention independent Θ = Θ + Arg DetM q , which includes both, Θ from all sources and the overall phase of the quark matrix M q . Since, the argument of the determinant is ill-defined when it is zero, all physical effects of Θ vanish in the presence of even a single massless quark flavor.
If the overall Θ is nonzero, then this operator would induce an nEDM d n of size Here X is obtained from theCP part of the matrix element of the electromagnetic vector current within the neutron state in the presence of the Θ-term and F 3 is theCP violating form factor defined in Eq. (6). This is obtained, at the leading order, from theCP part of the matrix element where we have assumed that the Θ-term is the only source ofCP. In other words, X provides the connection between theCP coupling (Θ) and the nEDM (d n ). At present, the upper bound on the nEDM, |d n | < 1.8 × 10 −26 e cm (90% CL) [13], is used along with an estimate X ∼ (2.50±1.25)×10 −16 e cm [1] to set a limit on the size of Θ 10 −10 . This is an unnaturally small number! One solution to this unnaturalness is the dynamical tuning of Θ = 0 using the Peccei-Quinn mechanism 3 [23].
Our goal is to calculate X using lattice QCD, which multiplied by the cumulative value, Θ, from all sources (SM or BSM), gives the full contribution to nEDM from the dimension-4 GG operator in Eq. (1). Knowing X will allow current and future bounds on (or measured value of) d n to more stringently constrain or pin down Θ.
In the rest of the paper, all the analyses are carried out assuming that the onlyCP coupling arises from the Θ-term, whose strength is Θ. Results are presented for Θ = 0.2, which we have checked is small enough so that O(Θ 2 ) corrections are negligible for all quantities of interest (α and F 3 defined later). The lattice calculation consists of the evaluation of the connected and disconnected diagrams shown in Fig. 1. The disconnected diagram gets contributions from all The connected (left) and disconnected (right) diagrams with the insertion of the bilinear vector current (red filled circle) in the nucleon two-point function. The signal is given by the correlation between this 3-point function and the topological charge shown by the filled yellow circle.
3 The Peccei-Quinn mechanism relaxes Θ dynamically to Θ ind , the point where the effective potential achieves its minimum. In the absence of other sources of CP violation in the theory, Θ ind = 0.
quark flavors in the loop-but their contributions to the CP-conserving form-factors of the vector current are small [50]. In this work, we assume the same holds for the CP-violating ones and neglect the contribution to the electric dipole moment coming from these diagrams.

III. FORM FACTOR OF THE ELECTROMAGNETIC CURRENT
The parameterization of the matrix element of the electromagnetic current, J EM µ (q), defined in Eq. (5), within the nucleon state in terms of the most general set of form factors consistent with the symmetries of the theory is where M N is the nucleon mass, q = p − p is the Euclidean 4-momentum transferred by the electromagnetic current, σ µν = (i/2)[γ µ , γ ν ], and u N (p, s) represents the free neutron spinor of momentum p and spin s obeying (i / p + M N )u N (p, s) = 0, with γ 4 implementing the asymptotic (i.e., free) parity operation. Throughout, we work in Euclidean space and refer the reader to Appendix A for details on our conventions. F 1 and F 2 are the Dirac and Pauli form factors, in terms of which the Sachs electric and magnetic form factors are G E = F 1 − (q 2 /4M 2 N )F 2 and G M = F 1 + F 2 , respectively. 4 The anapole form factor F A and the electric dipole form factor F 3 violate parity P; and F 3 violates CP as well. The zero momentum limit of these form factors gives the charges and dipole moments: the electric charge is G E (0) = F 1 (0), the magnetic dipole moment is G M (0)/2M N = (F 1 (0) + F 2 (0))/2M N , and the EDM is defined in Eq. (4).
In all the discussions in this paper, the current J EM µ used is the renormalized local vector current Z V i e i ψ i γ µ ψ i , where e i is the electric charge of a quark with flavor i. The renormalization is carried out by taking ratios of all three-point fermion correlators with the lattice estimate of the vector charge, g V ≡ 1/Z V , which is given by the forward matrix element of ψ i γ 4 ψ i . These ratios are constructed with identical source, sink, and current insertion positions and within the single jackknife loop used for the statistical analysis of the data to  take advantage of error reduction due to correlated fluctuations. 5

IV. LATTICE PARAMETERS
We present results on seven ensembles, whose parameters are defined in Table I. These were generated by the MILC collaboration [51] using 2+1+1-flavors of highly improved staggered quarks (HISQ) action. For the construction of the nucleon correlation functions we use the clover-on-HISQ formulation that has been used extensively by us in the calculation of the nucleon charges and form factors as described in Refs. [52,53]. These ensembles cover three values of the lattice spacing, a ≈ 0.12, 0.09 and 0.06 fm and three values of the pion mass M π ≈ 315, 220 and 130 MeV. Further details of the lattice parameters and methodology, statistics, and the interpolating operator used to construct the nucleon 2-and 3-point correlation functions can be found in Refs. [52,53].

V. TOPOLOGICAL CHARGE UNDER GRADIENT FLOW
We calculate the topological charge using the gradient flow scheme to implement operator renormalization and to reduce lattice discretization effects [41,54]. The primary advantage of the scheme is that at finite flow times 6 , i.e., for τ gf > 0, the flow time provides an ultraviolet cutoff, and the continuum limit, a → 0, of all operators built solely from gauge fields is finite. Moreover, since topological sectors arise dynamically as we take the continuum limit, the gradient flowed topological charge takes on integer values, and no renormalization is needed to convert it to a scheme that preserves this property; in particular, correlators of the topological charge are flow-time independent [54].
These statements are, however, not true at finite lattice spacing and volume. At small τ gf , we get O(a 2 /τ 2 gf ) artefacts. In Fig. 2, we show the distribution of the topological charge Q as a function of the flow time τ gf in physical units. Its distribution has stabilized by τ gf = 0.24 fm for the 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 are smoothed out, indicating that they are lattice artifacts.
In Fig. 3, we show the distribution of the difference from the nearest integer. This distribution stabilizes more slowly and it is only by τ gf = 1.31 fm (τ gf = 0.76 fm) on the a ≈ 0.12 fm (a ≈ 0.09 and 0.06 fm) ensembles that the charges are close to integers. The relevant distribution important for the calculation of the nucleon correlation functions is, however, likely to be the distribution of Q shown in Fig. 2. To explore this, we show in Fig. 4 the value of F 3 as a function of τ gf for the a ≈ 0.12 and 0.09 fm ensembles, and find that indeed the correlation functions, and thus F 3 , do stabilize early but the 6 We use the notation τ gf ≡ √ 8t for the flow time, where t is the parameter in the flow equations in Ref. [54]. We used the Runge-Kutta integrator given in that reference for integrating the flow equations, with a step size of 0.01. Changing the step size to 0.002 changed the results on topological susceptibility by less than 0.2%. 75   τ gf required for the coarser lattices is longer. Thus, to be conservative, the results presented below are obtained with flow times τ gf (a06) = 0.68 fm, τ gf (a09) = 0.68 fm and τ gf (a12) = 0.86 fm respectively. In Fig. 5, we show the distribution of the nearest integer, Q int , to the topological charge at τ gf ≈ 1.4 fm (τ gf = 0.76 fm) on the a ≈ 0.12 fm (a ≈ 0.09 and 0.06 fm) ensembles, by which time the Q int identified with a given configuration has stabilized. This distribution is approximately symmetric about zero as expected since Q = 0, and no gaps are visible in the distribution. In Fig. 6, we show the autocorrelation function of Q versus    the flow time. The data show no significant change after τ gf 0.3 fm, so we can determine the autocorrelation from these data. We do not observe a long time freeze in Q in any of the ensembles analyzed as illustrated using the a09m130 and a06m135 ensembles at flow time τ gf = 0.68 fm in Fig. 7. The autocorrelation is less than about 10 configurations for all but the a06m135 ensemble. Based on this study, the bin size used in the single elimination jackknife procedure is given in Table I.

VI. TOPOLOGICAL SUSCEPTIBILITY
The topological susceptibility χ Q is defined as Its value in the pure gauge theory, χ quenched Q , is related to the mass of the η meson in a theory with N f light flavors in the chiral limit via the axial anomaly, viz., the

Witten-Veneziano relation [55, 56]
where F π is the pion decay constant in the convention where its physical value is about 93 MeV. Following Ref. [57], we can include the effects of the quark masses. Including SU (3) breaking at leading order in χPT but neglecting the heavier quarks gives With dynamical fermions, however, the susceptibility should vanish in the chiral limit. For SU (N f ) flavor group with finite but degenerate quark masses, it should behave as [58][59][60]: For N f = 2 light flavors and the strange quark, but neglecting the heavier quarks that give negligible corrections, leading order chiral perturbation theory (χPT) modifies this to We calculate χ Q on the 2+1+1 flavor HISQ ensembles, which are O(a) improved. The results are given in Table I. In addition to the seven ensembles used to calculate F 3 , we include data from the a06m310 and a06m220 ensembles. We remind the reader that the MILC collaboration has previously highlighted the issue of frozen topology on these ensembles [61], which is why we do not use them in the calculation of F 3 .
As discussed in Section V, the topological susceptibility at finite flow time needs no renormalization, and should be independent of flow time up to O(a 2 /τ 2 gf ) effects. 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 a12m220 and a12m220L ensembles, and show that this is a τ 2 gf /L 2 effect, where L is the lattice size. 7 At the flow times and volumes we use in the calculation, this is a small effect and therefore neglected.
To obtain χ Q at M π = 135 MeV and a = 0, we use the fit ansatz which assumes χ Q is zero in the chiral-continuum limit. We do not find a viable χ 2 /dof on including all nine data points. Reasonable fits are found on neglecting (i) all three a ≈ 0.12 fm points and (ii) all three a ≈ 0.12 fm and the a06m310 point. These two fits give χ Q = [70 (6) MeV] 4 and χ Q = [63 (9) MeV] 4 , respectively, at M π = 135 MeV. We take the average χ Q = [66(9)(4) MeV] 4 as our best estimate, the larger of the two errors and an additional systematic uncertainty, which is half the difference. These results are in good agreement with the expected value, (79 MeV) 4 , obtained using the physical meson masses and decay constants in Eqs. (9) and (11). The data and the fit case (i) are shown in Fig. 10.

VII. CALCULATION OF THECP PHASE α
In a field theory in which parity is not conserved, the definition of parity of a composite state, e.g., the neutron state, needs care [22,33,37]. To explain this, we start with the most general spectral decomposition of the timeordered 2-point nucleon correlator where A i is the amplitude for creating state i, E i is its energy, the Euclidean time τ is the separation between the source and the sink, and, for notational convenience, we are assuming a discrete spectrum. A common choice on the lattice of the neutron interpolating operator N is where C = γ 2 γ 4 (the sign is conventional and does not affect the nucleon correlators we study; see Appendix A for details of our convention) is the charge conjugation matrix, a, b, c are the color indices and u, d are the quark flavors. The 4 × 4 spinor matrix M s i in Eq. (13) depends on the state and the momentum p. Its most general form consistent with Lorentz covariance is 8 where p 4 i ≡ iE i . It is clear that because of the presence of the phases α i , the parity operator that transforms the spinor associated with the i th asymptotic state is P αi ≡ e iαiγ5 Pe −iαiγ5 , where P ≡ ηγ 4 is the usual parity operator for a particle with intrinsic parity η. The phases α i depend on the realization of discrete symmetries: If the interpolating field is chosen such that P implements parity in the free theory, Im α i = 0 for a PT symmetric theory, Re α i = 0 for the CP symmetric theory, α i = 0 for a P symmetric theory. For our case of onlyCP, all α i are, therefore, real, which will be implicit except in Appendix B. It is important to note that the value of α i depends on the interpolating operator N , the state, and the source ofCP. Its value for the ground state can be extracted from the large τ behavior of the imaginary part of the nucleon 2-point function. Consider Keeping only the first two states one gets The data for r α versus τ are shown in Fig. 11 for all seven ensembles. The α 0 for the ground state obtained from the two-state fit agrees with the plateau at large τ , where the lowest state dominates, and is independent of momentum.

VIII. THREE-POINT FUNCTIONS IN THE PRESENCE OF THE PHASE α
In the presence of the phase α N ≡ α 0 corresponding to the ground-state nucleon [47], the most straightforward way to extract the matrix element of the electromagnetic current J EM µ within the neutron ground state in the presence ofCP is to calculate the correlation function a12m310 a12m220 a12m220L a09m310 a09m220 a09m130 a06m310 a06m220 a06m135 Extrap FIG. 10: Fits to the data for the topological susceptibility, χQ, using the ansatz given in Eq. (12). . It is a Lorentz scalar and independent of the momentum as confirmed by the lattice data. The χ 2 /dof values presented are from fully correlated fits, except for the case of the a09m220 and a06m135 ensembles on which we use uncorrelated fits to avoid instabilities.
where p ≡ p + q and Here, the current J EM µ is inserted at times t between the neutron source and sink operators located at time 0 and τ , and a sum over the spin labels is implicit. We also assume that t and τ are large enough that only the ground state dominates the correlation function. This form results from the realization that γ 4 remains the parity operator for the ground state nucleon when working with the interpolating field defined to be e −iα N γ5 N instead of N in all correlation functions.
This approach, however, requires, evaluating the full 4 × 4 matrix of 3-point correlation functions. In our calculation, we have implemented the spin projection using so the contribution of a nonzero α N has to be incorporated at the time of the decomposition of the matrix element 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 nucleon ground state contribution to the matrix element of the electromagnetic current, Tr e iα N γ5 P 3pt e iα N γ5 where O µ is given in Eq. (23). The full expressions for V 1,2,3,4 , along with a general strategy for extracting F 3 , from the four coupled complex equations is given in Appendix B. To extract F 3 , theCP part of the three-point functions, a very significant simplification of the analysis and improvement in the signal is achieved by subtracting the Θ = 0 contribution from each component of the current in Eq. (25) before making the excited state fits and decomposing the resulting ground state matrix element in terms of form factors. This is implemented by analyzing the ground state contribution in terms of the combina-tionV µ = V µ (Θ) − V µ (0). Working to first order in Θ, and recalling that s α N ≡ sin α N cos α N ∼ α N ∼ O(Θ), and F 3 ∼ O(Θ), the expressions for the ground state contributions of the three-point functionsV 1,2,3,4 in terms of form factors simplify tō where G 1 = F 1 + F 2 and G 3 = F 3 + s α N F 2 . We solve the above system for G 1 and G 3 . At q 2 = 0 there is a further simplification because G 1 (0) = Q N + F 2 (0) where Q N is the nucleon charge. With this, we get Though the nucleon anomalous magnetic moment G 1 (0) − Q N = F 2 (0) ≡ κ N has been measured very precisely, the largest contribution to G 3 comes from s α N F 2 , and the statistical error is much smaller when extrapo- , rather than extrapolating only G 3 (q 2 ) and then combining it with s α N µ N to get the right hand side of Eq. (27). Also, note that G 3 (q 2 ) can be obtained uniquely fromV 1 andV 2 for a number of values of q 2 , which provides a useful check.
One can extend Eq. (27) to definẽ To get F 3 (0) =F 3 (0), we find better control by extrapo-latingF 3 (q 2 ) to q 2 → 0. The subtraction of the Θ = 0 contribution also allows averaging of the three point functions over momenta related by cubic invariance, as seen by comparing the simpler Eqs. (26) with Eqs. (B8). We illustrate the improvement in the signal in Fig. 12. The averaging over equivalent cases (over momenta related by cubic symmetry and over V 1 and V 2 ) significantly reduces the statistical errors and improves the analysis of excited state contamination (ESC) discussed next.

IX. REMOVING ESC IN F3
In order to extract the ground state contributionV µ from lattice data on the ratio R µ (τ, t, q) of three-and two-point functions defined in Eq. (B7), we need to remove all excited states that make a significant contribution.
We have analyzed data on R µ (τ, t, q) in terms of a twostate fit, following two strategies. In the first, we have taken the first excited-state energies from a three-state fit to the two-point function. In the second strategy, we have set the first excited-state energy to the non-interacting energy of the N π state, motivated by the χPT expectation that the leading excited state is the N π state, with amplitude of the same size as the ground state contribution (see Appendix D for more details). In Fig. 13 we compare the two strategies for Im(R 4 (τ, t, q)). The χ 2 /dof of the fits are similar for the two cases on all three ensembles, but the ground state estimate is vastly different and thus the contribution to the nEDM. With  q=(1,1,1)  the current data, picking between them is the key unresolved challenge for this calculation. The very large extrapolation for τ → ∞ in the N π case, however, leads us to question whether a two-state fit is sufficient if the N π state is included and whether a similar effect might contaminate our extraction of α N . We therefore first perform the analysis taking the excited state energy, E 1 , from a three-state fit to the two-point function and return to an analysis including a N π state in Sec. XIII. A second issue arising from the small signal in F 3 is that two-state fits to many of the correlation functions with the full covariance matrix are unstable with respect to variations in the values of τ and t skip , the number of points skipped in the fits adjacent to the source and sink for each τ . Examples of this are shown in Fig. 14 for Re (R 1 (τ, t, q)). This has two consequences for the analysis. First, we have carried out the final analysis using only the diagonal elements of the covariance matrix. We have, however, checked that in cases where fully covariant fits are possible, the two results are consistent. Since we use uncorrelated fits for removing excited-state contamination, we do not quote a χ 2 /dof for these fits. Second, the system of four equations, Eqs. (26), over determines G 3 and G 1 . While we solve the full set of equations as explained in appendix B, the data from Re(R 1,2 (τ, t, q)), which have poor signal, do not make a significant contribution. We have checked this by removing them from the analysis and the results are essentially unchanged, i.e., the results are dominated by Re(R 3 (τ, t, q)) and Im (R 4 (τ, t, q)).
The ansatz used to extrapolate F 3 (q 2 ) to q 2 → 0 is given in Eq. (C1) with one caveat. We useF 3 (q 2 ), defined in Eq. (28), instead of F 3 (q 2 ) as they are consistent to leading order and the extraction ofF 3 (q 2 ) is better controlled. We examine three fits based on Eq. (C1): • Linear: the quantities d i and S i are free parameters and H i is set to zero. • χPTg0: Same as χPT except g 0 is left as a free parameter.
The data and fits for the neutron and proton are presented in Figs. 15 and 16. The data are, within errors, flat in all cases and the extrapolated values from the three types of fits are consistent. Since in most cases, we have reliable data at only three values of q 2 , we take the final result from the χPT fit. At the end, we will take the difference between the Linear and χPT fits to estimate the associated systematic uncertainty.

XI. ADDITIONAL O(a) ARTIFACTS
Before performing a chiral-continuum extrapolation of the results, in this section we justify our continuum extrapolation 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 0 q . There are multiple sources of O(a) corrections that we need to consider. First, since our clover coefficient c SW is set to its tadpole-improved tree-level value, the action, 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 approximation is extremely good, and these are expected to be tiny effects. Second, the vector current we insert is not improved [62], and, hence, we expect its renormalization coefficient to have O(am q ) corrections. Such multiplicative terms, however, are unimportant near the chiralcontinuum limit, where theCP form factors vanish. A third source of O(a) effects is the required improvement of the vector current by an O(am 0 q ) mixing with the derivative of the tensor current, which can give rise to a nonzero F 3 , but only in the presence of CP violation in the theory. Since the topological charge does not introduce CP violation in the chiral limit, we would expect the behavior of d n to be dominantly O(a 2 ) in the chiral limit, if these were the only O(a) effects.
In Appendix E, we analyze the Wilson-clover theory based on the framework of a continuum EFT for the lattice action and the axial Ward Identities. Following Refs. [63][64][65], we show that the topological charge gives O(a)CP corrections, and identify this as effectively due to the insertion of the isoscalar quark chromo-EDM operator, which the topological term can mix with. Since this term is expected to survive in the chiral limit, we include an O(am 0 q ) term in our chiral continuum fits.

XII. CHIRAL-CONTINUUM EXTRAPOLATION AND RESULTS
In this section, we present the chiral-continuum (CC) extrapolation of data for d n (and, similarly, d p ) obtained on the seven ensembles. For each, we examine four cases.
These consist of two CC fits, Linear and χPT, using the leading order terms FIG. 16: The extrapolation ofF3(q 2 ) to q 2 → 0 using Eq. (C1) for the proton. The rest is the same as in Fig. 15. where the term c 3 a is the O(a) effect discussed in Sec. XI, because of which d n,p do not vanish in the chiral limit at finite a. The ansatz are distinguished by the terms proportional to c 2 (Linear) and c 2L (χPT). In these fits, M N is set to its physical value 940 MeV. We make these two fits to the data for d n,p obtained using (i) the linear and (ii) χPT extrapolation in q 2 , which leads to four estimates. These four CC fits for the neutron and the proton are shown in Figs. 17 and 18. The results and the fit coefficients c i are given in Table II. As discussed in Appendix C, at NLO in χPT the coefficient of the chiral logarithm c 2L is fixed in terms of the isovector scalar charge, the quark condensate and the pion decay constant, leading to (c 2L ) n = −(c 2L ) p = 0.033 fm-GeV 2 . Although the central values of the fits are approximately one order of magnitude larger, our results are compatible with this estimate at the 1σ-2σ level.
For the central value we take the χPT(q 2 )|χ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 = −0.003(7)(20)Θ e · fm (31) d p = 0.024(10)(30)Θ e · fm (32) where the second systematic error is the spread in the four estimates given in Table II.

XIII. ANALYSIS INCLUDING THE N π EXCITED STATE
In this section, we describe how all ground state quantities change when the N π excited state is included. This analysis should be considered exploratory because (i) the extrapolations in the fits to remove ESC (see Fig. 13), (ii) the errors, and (iii) the cancellations when combining different terms to get F 3 using Eqs. (26) are all large.
In Fig. 19, we show the increase in the value of α for the two physical mass ensembles as compared to the data presented in Fig. 11. The q 2 behavior is similar to that shown in Figs. 15 and 16, and the final results for the four strategies are given in Table II. The CC fits for the neutron and the proton using the χPT(q 2 )|χPT(CC) strategy are shown in Fig. 20 For the central value we again take the χPT(q 2 )| χPT(CC) result and the full spread for the error. This       The chiral-continuum extrapolation of dn using the ansatz given in Eq. (30). The four rows show (i) a linear CC fit to the data obtained using a linear extrapolation in q 2 discussed in Sec. X; (ii) a linear CC fit to the data obtained using the χPT extrapolation in q 2 ; (iii) a χPT CC fit to the data obtained using a linear extrapolation in q 2 ; and (iv) a χPT CC fit to the data obtained using the χPT extrapolation in q 2 . All data are with Θ = 0.2.      The chiral-continuum extrapolation of dp using the ansatz given in Eq. (30). The rest is the same as in Fig. 17.  gives d n | N π = −0.028(18)(54)Θ e · fm (33) d p | N π = 0.068(25)(120)Θ e · fm (34) where the second systematic error is the spread in the four estimates given in Table II.

XIV. COMPARISON TO PREVIOUS WORK
There are two estimates [44,66] of the contribution of the Θ-term to the nEDM since the clarification of the impact of the phase α that arises in the nucleon spinor in a theory withCP in Ref. [47]. That work also contains a review of previous results, which after correction were consistent with zero. No estimate is given in Ref. [47], but there is a preliminary value in a subsequent conference proceedings, Ref. [67]. All three of these calculations use the small Θ expansion and gradient flow method for topological charge renormalization as in this work. All results are summarized in Table III.
The ETM collaboration [66] has performed the calculation on one 2+1+1-flavor twisted mass clover-improved ensemble with a = 0.0801(4) fm, 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-size effects. They also implicitly implement the Θ = 0 subtraction (see Eqs. (26)) that we find reduces the statis-tical noise by using the spin projector (1 + γ 4 )iγ 5 γ k /4. They determine F 3 (0) by making a constant fit to the lowest three q 2 points. Their final result is taken using the spectral projectors method, which they find reduces the errors by a factor of about two compared to the fieldtheoretic definition of the topological charge used in this work. They do not, however, assess a systematic error associated with excited-state effects, extrapolation in q 2 , or the chiral-continuum fit.
The calculation presented in Ref. [44] uses six 2+1flavor Wilson-Clover ensembles but only one below M π = 567 MeV, with M π = 410 MeV. The values of lattice spacings range between 0.068 < a < 0.11 fm. A linear fit in q 2 is made to obtain F 3 (0). Also, artifacts due to ESC are not analyzed and, in any case, data with the heavy pion masses studied, M π > 410 MeV, would not provide sensitivity to analyses with or without including a N π state. This is the only other calculation that has presented a chiral extrapolation using the χPT ansatz (Eq. (30) but with a O(a 2 ) discretization correction instead of our c 3 a term). As shown in the bottom right panels in Figs. 17 and 18, such chiral fits have an inflection point close to the smallest M π data point in order to satisfy the constraint F 3 = 0 at M π = 0. In the case of Ref. [44], this occurs around M π = 400 MeV, raising questions on the reliability of the extrapolation.    The chiral-continuum extrapolation of dn (top) and dp (bottom) using the ansatz given in Eq. (30), using N π as the excited state fits, and with the χPT(q 2 )|χPT(CC) strategy. All data are with Θ = 0.2.

XV. CONCLUSIONS
This paper presents a calculation of the contribution of the Θ-term to the nucleon electric dipole moment using 2+1+1-flavor HISQ ensembles and Wilson-clover valence quarks. Two of the seven ensembles analyzed are at the physical pion mass, which anchor our chiral fits. The calculation has been done using the small Θ expansion method. Significant effort has been devoted to getting a reliable signal in theCP violating form factor F 3 . The gradient flow scheme has been used to renormalize the Θterm and the results are shown to be independent of the flow time. Our estimate of the topological susceptibility for the 2+1+1 theory is χ Q = (66(9)(4) MeV) 4 in the continuum limit at M π = 135 MeV.
We also present two technical issues. First, in Appendix D, we show that, in chiral perturbation theory, the N π excited state should provide the dominant contamination. We have, therefore, used two strategies for removing excited state contamination. In the first, the mass gaps are taken from fits to the spectral decomposition of the nucleon two-point function, and in the second we assume they are given by the non-interacting energy of the N (0)π(0) state. We find a very significant difference between the two as shown in Secs. IX and XIII, and by the results summarized in Tables II and III. The second technical issue discussed in Sec. XI and appendix E is that lattice artifacts introduce a term proportional to am 0 q , because of which d n does not vanish in the chiral limit at finite a. Our chiral-continuum fits have been made including this term.
The analysis of the q 2 dependence of F Θ 3 has been carried out using both a linear and the leading order χPT expression as described in Sec. X. The current data do not distinguish between the two. Similarly, the chiral fit is also carried out using a linear and the leading order χPT expression as described in Sec. XII. The results from these four sets of fits and the two strategies to remove excited-state contributions are summarized in Table II. Our preferred values are obtained using the leading order χPT expressions. The analysis using excited states from fits to the two-point function indicate that d Θ n is small, |d Θ n | 0.01 Θ e · fm, whereas for the proton we get |d Θ p | ∼ 0.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 | ∼ 0.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 Table III, which also includes estimates from previous works, it is clear that, at present, lattice calculations do not yet provide a reliable estimate. To improve the current 100% uncertainty to a 3σ result will require a factor of at least ten improvement in statistics.
To make our conventions explicit, we present the connection between Minkowski and Euclidean variables in Table IV.
To connect the Lagrangian density for the Θ term in Minkowski and Euclidean spaces, we take the Minkowski action associated with the QCD Θ term to be (A2) The factor of +i arises from the transformation of the field strength and because each term in the sum has one factor of G 0i (or G i0 ) and one factor of G jk . Moreover, we used which implies and hence ( E ) 1234 = +1.
Minkowski gamma matrices are unitarily transformed from the standard chiral basis, γ µ c [70]  Putting together the change in the measure and the change in the Lagrangian density we have consistently with Eq. (1).

Appendix B: Extraction of F3
The Euclidean four-vector V µ (q) defined in Eq. (25) can be determined from lattice data by taking appropriate ratios of 3-pt function and 2-pt functions. This is achieved by defining the projected 2-and 3-point functions as follows, with q = p − p, p = 0, P 3pt given in Eq. (24), and, neglecting the contributions of heavier quarks, The ratiõ becomes independent of t and τ if t, τ are sufficiently large that excited state effects can be neglected, and takes the form .
(B6) In our plots to demonstrate the signal and excited states, we, therefore, choose to show the quantity where g V ≡ C µ 3pt (τ, t, 0)/C 2pt (t, 0), and α N is calculated from fits to the 2-pt functions with momentum p or p as discussed in Section VII.
From the above expressions we want to extract F 3 (q 2 ), that gives the neutron EDM. It turns out that the RHS of Eqs. (B8) is most naturally expressed in terms of G 1,2,3 given by where q 2 E = q 2 + q 2 4 and s ≡ sin α cos α, c ≡ cos 2 α. For a given momentum transfer q = (q 1 , q 2 , q 3 ), Eqs. (B8) thus represents eight equations for G 1,2,3 . They can be written in a compact form as follows: where K(q) is an 8 × 3 matrix given in block form by and V (q) is an eight-dimensional array given by .
To solve for G 1,2,3 (q 2 ), for a given three-momentum transfer q = (q 1 , q 2 , q 3 ) we can use a least squares esti-mator. Namely, we minimize the function where the weights matrix is the inverse of the covariance matrix of lattice "measurements" V i (q): For independent variables V i (q) the covariance matrix C V and it inverse are positive definite. 9 This guarantees that F (G 1,2,3 ) is minimized if and only if E i (q) = 0 for all i.
(B18) or even more explicitly which is a system of three equations for G 1,2,3 . The extremum condition for F (G 1,2,3 ) implies the following linear equation for G 1,2,3 (q 2 ): where the 3×3 matrix A and the three dimensional array B are given by K jα (q) w ji (q) K iβ (q) (B21a) 9 For ease of notation, we are ignoring current conservation, which relates the various components V i (q). Strictly speaking, we need to eliminate the dependent components of V i (q) when using a conserved current to get an invertible covariance matrix.
So from the lattice data on V i (q), their covariance matrix, and the explicit form of the matrix K iα (q) given in Eq. (B11) one can construct A αβ and B α and solve for G 1,2,3 . Error on G 1,2,3 can be assigned with the bootstrap method.
with g S = 1.0 and B = 2.8 GeV.d 0,1 are two low-energy constants, which, by naive-dimensional-analysis, scale as The first derivative of the form factor is [28,30,32] S 0 = 0, (C10) At N 2 LO there are additional long-and short-distance contributions to both isoscalar and isovector components. The remaining momentum dependence of the EDFF is given by the functions H i (q 2 ) introduced in Eq. (C1), . h a appears at leading order, Since these behave as h In this appendix, we show that, in χPT, the gap between the ground state and excited state contributions to the CP-odd components of the three-point function C µ 3pt is expected to be of order of the pion mass M π . This can be intuitively understood from the fact that the nucleon EDM induced by the QCDΘ term receives a LO contribution from a long-range pion loop [24], shown in Figure  21. In Minkowski space, this diagram has a branch cut when the intermediate pions and nucleon go on-shell. In Euclidean space, this translates into a N π excited state, whose amplitude is of the same size as the ground state contribution. For simplicity, we focus only on the diagram shown in Figure 21, and assume that the nucleon interpolating field does not couple to nucleon plus pions. We start from the 4 th component of the three-point function. Carrying out the Dirac traces in Eq. (25), in the limit M N q, we find where π + q 2 and . . . denotes terms with a gap with two or more units of momentum. f 0 (M π , q, L) denotes the ground state loop function, which we write as an infinite volume term f ∞ 0 and a correction ∆ f 0 (M π , q, L) = f ∞ 0 (M π , q) + ∆(M π , q, L). (D2) In the non-relativistic limit, f ∞ 0 (M π , q) is given by and is ultraviolet divergent. In dimensional regularization and in the MS scheme f ∞ 0 (M π , q) = log with x = q 2 /(4M 2 π ), which is of course the same function as in Section C. The finite volume correction is given by ∆(M π , q, L) = (4π) 2 dk 0 2π which can be written in terms of Bessel functions as [72] ∆(M π , q, L)= At q = 0, for M π L ∼ 4, ∆ amounts to a 0.1% correction. Eqs. (D1) and (D4) thus show that the excited states have a gap of O(M π ). The ratio of the ground and excited state contributions is determined by the quantity (4π) 2 /(LM π ) 3 , which is a number of order 1 for LM π = 4. We thus do not expect a significant suppression of the excited states. A similar calculation can be performed for the spatial components C i 3pt , yielding a result similar to Eq. (D1), but with a sinh rather than cosh behavior. terms in the continuum limit of axial WIs involving composite fields [63,64]. It is, however, essential in order to identify the O(a) corrections to d n (Θ). Using the above expression in (E6), and taking into account the mixing between (GG) and ∂ µ A µ (which involves the renormalization constant Z C ) one arrives at [64,65] O(x 1 , ..., where is the quark mass free of power divergences as we take the continuum limit. Here, and henceforth, the O(ma) dependence of the coefficients of the operators are suppressed. Finally, upon integrating over d 4 x we arrive at .
(E12) Ref. [74] performed a detailed diagrammatic analysis of Eq. (E12), with O(x 1 , x 2 .x 3 ) = N (x 1 ) J EM µ (x 2 )N (x 3 ) in the a → 0 case, showing that the δO terms cancel the connected insertions of 2mψγ 5 ψ. Their analysis shows that insertions of the operator GG can be replaced by 2m times the disconnected insertions of the isosinglet pseudoscalar densityψγ 5 ψ. Since the disconnected matrix elements of the isoscalar density do not diverge in the chiral limit, this implies as a corollary that the neutron EDM should vanish as m → 0. O(a) effects would modify the result of Ref. [74] by modifying the RHS of their Eqs. (2.11) and (3.5). In the context of our analysis, the term proportional to aX in Eq. (E12) provides O(a) effects, which we discuss next.
First, we project the subtracted operatorX on the basis of (subtracted) dim-5 operators, given in Ref. [75], whereσ µν ≡ 1 2 (σ µν γ 5 + γ 5 σ µν ) and ψ E = ( / D +m)ψ. Keeping in mind that O(x 1 , ..., x n ) has the structure N (x 1 ) J EM µ (x 2 )N (x 3 ), in terms of the neutron source and sink operator and the electromagnetic current, the various O (5) n contribute to Eq. (E12) as follows: • O • O 3,4 involve one and two powers of the electromagnetic field strength. In order to eliminate the photon field in the correlation functions in Eq. (E12), one needs electromagnetic loops, making the contribution of O (5) 3,4 to Eq. (E12) of O(a α EM /π), and thus negligible to the order we are working.
• O 11 contains two equation of motion operators. Therefore, when inserted in Eq. (E12), it will always involve a contraction with a quark field in the neutron source or sink operator, and thus it will not contribute to the residue of the neutron pole. O (5) 12 is a total derivative and drops out of Eq. (E12). O (5) 13 is gauge-variant operator and drops out of Eq. (E12) as long as O(x 1 , ..., x n ) is a gauge singlet, which is the case for O(x 1 , x 2 , x 4 ) ∝ N (x 1 ) J EM µ (x 2 )N (x 3 ). O