Neutron electric dipole moment using lattice QCD simulations at the physical point
NNeutron electric dipole moment using lattice QCDsimulations at the physical point
C. Alexandrou ( a,b ) , A. Athenodorou ( c,b ) , K. Hadjiyiannakou ( a,b ) A. Todaro ( a,d,e ) ( a ) Department of Physics, University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus ( b ) Computation-based Science and Technology Research Center, The Cyprus Institute, 20 Kavafi Str., Nicosia 2121, Cyprus ( c ) Dipartimento di Fisica, Universit`a di Pisa and INFN, Sezione di Pisa, Largo Pontecorvo 3, 56127 Pisa, Italy ( d ) Institut f¨ur Physik, Humboldt-Universit¨at zu Berlin, Newtonstr. 15, 12489 Berlin, Germany ( e ) Dipartimento di Fisica, Universit`a di Roma “Tor Vergata”, Via della Ricerca Scientifica 1, 00133 Rome, Italy
Abstract
We extract the neutron electric dipole moment | (cid:126)d N | within the lattice QCD formalism. We analyse oneensemble of N f = 2 + 1 + 1 twisted mass clover-improved fermions with lattice spacing of a (cid:39) .
08 fm andphysical values of the quark masses corresponding to a pion mass m π (cid:39)
139 MeV. The neutron electricdipole moment is extracted by computing the CP -odd electromagnetic form factor F ( Q →
0) throughsmall θ -expansion of the action. This approach requires the calculation of the topological charge for whichwe employ a fermionic definition by means of spectral projectors while we also provide a comparison withthe gluonic definition accompanied by the gradient flow. We show that using the topological charge fromspectral projectors leads to absolute errors that are more than two times smaller than those providedwhen the field theoretic definition is employed. We find a value of | (cid:126)d N | = 0 . θ e · fm when usingthe fermionic definition, which is statistically consistent with zero.November 3, 2020 a r X i v : . [ h e p - l a t ] N ov Introduction
The discrete symmetries of parity P , charge conjugation C and time-reversal T play a significant role in thephenomenological structure of the Standard Model (SM). Although the neutron has an overall zero electriccharge, a conundrum that physicists are trying to understand for decades is whether the distribution of thepositive electric charge could possibly not coincide with the distribution of the negative electric charge. Thatis to say whether there is no invariance under CP parity, which can address the well known unsolved puzzleof the origin of the imbalance of matter and antimatter in the universe. Such an imbalance between positiveand negative electric charge would manifest it self as the non-vanishing of the neutron electric-dipole moment(nEDM).Up to the present, no finite nEDM value has been measured in experiments. In addition, reported currentbounds are still several orders of magnitude below the SM prediction on CP violation induced by the weakinteractions. A finite nEDM value would point towards physics beyond the standard model (BSM) [1] andit is thus an interesting quantity to study theoretically.Experimental measurements are under way to improve the upper bound of the value of nEDM denoted as (cid:126)d N provided by a number of experiments, such as those given in Refs. [2–4]. Until recently, the best measuredupper bound was that given in Refs. [3, 4] as | (cid:126)d N | < . × − e · fm (90% CL) . (1)This result is extracted using stored “ultra-cold” neutrons, applying a weak magnetic field when a strong,parallel background electric field is reversed and measuring the alternation on the neutron spin precessionfrequency. The experiment has been carried out at the Institut Laue-Langevin (ILL) reactor in Grenoble.Very recently [5], an experiment measured the nEDM at the Paul Scherrer Institute (PSI) in Switzerlandusing Ramsey’s method of separated oscillating magnetic fields with ”ultra-cold” neutrons. The novelty ofthis experiment lies on using the Hg-199 co-magnetometer and an array of optically pumped cesium vapormagnetometers to cancel and correct for magnetic field changes. The result of this experiment yields animproved upper limit of | (cid:126)d N | < . × − e · fm (90% CL) . (2)The nEDM value has also been extracted from different theoretical perspectives, such as model dependentstudies [6–14], as well as effective field theory calculations [15–23], reporting values for d N in the range of | (cid:126)d N | ∼ θ · O (cid:0) − − − (cid:1) e · fm. Using either the result given in Eq. (1) or Eq. (2) one derives a bound ofthe order θ (cid:46) O (cid:0) − − − (cid:1) .In this work, we investigate the nEDM induced by strong interactions, namely by the topological term.We neglect contributions which might come from higher dimensional components in the action, such as thelowest-dimensional ( d = 5) effective quark-gluon CP -odd interaction [24, 25]. In this approach, we include a θ -topological term in the QCD Lagrangian density: L QCD ( x ) = 12 g Tr [ F µν ( x ) F µν ( x )] + (cid:88) f ψ f ( x ) ( γ µ D µ + m f ) ψ f ( x ) − iθq ( x ) , (3)written in Euclidean time. The first two terms are CP -conserving while the θ -term is CP -violating that cangive rise to a non-zero nEDM. In the above expression, ψ f denotes a fermion field of flavor f with bare mass m f , F µν is the gluon field tensor and q ( x ) is the topological charge density, which in Euclidean space, isdefined as q ( x ) = 132 π (cid:15) µνρσ Tr [ F µν ( x ) F ρσ ( x )] . (4)1y considering the electroweak (EW) sector of the SM, the action in Eq. (3) receives a contribution fromthe quark mass matrix M , arising from the Yukawa couplings to the Higgs field. Hence, the parameter θ shifts to θ = θ + arg det M where now θ describes the CP -violating parameter of the extended strong andEW symmetry. Here a delicate issue arises, namely that given the smallness of the total value of θ , either θ and arg det M are both tiny or they cancel each other at the level that satisfies the experimental boundon the nEDM value. This is referred to as the “strong CP problem”. We assume that a near cancellationof θ and arg det M is extremely difficult to occur, thus, we take θ to be small and we perturbatively expandaround it.At low momentum transfer, the nucleon effective Lagrangian gives rise to the expression for the nEDM [1]: | (cid:126)d N | = θ lim Q → | F ( Q ) | m N , (5)at leading order in θ , where m N denotes the mass of the neutron, Q = − q the four-momentum transfer inEuclidean space ( q = p f − p i ) and F ( Q ) is the CP -odd neutron form factor. We can, therefore, calculate theelectric dipole moment by extracting the zero momentum transfer limit of the CP -odd form factor. This isthe framework on which our work is based on. In practice, it is impossible to extract F ( Q ) at Q = 0 dueto the fact that the CP -violating nucleon matrix element, decomposes to Q k F ( Q ) ( k =1 , ,
3) and not justto F ( Q ). Hence, a direct extraction of F (0) is prohibited by the theory. We then need to parametrize the Q -dependence of F ( Q ) and then take the limit Q = 0.Since we expect that the θ -parameter is very small we can expand the exponential of the topological termin powers of θ . This enables us to integrate out (over space-time) the topological charge density, which givesthe total topological charge Q : e iθ (cid:82) d x q ( x ) ≡ e iθ Q = 1 + iθ Q + O ( θ ) ; Q = (cid:90) d x q ( x ) . (6)A consequence of this expansion is that the value of F ( Q ) depends on the nucleon two- and three-point functions correlated with the topological charge Q , namely with (cid:104) J N ( (cid:126)p f , t f ) J N ( (cid:126)p i , t i ) Q(cid:105) as well as (cid:104) J N ( (cid:126)p f , t f ) J em µ ( (cid:126)q, t ) J N ( (cid:126)p i , t i ) Q(cid:105) , respectively. While the computation of nucleon two- and three-point func-tions within lattice QCD in the absence of the parity- violating term is well-known and rather precise [26–28],when correlating with Q , introduces large statistical fluctuations. The topological charge, if sampled ade-quately, has a gaussian distribution centered at zero with a width that increases with the volume of the latticeand would by itself average to zero. If there is no correlation, or mild correlation, between the correlationfunctions and the topological charge then no signal can be obtained for F ( Q ). Different definitions of thetopological charge may be more suitable to pick up the correlations with the nucleon two- and three-pointfunctions [29, 30].Recently, a comparison of lattice QCD determinations of the topological susceptibility using a gluonic def-inition with gradient-flow, cooling and other equivalent smoothing schemes, as well as, a fermionic definitionusing spectral projectors [31–34] has been carried out [29, 30]. Although both these gluonic and fermionicdefinitions lead to compatible results in the continuum limit, the gluonic definitions are much more affectedby cut-off effects. For the spectral projectors one can choose a value of the spectral cutoff that would com-pletely eliminate discretization effects in the topological susceptibility. Therefore, using spectral projectorsis preferable since the discretization error on the topological susceptibility as well as the statistical error aresuppressed. Given the statistical superiority of spectral projectors in calculation of quantities that dependon the topological charge, one needs to consider the fermionic definition as compared to the commonly usedgluonic definition. One has to bear in mind that the computational cost of spectral projectors is orders2f magnitude larger than that of the gluonic. However, for quantities calculated at a given lattice spacingwithout a continuum extrapolation, the fermionic definition can be the only option to avoid large cut-offeffects. Since in this work we investigate the nEDM at physical pion mass we are, at the moment, restrictedto use one N f = 2 + 1 + 1 ensemble [35], and therefore a continuum extrapolation is not possible. In thefuture, we plan to include two more ensembles with smaller lattice spacings in order to be able to take thecontinuum limit.An alternative way to compute the nEDM is to use an expression from chiral perturbation theory toextrapolate the nEDM obtained at larger pion masses to the physical point. However, such an extrapolationin the pion mass carries its own uncontrolled systematic errors. Current values of F ( m π ) computed for m π >
135 MeV within lattice QCD, as can be seen in Fig. 10, are non-zero only within a few standarddeviations or even consistent with zero within the given statistical accuracy. Furthermore, F decreases asthe pion mass decreases vanishing at the chiral limit. Thus, the smaller the pion mass the more difficult itis to compute F ( Q ) with a statistical significance that would exclude a zero value. Therefore, performinga chiral extrapolation from such data can be difficult. In a recent study [36], the authors performed a chiralfit and obtained a statistically significant non-zero value for F (0) at pion mass m π = 135 MeV, concludingthat there is a signal for CP violation induced by the θ -term. However, large uncertainties on the few datapoints used in the fit, cast in doubt the reliability of the result. It is, therefore, very important to perform afirst-principles study directly at the physical pion mass.This paper is organized as follows: In Section 2 we explain our methodology, which enables us to extract F ( Q ) from two and three-point correlation functions. Subsequently, in Section 3 we provide the details ofthe lattice QCD setup and the calculation, including the expressions of the correlation functions in terms ofthe topological charge and discussion on the topological properties (Section 4) of the ensemble used in thiswork. In Section 5 we present our results for the nucleon mixing angle and the CP -odd form factor F ( Q ) aswe take the Q = 0 limit. We compare with other studies in Section 6 and finally, in Section 7 we conclude. The precise computation of the nEDM from first principles is one of the long-standing challenges and an activetopic of research within lattice QCD. The lattice QCD formulation provides an ideal framework to accessnon-perturbatively the nEDM. The first pioneering attempt was reported nearly three decades ago [37], andwas based on the introduction of an external electric field and the measurement of the shift in the associatedenergy. This approach of extracting the nEDM within lattice QCD, however, breaks unitary since it requiresthe introduction of an external electric field. Subsequently, two new approaches have been proposed. Themost commonly used method involves the calculation of the CP -odd F ( Q ) form factor by treating the θ -parameter as a small perturbation [36, 38–42]. Another approach is to extract the CP -odd F ( Q ) formfactor by simulating the theory with an imaginary θ [43, 44] to avoid the sign issue that a real θ introduces.This can be achieved either by using a field theoretical definition of the topological charge density or byreplacing the topological charge operator with the flavour-singlet pseudoscalar density employing the axialchiral Ward identities [38, 39]. Although this approach provides a well-defined framework, it requires theproduction of new ensembles at various values of θ and analytic continuation. The cost of simulations ofensembles at the physical point for several values of θ and different values of lattice spacing is currentlyprohibitively high. Therefore, in this work we use the former approach, keeping θ real and expanding inpowers of θ . 3llowing for the CP -violating θ term, the matrix element of the electromagnetic current (cid:104) N ( p (cid:48) , s (cid:48) ) |J µe.m. | N ( p, s ) (cid:105) θ = ¯ u θN ( p (cid:48) , s (cid:48) ) [Γ µ ( q )] u θN ( p, s ) , (7)can be decomposed in four form factors as followsΓ µ ( q ) = F ( Q ) γ µ + (cid:0) F ( Q ) + iγ F ( Q ) (cid:1) iσ µν q ν m θN + F A ( Q ) ( /qq µ − q γ µ ) γ m θ, N . (8)The electromagnetic current is given by J µe.m. = (cid:80) f e f ¯ ψ f γ µ ψ f , where e f is the electric charge of the quarkfield ψ f , ¯ u θN ( p (cid:48) , s (cid:48) ) is the nucleon spinor in the presence of the θ -term and q = p (cid:48) − p is the four-momentumtransfer. In the above expression F ( Q ) and F ( Q ) are the Dirac and Pauli electromagnetic form factors, F ( Q ) is the CP -odd form factor and F A ( Q ) is the anapole form factor, that vanishes for C -preservingactions. They are all expressed as function of the Euclidean four-momentum transfer squared Q . Thepresence of (cid:104) (cid:105) θ in Eq. (7) indicates that it is the state with the action that includes the θ -term. Once thedependence of the form factors on the transferred momentum is known, one can extract the electric dipolemoment from F ( Q ) in the limit Q → F ( Q ) and F ( Q ) form factors occursif these quantities are not carefully defined. In particular, the contribution to the J µe.m. matrix elementscoming from F ( Q ) transforms as an axial 4-vector as expected, only if the spinor u θN ( p, s ) appearing inEq. (7) transforms as a regular Dirac spinor under parity, i.e. only if u θN ( p (cid:48) = ( p , − (cid:126)p )) = γ u θN ( p ). Thisholds if spinors satisfy the Dirac equations: (cid:0) i/p + m θN (cid:1) u θN ( p, s ) = 0 , ¯ u θN ( p, s ) (cid:0) i/p + m θN (cid:1) = 0 , (9)where the real-valued m θN is the mass of the nucleon in the θ (cid:54) = 0 vacuum.On the lattice, the above matrix elements can be extracted from the Euclidean three-point function givenby G µ, ( θ )3 pt ( (cid:126)p f , (cid:126)q, t f , t ins ) ≡ (cid:104) J N ( (cid:126)p f , t f ) |J µe.m. ( (cid:126)q, t ins ) | ¯ J N ( (cid:126)p i , t i ) (cid:105) θ , (10)where J N ( (cid:126)p f , t f ), ¯ J N ( (cid:126)p i , t i ) are the nucleon interpolating operators that respectively create a nucleon at time t i (source) with momentum (cid:126)p i and annihilate it at time t f (sink) and momentum (cid:126)p f .Inserting unity as a complete set of energy and momentum eigenstates in Eq.(10) leads, in the largeEuclidean time limit, to the ground state contribution given as G µ, ( θ )3 pt ( (cid:126)p f , (cid:126)q, t f , t ins , t i ) = e − E fN ( t f − t ins ) e − E iN ( t ins − t i ) (cid:88) s,s (cid:48) (cid:104) J N | N ( p f , s (cid:48) ) (cid:105) θ (cid:104) N ( p f , s (cid:48) ) |J µe.m. | N ( p i , s ) (cid:105) θ (cid:104) N ( p i , s ) | ¯ J N (cid:105) θ , (11)with E iN ≡ E N ( (cid:126)p i ) = (cid:113) (cid:126)p i + ( m θN ) and E fN ≡ E N ( (cid:126)p f ) = (cid:113) (cid:126)p f + ( m θN ) . Excited states contribution aresuppressed in the above expressions. The overlap between the interpolating operators and the nucleon stateof a given momentum can be expressed as (cid:104) J N | N ( (cid:126)p, s ) (cid:105) θ = Z θN ˜ u θN ( (cid:126)p, s ) , (cid:104) N ( (cid:126)p, s ) | ¯ J N (cid:105) θ = ( Z θN ) ∗ ¯˜ u θN ( (cid:126)p, s ) , (12)where the spinors ˜ u θN ,¯˜ u θN satisfy the following Dirac equations (cid:16) i/p + m θN e − iα N ( θ ) γ (cid:17) ˜ u θN ( (cid:126)p, s ) = 0 , ¯˜ u θN ( (cid:126)p, s ) (cid:16) i/p + m θN e − iα N ( θ ) γ (cid:17) = 0 . (13)4he imaginary phase in the mass terms that arises due to the CP -violation induced by the θ -term, isparametrized by the so called mixing angle α N ( θ ). These spinors are related to the u θN , ¯ u θN that are well-behaved under parity by an axial rotation, namely˜ u θN = e iα N γ u θN , ¯˜ u θN = ¯ u θN e iα N γ . (14)This can be verified by substituting Eq. (14) in Eq. (13) and noticing that one recovers Eq. (9). With thisin mind, one can rewrite Eq. (11) as G µ, ( θ )3 pt ( (cid:126)p f , (cid:126)q, t f , t ins , t i ) (cid:39) | Z θN | e − E fN ( t f − t ins ) e − E iN ( t ins − t i ) e iα N γ (cid:32) − i/p f + m θN E fN (cid:33) Γ µ ( q ) (cid:32) − i/p i + m θN E iN (cid:33) e iα N γ , (15)where we have used the summation property of the spinors.For small θ , the quantities in the r.h.s of Eq.(15) can be safely replaced by their leading-order terms, asfollows: m θN (cid:39) m N + O ( θ ) , Z θN (cid:39) Z N + O ( θ ) ,α θN (cid:39) α (1) N θ + O ( θ ) , F ( Q ) (cid:39) F (1)3 ( Q ) θ + O ( θ ) , (16)while higher order contributions can be neglected, therefore in the following sections we will simply refer tothe mixing angle and the CP -odd form factor as α N and F ( Q ) correspondingly. The expectation values (cid:104) . . . (cid:105) θ can also be expanded by using Eq.(6). This gives (cid:104)O(cid:105) θ = 1 Z θ (cid:90) [dU][d ¯ ψ ][d ψ ] O e − S QCD + iθ Q (cid:39) (cid:104)O(cid:105) θ =0 + iθ (cid:104)OQ(cid:105) θ =0 + O ( θ ) . (17)The left hand side (l.h.s.) of Eq.(15) is now written in terms of expectation values of states with θ = 0and can be computed using gauge configurations extracted with the standard CP -even QCD action. Thethree-point function reads: G µ, ( θ )3 pt ( (cid:126)p f , (cid:126)q, t f , t ins ) = G µ, (0)3 pt ( (cid:126)p f , (cid:126)q, t f , t ins ) + iθG µ, (0)3 pt, Q ( (cid:126)p f , (cid:126)q, t f , t ins ) + · · · , (18)with G µ, (0)3 pt ( (cid:126)p f , (cid:126)q, t f , t ins ) ≡ (cid:104) J N ( (cid:126)p f , t f ) J µe.m. ( (cid:126)q, t ins ) ¯ J N ( (cid:126)p i , t i ) (cid:105) , (19) G µ, (0)3 pt, Q ( (cid:126)p f , (cid:126)q, t f , t ins ) ≡ (cid:104) J N ( (cid:126)p f , t f ) J µe.m. ( (cid:126)q, t ins ) Q ¯ J N ( (cid:126)p i , t i ) (cid:105) . (20)One can, thus, divide the right hand side (r.h.s.) of Eq. (15) into a CP -even and a CP -odd part andrelate them respectively to the three-point functions of Eq. (19) and Eq. (20). The unknown normalizationcoefficients Z N can be canceled by taking an appropriate ratio with the two-point function, while the nucleonmixing angle α N can be measured by using the relation G ( θ )2 pt ( (cid:126)p f , t f ) ≡ (cid:104) J N ( (cid:126)p f , t f ) ¯ J N ( (cid:126)p i , t i ) (cid:105) θ = | Z θN | e − E N ( t f − t i ) − i/p f + m θN e i α θN γ E N , (21)and expanding in powers of θ . The result is that the mixing angle α N can be determined from the computationof the following two-point functions G (0)2 pt ( (cid:126)p f , t f ) ≡ (cid:104) J N ( (cid:126)p f , t f ) ¯ J N ( (cid:126)p i , t i ) (cid:105) , (22) G (0)2 pt, Q ( (cid:126)p f , t f ) ≡ (cid:104) J N ( (cid:126)p f , t f ) Q ¯ J N ( (cid:126)p i , t i ) (cid:105) . (23)5 Lattice setup
We use one gauge ensemble of twisted mass fermions produced with 2 degenerate light flavours, a strangeand a charm quark ( N F = 2 + 1 + 1), with all quark masses tuned close to their physical values. We use theIwasaki improved gauge action [46], given by S G = β (cid:88) x c (cid:88) µ,ν =1 µ<ν (cid:2) − Re Tr (cid:0) U × x,µν (cid:1)(cid:3) + c (cid:88) µ,ν =1 µ (cid:54) = ν (cid:2) − Re Tr (cid:0) U × x,µν (cid:1)(cid:3) , (24)where β = 6 /g , U × the plaquette and U × the rectangular Wilson loops. The Symanzik coefficients areset to c = 3 .
648 and c = (1 − c ) /
8. The fermionic sector is implemented using the twisted mass formulationof lattice QCD [47, 48], which for the degenerate light quark double up and down has the form S lF = a (cid:88) x ¯ χ ( l ) ( x ) (cid:18) D W [ U ] + i c SW σ µν F µν [ U ] + m ,l + iµ l γ τ (cid:19) χ ( l ) ( x ) . (25)In the equation above χ ( l ) is the field representing the light quarks doublet, expressed in the twisted basis, m ,l and µ l are respectively the untwisted and twisted mass parameters, τ is the third Pauli matrix actingin flavor space and D W is the massless Wilson-Dirac operator. The clover term ∝ σ µν F µν is included inthe action to suppress cut-off effects reducing the difference between the mass of the charged and neutralpions [35]. The strange and charm quarks are included as a non-degenerate twisted doublet χ ( h ) = ( s, c ) t ,with the action [49] S hF = a (cid:88) x ¯ χ ( h ) ( x ) (cid:18) D W [ U ] + i c SW σ µν F µν [ U ] + m ,h − µ δ τ + iµ σ γ τ (cid:19) χ ( h ) ( x ) , (26)where m ,h is the bare untwisted quark mass for the heavy doublet, µ δ the bare twisted mass along the τ direction and µ σ the mass splitting in the τ direction.We tune the partial conserved axial current (PCAC) mass to zero in order to achieve maximal twist. Thisensures automatic O ( a ) improvement for the expectation values of the observables of interest [50]. Moredetails about the generation of this ensemble can be found in Ref. [35]. The lattice size is 64 × a = 0 . m π = 139(1)MeVand Lm π = 3 .
62. We will refer to this ensemble as cB211.72.64. For the analysis, we used 750 gaugeconfigurations, separated by 4 trajectories each. The parameters of the simulation are summarized in Table 1. β = 1 . , c SW = 1 . , a = 0 . × m π = 139(1)MeV , m π L = 3 . L = 5 . m N = 940(2)MeVTable 1: Simulation parameters for the cB211.72.64 ensemble [28, 35] used in this work. For the computation of the nucleon two- and three-point functions we employ the standard proton interpo-lating field, namely J N ( x ) = (cid:15) abc (cid:2) u a,T ( x ) C γ d b ( x ) (cid:3) u c ( x ) , (27)6here u ( x ) and d ( x ) are up and down quark fields in the physical base, and C = iγ γ is the charge conjugationmatrix. Since up and down quarks are degenerate in our formulation, the proton and neutron are degenerate.In order to improve the overlap with the neutron ground state and suppress excited states, we use Gaussiansmeared quark fields [51, 52], with 125 smearing steps and parameter α G = 0 .
2. The smearing is tuned inorder to approximately reproduce a mean square radius for the neutron of 0.5 fm. We apply 50 steps ofAPE-smearing [53] on the gauge links that enter the smearing operator with α APE = 0 . J µe.m. ( x ) we use the symmetrized lattice conserved vector current definedas j µf ( x ) = 14 (cid:2) ¯ ψ f ( x + ˆ µ ) U † µ ( x )(1 + γ µ ) ψ f ( x ) − ¯ ψ f ( x ) U µ ( x )(1 − γ µ ) ψ f ( x + ˆ µ )+ ¯ ψ f ( x ) U † µ ( x − ˆ µ )(1 + γ µ ) ψ f ( x − ˆ µ ) − ¯ ψ f ( x − ˆ µ ) U µ ( x − ˆ µ )(1 − γ µ ) ψ f ( x ) (cid:3) , (28)which, in contrast to the local current ¯ ψ f ( x ) γ µ ψ ( x ), does not need renormalization.Figure 1: Diagrammatic representation of the two-point function (left) and connected three-point function(right).The two- and three-point functions are then given by G pt (Γ , (cid:126)p f , t f , t i ) ≡ (cid:88) (cid:126)y Tr (cid:2) Γ (cid:104) J N ( (cid:126)y, t f ) ¯ J N ( (cid:126)x, t i ) (cid:105) (cid:3) e − (cid:126)p f ( (cid:126)y − (cid:126)x ) , (29) G pt, Q ( γ , (cid:126)p f , t f , t i ) ≡ (cid:88) (cid:126)y Tr (cid:104) γ (cid:104) J N ( (cid:126)y, t f ) Q ¯ J N ( (cid:126)x, t i ) (cid:105) (cid:105) e − (cid:126)p f ( (cid:126)y − (cid:126)x ) , (30) G µ pt (Γ k , (cid:126)q, (cid:126)p f , t f , t ins , t i ) ≡ (cid:88) (cid:126)y,(cid:126)z Tr (cid:2) Γ k (cid:104) J N ( (cid:126)y, t f ) J µe.m. ( z, t ins ) ¯ J N ( (cid:126)x, t i ) (cid:105) (cid:3) e − (cid:126)p f ( (cid:126)y − (cid:126)x ) e (cid:126)q ( (cid:126)z − (cid:126)x ) , (31) G µ pt, Q (Γ k , (cid:126)q, (cid:126)p f , t f , t ins , t i ) ≡ (cid:88) (cid:126)y,(cid:126)z Tr (cid:2) Γ k (cid:104) J N ( (cid:126)y, t f ) J µe.m. ( z, t ins ) Q ¯ J N ( (cid:126)x, t i ) (cid:105) (cid:3) e − (cid:126)p f ( (cid:126)y − (cid:126)x ) e (cid:126)q ( (cid:126)z − (cid:126)x ) , (32)where the projectors Γ and Γ k are given byΓ = 14 ( + γ ) , Γ k = i Γ γ γ k . (33)For the computation of the three-point functions, we consider only the connected contribution as shownin Fig. 1. We use the standard method of sequential inversions through the sink, taking the final momentum (cid:126)p f = (cid:126)
0. The time slice of the sink relative to the source (sink-source time separation) is kept fixed to7 f − t i = 12 a . From our previous investigation using an ensemble of N f = 2 + 1 + 1 twisted mass fermionswith pion mass m π ≈
370 MeV we showed that such a time separation is sufficient to suppress excited statecontributions to the accuracy of the present study [54]. Larger values of the sink-source time separation, evenif desirable for a better suppression of the excited states, lead to large statistical uncertainties, that for thecurrent investigation would require a prohibitively high statistics. We compute the three-point functions using54 randomly distributed source positions per configuration, that can be treated as statistically independentmeasures, leading to a total statistics of ∼
40k data. The two-point functions entering the F determinationboth explicitly, appearing in the ratio with the three-point function (see Eq. 40 in Section 5.2), and implicitlythrough the mixing angle α N . In the former case, we employ the same 54 source positions used for the three-point functions, while for the computation of the mixing angle, we had higher statistics available, namely200 source positions per configuration. A summary of the statistics used is given in Table 2.Correlation functions N src N cnfs N tot α N G pt
200 750 150000 F G pt , G pt
54 750 40500Table 2: Statistics employed for the evaluation of α N and F . As already mentioned, for the computation of the three- and two- point functions in Eqs. (30),(32) oneneeds the topological charge. In the continuum, the topological charge is defined as the integral over thefour-dimensional volume of the topological charge density given in Eq. (4), namely Q = 132 π (cid:90) d x (cid:15) µνρσ Tr [ F µν ( x ) F ρσ ( x )] . (34)The discrete counterpart of the above quantity can be obtained by replacing the gluonic field tensor with alattice operator that reproduces the correct continuum limit. The choice is not unique, and operators withbetter finite-size effects can be obtained by using O ( a )-improved discretizations of F µν . In particular, one ofthe definitions of Q we used in this work, is the symmetric or ’clover’ definition, firstly introduced in Ref. [55].It has the form Q L = 132 π (cid:88) x (cid:15) µνρσ Tr [ C µν ( x ) C ρσ ( x )] , (35)and uses a discretization of the gauge strength tensor in terms of a ”clover leaf” path C µν , made by the sumof the plaquettes P µν ( x ) centered in x and with all the possible orientations in the µν -plane, i.e. C µν ( x ) = 14 Im [ P µν ( x ) + P ν, − µ ( x ) + P − µ, − ν ( x ) + P − ν,µ ( x )] . (36)This operator is even under parity transformations and exhibits O ( a ) discretization effects. We use the gra-dient flow [56] in order to suppress the UV fluctuations of the gauge field defining the topological charge. Thesmoothing action employed in the flow equation is the standard Wilson action. The elementary integrationstep is (cid:15) = 0 .
01 and the topological charge is computed on the smoothed fields at multiples of ∆ τ flow = 0 . τ flow searching for a plateau region. Accordingly to Ref. [56] we expect that this happens for a √ τ flow ∼ O (0 . D † W D W , by employing the relation Q = λ i 65 MeV.In the rest of the paper, we will refer to the definition of Eq. (35) as “gluonic” or “field theoretic” definitionof the topological charge, while the one defined by Eq. (38) will be referred to as the “fermionic” or “spectralprojectors” definition. Before presenting results on the CP -odd form factors, we investigate the topological properties of thecB211.72.64 gauge ensemble. In the left panel of Fig. 2 we show the Monte Carlo (MC) history of thetopological charge using the two definitions. One can observe that the autocorrelation time of Q is notcritical, i.e. no topological freezing occurs at this lattice spacing, and the configurations explore several We would like to thank M. Constantinou for providing this value. τ flow = 3 . 5, and M thr = 64 . 98 MeV but similarobservations hold also for other values of the smoothing and cut-off scales. traj − − Q field theo. spectral proj. − − 20 0 20 40 Q . . . . Q field theo.spectral proj. Figure 2: Monte Carlo history of Q (left panel), and resulting histogram (right panel) computed on 750configurations (selected every fourth trajectory), using the gluonic τ flow = 3 . M thr = 64 . 98 MeV (red) definitions. τ flow . . . . . χ t o p [ f m − ] field theo. 20 30 40 50 60 M thr [MeV] . . . . . . . χ t o p [ f m − ] spectral proj. Figure 3: Topological susceptibility with the two definitions of Q using the gluonic definition (left) as afunction of the time flow and the fermionic (right).In the right panel of Fig. 2, histograms of Q are compared. Both definitions lead to a Gaussian-likedistribution with a mean value of the topological charge compatible with zero within errors with, however,different widths. In particular, the value of the topological charge extracted via spectral projectors shows lessfluctuations and, thus, exhibits a narrow distribution. This impacts the values of the topological susceptibilitydefined as (cid:104)Q (cid:105) /V that differs between the two definitions, as can be seen from Fig. 3. This mismatch is10ctually expected, due to the fact that the spectral projectors definition shows milder cut-off effects withrespect to (cid:104)Q (cid:105) /V computed via the field theoretical definition. This was already observed in Ref. [29].Another observation that emerges from Fig. 3 is that, while the susceptibility computed using the fieldtheoretical definition is essentially flat as a function of the time flow, the one computed using spectralprojectors depends strongly on the threshold M thr . The absence of a plateau for the window of eigenvaluescomputed makes it harder to properly justify a good choice of M thr . In principle, there is no reason to expectto obtain a plateau for a large value of the renormalized cut-off. Actually, the results from Ref. [30] suggestthat for the particular lattice spacing we cannot observe a plateau for cut-offs up to 180 MeV.However, the fact that we do not observe a plateau is not a problem di per se . Indeed, the particularchoice of M thr becomes irrelevant once one takes the continuum limit, as long as the renormalized value ofthe threshold in physical units is kept fixed. Since in this work we use one gauge ensemble, the question thatarises is what are the finite-size effects that result from taking a threshold that is not in the plateau region. . 000 0 . 002 0 . 004 0 . 006 0 . a [fm ] . . . . . . . r χ field theo.160MeV120MeV 80MeV60MeV Figure 4: Continuum extrapolation of topological susceptibility taken from Ref. [29], with the adjoint ofpoints at M thr = 60 MeV. Lattice spacing used in that work corresponds to a/r = 0 . M thr are negligible in respect to the ones that affect the field theoreticaldefinition.To give a qualitative answer we use the results obtained in the work of Ref. [29]. In Fig. 4, we show thecontinuum extrapolation of the topological susceptibility for three values of the threshold used in Ref. [29],plus a fourth value, at M thr = 60MeV. What can be seen is that, even if the slope of the continuum limitincreases when we decrease the threshold, it is still milder than the slope of the susceptibility computed viathe gluonic definition. Thus, the error in taking the value at finite a instead of its continuum extrapolation,is negligible with respect to the one coming from cut-off effects of the gluonic definition. Therefore, based onthis observation, we use the maximum value of the threshold accessible with the current number of eigenvaluescomputed, i.e. M thr = 64 . 98 MeV. We can reasonably expect that similar considerations also hold for thequantities entering the determination of the nEDM, i.e. the two- and three-point function correlated with Q . Their dependence on the smoothing scale ( τ flow and cut-off M thr ) is investigated in the following section.11 Results For the determination of the CP -odd form factor F ( Q ), one requires the knowledge of the mixing angle α N . We extract it from the following ratio of two-point functions at zero-momentum α N = lim t f →∞ G pt, Q (Γ ,(cid:126) , t f ) G pt (Γ ,(cid:126) , t f ) . (39)We take t i = 0 and thus suppress the dependence on t i . We seek for an interval for which the ratio becomestime-independent as a function of t f (plateau region). The time evolution of the ratio of Eq. (39) is illustratedin Fig. 5 as a function of t f /a . We use both the gluonic (left panel) and the fermionic (right panel) definitionsfor the topological charge that enters in the computation of α N . In particular, the one in the left panel iscomputed at timeflow τ flow = 3 . 5, while for the fermionic one, we used a cut-off of M thr = 64 . 98 MeV. t f /a . . . . . . α N field theo. 5 10 15 20 t f /a . . . . . . α N spectral proj. Figure 5: Value of the ratio in Eq. (39), as a function of t f /a , using the two definitions of the topologicalcharge, gluonic definition of Eq. (35) (left) and fermionic definition of Eq. (38) (right). With the blue (left)and red (right) bands we show the result of a constant fit within the plateau. With the grey band we showthe corresponding fits when using for the fit a constant plus an exponential term which takes into accountthe first excited state.We identify a plateau region in the range t f /a ∈ [9 , 18] and t f /a ∈ [10 , 18] for the gluonic and fermionicdefinitions, respectively. Fitting to a constant, we find α N = 0 . α N = 0 . , 12] and the final time slice in the interval [17 , Q , it is negligible if compared to the statistical uncertainty. To account for a residualtime-dependence we consider the contribution of the first excited state too. As shown in Fig. 5, both Ans¨atzegive compatible results and further validate the choice of the plateau region.In Fig. 6, we show the dependence of α N on the smoothing scale of the gluonic definition of the topologicalcharge. Values reported are extracted using a constant fit in the plateau region. As can be seen, from τ flow = 3,our resulting value of α N does not depend on τ flow . In the same figure we also show the dependence of α N on the threshold M thr computed using the spectral projectors for the topological charge and a constant fit.12his shows a stronger dependence on M thr . Such a behaviour is reminiscent of the one exhibited by thetopological charge as discussed in Section 4.1. τ flow . . . . . α N field theo. 20 30 40 50 60 M thr [MeV] . . . . . . . α N spectral proj. Figure 6: Dependence of the nucleon mixing angle α N on the smoothing scale τ flow and cut-off M thr . Left,using the gluonic definition of the topological charge and right, using spectral projectors. In both cases α N is extracted using a constant fit within the plateau region of the ratio of Eq. (39). CP -odd form factor F ( Q ) In order to extract the CP -odd form factor F ( Q ), we construct an appropriate ratio of the relevant three-point function given in Eq.(32) and a combination of two-point functions so that unknown overlaps and theexponential time dependence cancel in the large time limit of t f and t ins . This ratio is given byΠ µk pt, Q ( (cid:126)q ) ≡ lim t f ,t ins →∞ G µ pt, Q (Γ k , (cid:126)q, t f , t ins ) G pt (Γ ,(cid:126) , t f ) (cid:115) G pt (Γ , (cid:126)q, t f − t ins ) G pt (Γ ,(cid:126) , t ins ) G pt (Γ ,(cid:126) , t f ) G pt (Γ ,(cid:126) , t f − t ins ) G pt (Γ , (cid:126)q, t ins ) G pt (Γ , (cid:126)q, t f ) , (40)where (cid:126)p f = (cid:126) 0. The form factor F ( Q ) is then extracted fromΠ k pt, Q ( (cid:126)q ) = iq k C m N (cid:18) α N G E ( Q ) − F ( Q )2 m N ( E N + m N ) (cid:19) , (41)where E N is the initial energy of the nucleon, C = (cid:112) (2 m N ) / ( E N ( E N + m N )) is a kinematic factor and G E ( Q ) = F ( Q ) + ( q / (2 m N )) F ( Q ) is the electric Sachs form factor. F ( Q ) and F ( Q ) are the Pauliand Dirac form factors. We note that Eq. (41) derives from Eq. (55) of Ref. [54] with the correction given inRef. [45]. The neutron electric form factor G E ( Q ) is extracted from Π pt (see Eq. (A4) of Ref. [28])Π pt = C E N + m N m N G E ( Q ) . (42)We note that Eq. (41) is not the only way to extract the F ( Q ) form factor. Similar relations can be writtenthat express the Π kk pt, Q and Π j (cid:54) = k pt, Q components as linear combinations of F , F and F . We verified thatthey lead to worse signal-to-noise ratios, and thus we opt to use Eq. (41) in this work.13 − ( t ins − t f ) /a − . − . . . . -iΠ k /q k Q = 0 . 058 GeV − − ( t ins − t f ) /a − . − . . . . -iΠ k /q k Q = 0 . 058 GeV − − ( t ins − t f ) /a − . − . . . . -iΠ k /q k Q = 0 . 113 GeV − − ( t ins − t f ) /a − . − . . . . -iΠ k /q k Q = 0 . 113 GeV − − ( t ins − t f ) /a − . − . . . . -iΠ k /q k Q = 0 . 168 GeV − − ( t ins − t f ) /a − . − . . . . -iΠ k /q k Q = 0 . 168 GeV Figure 7: Ratio of Eq. (40) as a function of insertion time t ins at fixed sink-source time separation t f = 12 a .The three smallest values of the momentum transfer squared are shown. In the left column, we show resultsusing the gluonic definition of Q and ( τ flow = 3 . Q and ( M thr = 64 . 98 MeV). The bands are the result ofa constant fit in the plateau region excluding symmetrically 3 and 4 time slices for the gluonic (left) andfermionic (right) definition of Q , respectively.In Fig. 7, we show the ratio defined in Eq. (40) as a function of the insertion time t ins at fixed sink-source time separation t f = 12 a , for the smallest three non-zero values of the momentum transfer squared,i.e. Q = 0 . 056 GeV , Q = 0 . 111 GeV and Q = 0 . 164 GeV . The results shown for Π k pt, Q are aver-14ged among momenta with non-zero k -component in all k -directions. Despite the large statistics employed( ∼ − t fit , t fit ] and vary the fit ranges, taking t fit = 2 , , 4. The resultingvalues are all compatible within our current statistical accuracy. It is worth to be noting that the resultsextracted using the fermionic definition of the topological charge show a significant error reduction withrespect to the gluonic counterpart. Nevertheless, a zero value cannot be excluded for all three momentumtransfers. This is made even more evident in Fig. 8, where the values of F at different Q are reported. Inorder to extract F (0) and thus d θN , we extrapolate to Q = 0 by considering the weighted average of thevalues at the three smallest Q values. Other fit forms, such as the dipole Ansatz and the z -expansion [62],or even the momentum elimination technique [63], are not viable with this level of uncertainty. Our finalresults are field theoretical or gluonic definition | d θN | = 0 . θ e · fm , (43)fermionic definition via spectral projectors | d θN | = 0 . θ e · fm . (44) . 00 0 . 05 0 . 10 0 . 15 0 . 20 0 . 25 0 . 30 0 . Q [GeV ] − . − . . . . . F m N [e · fm] field theo. 0 . 00 0 . 05 0 . 10 0 . 15 0 . 20 0 . 25 0 . 30 0 . Q [GeV ] − . − . . . . . F m N [e · fm] spectral proj. Figure 8: F ( Q ) as a function of Q , using a the field theoretical or gluonic definition of the topologicalcharge (left panel) and the fermionic definition based on spectral projectors (right panel). The blue (left)and red (right) bands represent the weighted average of the values at three smaller values of Q .As can be seen, the fermionic definition based on spectral projectors of the topological charge provides amore accurate determination than the gluonic definition. Considering the absolute error as a bound for themagnitude of the nEDM, what we have is an improvement by a factor ∼ Q , is compensated by the increased precision, leading to a large payoff in terms of the computational cost.Moreover, the choice of maximum cut-off in the number of eigenmodes, i.e. M thr = 64 . 98 MeV, employed inthe computation of the results presented, is in some sense the most conservative choice. Indeed, by looking atthe dependence of the extrapolated value of F (0) as a function of M thr , shown in the bottom row of Fig. 9,it can be seen that the mean value of the form factor does not depend on M thr and only the error increases15ith increasing M thr . This is in contrast with what is observed for the mixing angle α N , that instead changessignificantly with M thr . This is because the changes in α N are small as compared to the uncertainty of Π k pt, Q and are thus not visible in F at the current level of precision. τ flow − . . . . F m N Q = 0 . 058 GeV [e · fm] τ flow − . . . . F m N Q = 0 . 113 GeV [e · fm] τ flow − . . . . F m N Q = 0 . 168 GeV [e · fm] 10 20 30 40 50 60 M thr [MeV] − . . . . F m N Q = 0 . 058 GeV [e · fm] 10 20 30 40 50 60 M thr [MeV] − . . . . F m N Q = 0 . 113 GeV [e · fm] 10 20 30 40 50 60 M thr [MeV] − . . . . F m N Q = 0 . 168 GeV [e · fm] Figure 9: Dependence of the F (0) on the smoothing scale τ flow for the gluonic (upper row) and cut-off M thr for the fermionic (bottom row) definitions used in computation of the topological charge, for the three smallervalues of the momentum transfer squared. 16 . . . . . . . . | d θN | [e · fm] m π [ M e V ] J.Dragos et al , 2019, a = 0 . E.Shintani et al , 2005, a = 0 . F. − K.Guo et al , 2015, a = 0 . C.Alexandrou et al , 2016, a = 0 . a = 0 . Figure 10: Comparison with other lattice QCD determinations of nEDM present in literature. Values fromRefs. [40, 44, 54] (dashed error bars) are not the ones from their original papers but are taken from Table IIIof Ref. [45], where the spurious contribution coming from F ( Q ) is subtracted. See Ref. [45] for furtherdetails. In Fig. 10 we provide a comparison of our result with those of other lattice QCD studies for a similar latticespacing. We observe that our value has a statistical error that is comparable with the recent results inRef. [36] that were, however, computed at much larger values of the pion mass. Since the errors grow withdecreasing pion mass and so does the computational cost, achieving such an accuracy it is a major outcomeof this work. We note that the authors of Ref. [36] using their results at these heavy pion masses performa chiral extrapolation and find a non-zero value of | d θN | = 0 . θ e · fm at the physical point. Theaccuracy of their determination is due to the chiral extrapolation where systematic errors from using chiralexpressions cannot be determined. It is worth mentioning that their actual data have uncertainties similarto this work. Our result computed directly at the physical point does not exclude a zero value, but providea significant bound to the value of nEDM. We compute the neutron electric dipole moment using an ensemble of N f = 2 + 1 + 1 twisted mass fermionssimulated at the physical point and with a lattice spacing of a (cid:39) . 08 fm. The extraction of the CP -violating17orm factor F with the approach adopted in this work requires the calculation of the topological charge.For this reason we made use of a gluonic definition of the topological charge as well as a fermionic definitionby means of gradient flow and spectral projectors, respectively. This enables us to test what effect differentdefinitions of the topological charge may have on the nEDM. F ( Q ) cannot be extracted directly from thenucleon matrix element due to the presence of momentum appearing multiplicatively in front of the formfactor. The usual approach is to compute F ( Q ) for finite Q and then extrapolate to Q = 0 using somefitting form. We would like to highlight the following three crucial aspects of our investigation of the CP -oddform factor: • Determination of nEDM directly at the physical point : We perform the computation directly at thephysical value of the pion avoiding uncontrolled errors that a chiral extrapolation may result in. Westudy F ( Q ) as a function of the the momentum transfer Q , and show that there is no noticeabledependence on low values of Q within our current statistical accuracy. We, thus, perform a weightedaverage on the three lowest ones to extract the value of the form factors at Q = 0. • Definition of the topological charge : The approach is based on using spectral projectors to determinethe topological susceptibility yielding a statistical error that is twice smaller than when using a gluonicdefinition whether using gradient flow, link smearing, or cooling are applied [29, 30]. Furthermore, ithas been shown that the topological susceptibility calculated using the spectral projectors approachshows milder lattice artifacts than when using a gluonic definition [29]. Both these observations haveled us to investigate the approach based on spectral projectors for the definition of the topologicalcharge entering the evaluation of the nEDM since we expected that these features would also hold bysecondary quantities involving the topological charge. The necessity to reduce as much as possible thelattice cut-off effects stem from the fact that currently we do not have three ensembles with differentlattice spacings at the physical point to enable us to take the continuum limit. Therefore, while F (0) extracted using the fermionic definition is consistent with the value extracted using the gluonicdefinition, the statistical errors of the former are about half as compared to those when using thegluonic definition and the cut-off effects are expected from the above consideration to be milder. We,thus, quote as our final value of nEDM the one resulting from using spectral projectors to define thetopological charge. We find | d θN | = 0 . θ e, fm (45)which is compatible with zero.There is a subtlety that we would like to clarify at this point, namely the fact that we do not observea clear plateau for α N when using spectral projectors as a function of M thr . This is an artifact ofcarrying out the investigation at one lattice spacing. Only after a continuum extrapolation keepingthe renormalised cut-off M thr in physical units fixed one can extract physical observables that areindependent of the cut-off. The same holds for the field theoretic definition for which the gradient flowtime in physical units must be fixed as one takes the continuum limit. This is illustrated in Fig. 3for the susceptibility, where only after the continuum limit is taken there is an agreement among theapproaches.A question arises as of why when using the gradient flow with the Wilson smoothing action we obtainquantities such as the topological susceptibility and α N that have a very mild dependence on the scale τ flow leading to a fast convergence to a plateau. This can be understood within the following context.18n the semi-classical level, gradient flow smooths out small instantons corresponding to ultra-violet(UV) contributions. This happens even at small flow times and, thus, we are left with a topologicalcontent which remains unaltered until the smoother starts affecting the large instantons contributingto the topological structure of the theory. On the other hand, the same quantities extracted usingthe spectral projectors reveal a behaviour that depends much stronger on the associated cut-off. Thisis because, for a theory that breaks chiral symmetry, as we sum up eigenvalues of the squared DiracMatrix, the summation will keep increasing. This observation is fully supported by the results on thetopological susceptibility demonstrated in Ref. [29]. One may argue that the existence of a plateauenables to reliably quote a value of a quantity at a given lattice spacing; this is of course true buthas no physical importance since in the end a continuum extrapolation is needed in order to get ridof discretization effects. Nevertheless, we observe that, when using spectral projectors, F appears toexhibit a plateau versus M thr coupled with the extra benefit of reduction in the statistical uncertainty. • Alternative approaches of extracting the θ -induced nEDM. The current investigation reveals the difficultyof this particular method to deliver a statistically significant result. Practically, one needs to increasemeasurements of at least an order of magnitude to reduce the error significantly hopping to exclude azero value. Employment of techniques, such as volume clustering [24, 64] or spectral projectors as adefinition of the topological charge as done in this work, help but do not solve the problem.Another possibility is the investigation of the nEDM using configurations generated with an imaginary θ -term. Since the nEDM depends on contributions from non-trivial topological sectors, introducing adynamical θ -term one may improve the importance sampling for the nEDM signal inducing nonzeroaverage topological charge. Exploratory investigations using this method with a field theoretic definitionof the topological charge are under way. We would like to thank all the members of the Extended Twisted Mass Collaboration for a most conducivecooperation. In particular, we express our gratitude to K. Cichy for providing the data on the susceptibility.We thank M. Constantinou for providing us the value of Z S . We also thank M. D’ Elia and C. Bonannofor valuable discussions. A.A. is financially supported by the European Union’s Horizon 2020 research andinnovation programme “Tips in SCQFT” under the Marie Sk(cid:32)lodowska-Curie grant agreement No. 791122.K.H. is financially supported by the Cyprus Research and Innovation foundation under contract numberPOST-DOC/0718/0100. A.T. is a Marie Sklodowska-Curie fellow funded by the European Union’s Horizon2020 research and innovation programme under grant agreement No 765048. Results were obtained usingPiz Daint at Centro Svizzero di Calcolo Scientifico (CSCS), via the project with id s702. We acknowledgePRACE for awarding us access to Marconi100 at CINECA, Italy, where part of our work is carried outwithin the project with Id Pra20 5171. This work also used computational resources from the John vonNeumann-Institute for Computing on the Juwels system at the research center in J¨ulich, under the projectwith id CHCH02. References [1] Maxim Pospelov and Adam Ritz. “Electric dipole moments as probes of new physics”. In: Annals Phys. 318 (2005), pp. 119–169. doi : . arXiv: hep-ph/0504231 .192] V. Helaine. “The neutron Electric Dipole Moment experiment at the Paul Scherrer Institute”. In: EPJWeb Conf. 73 (2014). Ed. by Marco Battaglieri et al., p. 07003. doi : .[3] C.A. Baker et al. “An Improved experimental limit on the electric dipole moment of the neutron”. In: Phys. Rev. Lett. 97 (2006), p. 131801. doi : . arXiv: hep-ex/0602020 .[4] C.A. Baker et al. “Reply to comment on ‘An Improved experimental limit on the electric dipole momentof the neutron’”. In: Phys. Rev. Lett. 98 (2007), p. 149102. doi : .arXiv: .[5] C. Abel et al. “Measurement of the permanent electric dipole moment of the neutron”. In: Phys. Rev.Lett. doi : . arXiv: .[6] Varouzhan Baluni. “CP Violating Effects in QCD”. In: Phys. Rev. D 19 (1979), pp. 2227–2230. doi : .[7] R.J. Crewther et al. “Chiral Estimate of the Electric Dipole Moment of the Neutron in QuantumChromodynamics”. In: Phys. Lett. B 88 (1979). [Erratum: Phys.Lett.B 91, 487 (1980)], p. 123. doi : .[8] Mikhail A. Shifman, A.I. Vainshtein, and Valentin I. Zakharov. “Can Confinement Ensure Natural CPInvariance of Strong Interactions?” In: Nucl. Phys. B 166 (1980), pp. 493–506. doi : .[9] Howard J. Schnitzer. “The Soft Pion Skyrmion Lagrangian and Strong CP Violation”. In: Phys. Lett.B 139 (1984), pp. 217–222. doi : .[10] E.P. Shabalin. “THE ELECTRIC DIPOLE MOMENT OF THE NEUTRON IN A GAUGE THEORY”.In: Sov. Phys. Usp. 26 (1983), p. 297. doi : .[11] M.M. Musakhanov and Z.Z. Israilov. “THE ELECTRIC DIPOLE MOMENT OF THE NEUTRONIN THE CHIRAL BAG MODEL”. In: Phys. Lett. B 137 (1984), pp. 419–421. doi : .[12] P. Cea and G. Nardulli. “A Realistic Calculation of the Electric Dipole Moment of the Neutron Inducedby Strong CP Violation”. In: Phys. Lett. B 144 (1984), pp. 115–118. doi : .[13] Michael A. Morgan and Gerald A. Miller. “The Neutron Electric Dipole Moment in the Cloudy BagModel”. In: Phys. Lett. B 179 (1986), pp. 379–384. doi : .[14] J.A. McGovern and M.C. Birse. “The neutron electric dipole moment in chiral quark-meson models”.In: International Workshop on Quark Cluster Dynamics . 1992, pp. 279–280.[15] Antonio Pich and Eduardo de Rafael. “Strong CP violation in an effective chiral Lagrangian approach”.In: Nucl. Phys. B 367 (1991), pp. 313–333. doi : .[16] B. Borasoy. “The Electric dipole moment of the neutron in chiral perturbation theory”. In: Phys. Rev.D 61 (2000), p. 114017. doi : . arXiv: hep-ph/0004011 .[17] W.H. Hockings and U. van Kolck. “The Electric dipole form factor of the nucleon”. In: Phys. Lett. B 605 (2005), pp. 273–278. doi : . arXiv: nucl-th/0508012 .[18] Stephan Narison. “A Fresh Look into the Neutron EDM and Magnetic Susceptibility”. In: Phys. Lett.B 666 (2008), pp. 455–461. doi : . arXiv: .[19] K. Ottnad et al. “New insights into the neutron electric dipole moment”. In: Phys. Lett. B 687 (2010),pp. 42–47. doi : . arXiv: .2020] J. de Vries et al. “The Nucleon Electric Dipole Form Factor From Dimension-Six Time-Reversal Vio-lation”. In: Phys. Lett. B 695 (2011), pp. 268–274. doi : . arXiv: .[21] E. Mereghetti et al. “The Electric Dipole Form Factor of the Nucleon in Chiral Perturbation Theory toSub-leading Order”. In: Phys. Lett. B 696 (2011), pp. 97–102. doi : .arXiv: .[22] J. de Vries et al. “The Effective Chiral Lagrangian From Dimension-Six Parity and Time-ReversalViolation”. In: Annals Phys. 338 (2013), pp. 50–96. doi : 10 . 1016 / j . aop . 2013 . 05 . 022 . arXiv: .[23] Feng-Kun Guo and Ulf-G. Meissner. “Baryon electric dipole moments from strong CP violation”. In: JHEP 12 (2012), p. 097. doi : . arXiv: .[24] Sergey Syritsyn, Taku Izubuchi, and Hiroshi Ohki. “Calculation of Nucleon Electric Dipole MomentsInduced by Quark Chromo-Electric Dipole Moments and the QCD θ -term”. In: PoS Confinement2018(2019), p. 194. doi : . arXiv: .[25] S. Syritsyn, T. Izubuchi, and H. Ohki. “Progress in the Nucleon Electric Dipole Moment Calculationsin Lattice QCD”. In: . Oct. 2018.arXiv: .[26] C. Alexandrou et al. “Complete flavor decomposition of the spin and momentum fraction of the protonusing lattice QCD simulations at physical pion mass”. In: Phys. Rev. D doi : . arXiv: .[27] C. Alexandrou et al. “The nucleon axial, tensor and scalar charges and σ -terms in lattice QCD”. In:(Sept. 2019). arXiv: .[28] C. Alexandrou et al. “Proton and neutron electromagnetic form factors from lattice QCD”. In: Phys.Rev. D doi : . arXiv: .[29] Constantia Alexandrou et al. “Topological susceptibility from twisted mass fermions using spectralprojectors and the gradient flow”. In: Phys. Rev. D doi : . arXiv: .[30] Constantia Alexandrou et al. “Comparison of topological charge definitions in Lattice QCD”. In: Eur.Phys. J. C doi : 10 . 1140 / epjc / s10052 - 020 - 7984 - 9 . arXiv: .[31] Martin Luscher. “Topological effects in QCD and the problem of short distance singularities”. In: Phys.Lett. B 593 (2004), pp. 296–301. doi : . arXiv: hep-th/0404034 .[32] Martin Luscher and Filippo Palombi. “Universality of the topological susceptibility in the SU(3) gaugetheory”. In: JHEP 09 (2010), p. 110. doi : . arXiv: .[33] Leonardo Giusti and Martin Luscher. “Chiral symmetry breaking and the Banks-Casher relation inlattice QCD with Wilson quarks”. In: JHEP 03 (2009), p. 013. doi : .arXiv: .[34] Claudio Bonanno et al. “Topology via spectral projectors with staggered fermions”. In: JHEP 10 (2019),p. 187. doi : . arXiv: .[35] Constantia Alexandrou et al. “Simulating twisted mass fermions at physical light, strange and charmquark masses”. In: Phys. Rev. D doi : . arXiv: . 2136] Jack Dragos et al. “Confirming the Existence of the strong CP Problem in Lattice QCD with theGradient Flow”. In: (Feb. 2019). arXiv: .[37] Sinya Aoki and Andreas Gocksch. “The Neutron Electric Dipole Moment in Lattice QCD”. In: Phys.Rev. Lett. 63 (1989). [Erratum: Phys.Rev.Lett. 65, 1172 (1990)], p. 1125. doi : .[38] D. Guadagnoli et al. “Neutron electric dipole moment on the lattice: A Theoretical reappraisal”. In: JHEP 04 (2003), p. 019. doi : . arXiv: hep-lat/0210044 .[39] P. Faccioli, D. Guadagnoli, and S. Simula. “The Neutron electric dipole moment in the instantonvacuum: Quenched versus unquenched simulations”. In: Phys. Rev. D 70 (2004), p. 074017. doi : . arXiv: hep-ph/0406336 .[40] E. Shintani et al. “Neutron electric dipole moment from lattice QCD”. In: Phys. Rev. D 72 (2005),p. 014504. doi : . arXiv: hep-lat/0505022 .[41] Eigo Shintani et al. “Neutron and proton EDM with N f = 2 + 1 domain-wall fermion”. In: PoS LATTICE2013 (2014), p. 298. doi : .[42] Andrea Shindler, Thomas Luu, and Jordy de Vries. “Nucleon electric dipole moment with the gradientflow: The θ -term contribution”. In: Phys. Rev. D doi : . arXiv: .[43] S. Aoki et al. “The Electric dipole moment of the nucleon from simulations at imaginary vacuum angletheta”. In: (Aug. 2008). arXiv: .[44] F.-K. Guo et al. “The electric dipole moment of the neutron from 2+1 flavor lattice QCD”. In: Phys.Rev. Lett. doi : 10 . 1103 / PhysRevLett . 115 . 062001 . arXiv: .[45] M. Abramczyk et al. “Lattice calculation of electric dipole moments and form factors of the nucleon”.In: Phys. Rev. D doi : . arXiv: .[46] Y. Iwasaki. “Renormalization Group Analysis of Lattice Theories and Improved Lattice Action: Two-Dimensional Nonlinear O(N) Sigma Model”. In: Nucl. Phys. B 258 (1985), pp. 141–156. doi : .[47] Roberto Frezzotti et al. “Lattice QCD with a chirally twisted mass term”. In: JHEP 08 (2001), p. 058.arXiv: hep-lat/0101001 .[48] R. Frezzotti and G.C. Rossi. “Chirally improving Wilson fermions. 1. O(a) improvement”. In: JHEP 08 (2004), p. 007. doi : . arXiv: hep-lat/0306014 .[49] R. Frezzotti and G.C. Rossi. “Twisted mass lattice QCD with mass nondegenerate quarks”. In: Nucl.Phys. B Proc. Suppl. 128 (2004). Ed. by A.C. Kalloniatis, D.B. Leinweber, and A.G. Williams, pp. 193–202. doi : . arXiv: hep-lat/0311008 .[50] R. Frezzotti et al. “Reducing cutoff effects in maximally twisted lattice QCD close to the chiral limit”.In: JHEP 04 (2006), p. 038. doi : . arXiv: hep-lat/0503034 .[51] S. Gusken. “A Study of smearing techniques for hadron correlation functions”. In: Nucl. Phys. B Proc.Suppl. 17 (1990). Ed. by N. Cabibbo et al., pp. 361–364. doi : .[52] C. Alexandrou et al. “The Static approximation of heavy - light quark systems: A Systematic latticestudy”. In: Nucl. Phys. B 414 (1994), pp. 815–855. doi : . arXiv: hep-lat/9211042 . 2253] M. Albanese et al. “Glueball Masses and String Tension in Lattice QCD”. In: Phys. Lett. B 192 (1987),pp. 163–169. doi : .[54] C. Alexandrou et al. “Neutron electric dipole moment using N f = 2 + 1 + 1 twisted mass fermions”.In: Phys. Rev. D doi : . arXiv: .[55] P. Di Vecchia et al. “Preliminary Evidence for U(1)-A Breaking in QCD from Lattice Calculations”.In: (May 1981). Ed. by J. Julve and M. Ram´on-Medrano, pp. 426–442. doi : .[56] Martin L¨uscher. “Properties and uses of the Wilson flow in lattice QCD”. In: JHEP 08 (2010). [Erratum:JHEP 03, 092 (2014)], p. 071. doi : . arXiv: .[57] G. Martinelli et al. “A General method for nonperturbative renormalization of lattice operators”.In: Nucl. Phys. B 445 (1995), pp. 81–108. doi : 10 . 1016 / 0550 - 3213(95 ) 00126 - D . arXiv: hep -lat/9411010 .[58] M. Gockeler et al. “Nonperturbative renormalization of composite operators in lattice QCD”. In: Nucl.Phys. B 544 (1999), pp. 699–733. doi : . arXiv: hep-lat/9807044 .[59] C. Alexandrou et al. “Renormalization constants for 2-twist operators in twisted mass QCD”. In: Phys.Rev. D 83 (2011), p. 014503. doi : . arXiv: .[60] C. Alexandrou et al. “Renormalization constants of local operators for Wilson type improved fermions”.In: Phys. Rev. D 86 (2012), p. 014505. doi : 10 . 1103 / PhysRevD . 86 . 014505 . arXiv: .[61] C. Alexandrou et al. “Moments of nucleon generalized parton distributions from lattice QCD simulationsat physical pion mass”. In: Phys. Rev. D doi : .arXiv: .[62] Bhubanjyoti Bhattacharya, Richard J. Hill, and Gil Paz. “Model independent determination of the axialmass parameter in quasielastic neutrino-nucleon scattering”. In: Phys. Rev. D 84 (2011), p. 073006. doi : . arXiv: .[63] Constantia Alexandrou et al. “Model-independent determination of the nucleon charge radius fromlattice QCD”. In: Phys. Rev. D doi : . arXiv: .[64] Keh-Fei Liu, Jian Liang, and Yi-Bo Yang. “Variance Reduction and Cluster Decomposition”. In: Phys.Rev. D doi : . arXiv:1705.06358 [hep-lat]