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

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


Introduction
A nonzero measurement of the electric dipole moment (EDM) of the nucleon in the foreseeable future would be a clear signal of new physics, since the known CP-violating phase of the CKM matrix leads to EDMs that lie orders of magnitude below current experimental limits. The source of a nonzero EDM could then either be the QCDθ-term or higher-dimension CP-violating quarkgluon operators that originate in beyond-the-Standard Model (BSM) physics, or a combination of these two. To interpret an EDM signal or lack thereof, and to possibly disentangle the source (e.g.θ-term or BSM), requires a non-perturbative calculation linking the CP-violating sources to the hadronic observables.
Lattice QCD can calculate the nucleon EDM directly in terms of CP-violating operators at the quark level. Various attempts have been made in this regard [1,2,3,4,5,6,7]. However, the renormalization of CP-violating operators within a lattice (discretised) formulation of QCD is very non-trivial, and for several operators presents large difficulties in interpreting lattice results. Further, theθ term itself introduces a complex phase in the determinant of the quark matrix, which produces a sign problem and precludes the use of standard stochastic methods. Several techniques have been used to address theθ-term contribution to the EDM and attempts have been made to solve the complicated renormalization patterns of the CP-violating operators [8]. We refer to the recent review [9] for a summary.
We proposed to use the gradient flow to calculate all CP-violating source to the EDM in refs. [10,11], and presented preliminary results [12,13,14,15]. In this paper we consider thē θ-term contribution to the EDM in a perturbative manner as discussed in [11]. This is well justified considering the stringent constraints onθ set by EDM experiments. In this way we avoid the problem of a complex fermionic determinant. To define the QCDθ term we use the gradient flow. The topological charge defined in this way has a finite and well defined continuum limit [16,17,18]. It is much faster to compute than using the Ginsparg-Wilson definition and it is theoretically more robust than definitions using cooling techniques. Another problem that hinders lattice calculations of the nucleon EDM is the very poor signal-to-noise ratio. In this respect we explore a novel technique to determine the space-time region where the signal in the relevant correlation functions is maximized. A first account of this technique has already been presented [12].
Additional insights into the EDM of the nucleon [19,20,21,22,23] and nuclei [24,25] can be provided by chiral effective field theory. The CP-violating quark operators are translated to effective CP-violating hadronic operators and EDMs depend on the unknown low-energy constants (LECs) of the theory. The LECs can be estimated from dimensional analysis or, preferably, be determined from experiments and/or by lattice-QCD calculations. We use these insights from chiral calculations to understand the pion mass dependence of our results, and to connect our nucleon EDM calculations to nuclear EDMs.
The remainder of the paper is organized as follows: section 2 gives a cursory discussion of the phenomenology of the nucleon EDM, followed by an overview of the lattice details and parameters in section 3, where we discuss the general lattice strategy used and we define the basic observables, including the gradient flow. Sanity checks are given in section 4, where we compute and display the topological charge using the gradient flow. The nucleon two-point correlation function is explored in section 5, then the setup of the computation of the EDM is derived and results shown in section 6. A discussion follows the results section in section 7, where we discuss the ramifications of our results and compare our results to the literature. Finally, we conclude in section 8.
2 Phenomenology of the QCD theta term.
The discrete space-time symmetries parity (P ) and time-reversal (T ), and hence via the CP T theorem also CP symmetry, are broken in QCD by theθ term. In the case of two quark flavors the QCD Lagrangian in Minkowski space is given by where q = (u , d) T denotes the quark doublet containing up and down quarks, G a µν is the gluon field strength tensor, µναβ ( 0123 = +1) is the completely antisymmetric tensor, D µ the gaugecovariant derivative, M the real 2 × 2 quark-mass matrix, andθ the coupling of the CP -odd interaction. In eq. (1) the complex phase of the quark-mass matrix has been absorbed in the physical parameterθ = θ + arg det(M ). For the application of chiral perturbation theory (χPT) it is useful to perform an anomalous axial U (1) transformation to replace the CP -odd gluonic term in favor of a complex mass term [19,26]. Under the assumption thatθ 1 the QCD Lagrangian can then be written as where we have defined the average quark massm = (m u + m d )/2, the quark-mass difference ε = (m u − m d )/(m u + m d ), and the reduced quark mass m * = m u m d /(m u + m d ) =m(1 − ε 2 )/2. The QCDθ term induces an EDM in hadrons and nuclei. The first EDM search occured with the neutron in 1957 [27,28]. Till this day, no signal has been found for the EDM, despite measurement sensivities having improved by six orders of magnitude. The current bound d n < 3.0 × 10 −26 e cm [29,30] sets strong limits on the size ofθ and sources of CP violation from physics beyond the SM [31].
In order to set a bound on theθ term, it is necessary to calculate the dependence of the neutron EDM onθ [19]. One way to do this is by using χPT. In the first step one derives interactions between the low-energy degrees of freedom, pion and nucleons (and heavier hadrons), that violate CP and transform the same way under chiral symmetry as the complex mass term in eq. (2). In the second step, one combines the chiral CP -odd interactions with the standard CP -even chiral Lagrangian to calculate the nucleon EDM. This calculation has been done up to next-to-leading order (NLO) in both SU (2) and SU (3) χPT [20,22] and gives in the two-flavored theory for the neutron (d n ) and proton (d p ) EDM: in terms of g A 1.27 the strong pion-nucleon coupling constant, F π 92.4 MeV the pion decay constant, m π and m N the pion and nucleon mass respectively, e > 0 the proton charge, and three low-energy constants (LECs) of CP -odd chiral interactionsḡ 0 andd p/n . The first term in brackets in eq. (3) arises from the leading-order (LO) one-loop diagram involving the CP -odd vertex L πN (θ) =ḡ 0N π · τ N (4) in terms of the nucleon doublet N = (p n) T and the pion triplet π. The LO loop is divergent and the divergence and associated scale dependence have been absorbed into the counter terms d p/n which signify contributions to the nucleon EDMs from short-range dynamics and appear at the same order as the LO loop diagrams. The second term in brackets in eq. (3) arises from finite next-to-leading-order (NLO) diagrams. Because theθ term breaks chiral symmetry as a complex quark mass, the LECḡ 0 can be related to known CP -even LECs using chiral symmetry arguments [19,32,33] where (m n − m p ) strong is the quark-mass induced part of the proton-neutron mass splitting for which we used the recent lattice results [34,35]. Inserting eq. (63) in eq. (3) we obtain d n (θ) =d n − 2.1(3) × 10 −3θ e fm , d p (θ) =d p + 2.5(3) × 10 −3θ e fm .
Under the assumption that the termsd p/n do not cancel against the calculable loop contributions, a comparison with the experimental bound gives the strong constraintθ ≤ 10 −10 . Clearly, a more reliable constraint onθ requires a direct nonperturbative calculation of the full nucleon EDMs. This is the main goal of this work. In the isoscalar combination d n + d p the loop contribution cancels out to a large extent. For observables sensitive to this combination, such as the deuteron EDM whose measurement is the goal of the JEDI collaboration [36], a first-principle calculation of the total nucleon EDM is important. EDMs of light nuclei have been calculated as a function ofθ in ref. [25]. Nuclear EDMs get contributions from the single-nucleon EDMs and from the CP-violating nucleonnucleon potential which is dominated by one-pion-exchange terms. The latter depend mainly onḡ 0 and are therefore relatively well under control. The dominant remaining uncertainty is the size of the nucleon EDMs as a function ofθ. With nonperturbative calculations of nucleon EDMs induced by theθ term, we immediately obtain predictions for EDMs of light nuclei. With future improvements of nuclear theory even EDMs of diamagnetic atoms such as 199 Hg and 225 Ra could be directly given as a function ofθ.

Lattice QCD action and numerical details
We discretize the QCD action on an hypercubic lattice with spacing a and volume L 3 × T . The fermionic part of our QCD lattice action is the non-perturbatively O(a)-improved Wilson action with N f = 2 + 1 dynamical quarks. The gauge part is the Iwasaki gauge action. For our calculation we have always used valence quarks with the same lattice action and the same bare parameters as the sea quark action, that is to say our framework is fully unitary.
We performed calculations using the publicly available PACS-CS gauge fields available through the ILDG [37]. We used 6 different ensembles that allow us to study discretization effects, finitesize effects and pion mass dependence. We studied the pion-mass dependence with 3 ensembles  at 3 different bare quark masses, at L/a = 32 and T /a = 64 and a lattice spacing a = 0.0907 (13) fm. The lattice spacing and the physical point were determined with the experimental input of m π , m K , and m Ω . More details on these ensembles are available in ref. [38] and are summarized in the first three M rows of tabs. 1, 2.
To study discretization effects we used 3 ensembles with 3 different lattice spacings but with the same volume, L 1.8 fm. The ratios of masses in the pseudoscalar and vector channels differ, between the 3 ensembles, at most by 1% in the light sector and at most of 3% in the strange quark sector. These very small mismatches are irrelevant for all purposes for our scaling violation study. The lattice spacings and quark masses in these ensembles are also determined using m π , m K , and m φ . Details for these ensembles can be found in ref. [39] and summarized in the last three A rows of tabs. 1, 2. Among the 6 ensembles described above there are 2 ensembles, M 1 and A 2 , with the same bare parameters, β = 1.9, κ l = 0.13700, κ s = 0.1364 and different lattice volumes with L/a = 20, L/a = 32 and T = 2L. These 2 ensembles allow us to investigate finite-size effects.
To improve the overlap with the ground state of the relevant matrix elements in the two-and three-point functions, we applied a Gaussian gauge-invariant smearing [40,41] at the source and at the sink of our quark propagators. Using the notation of refs. [40,41] we use 64 iterations of the smearing algorithm with a smearing fraction of α = 0.39 using the definition in ref. [40]. These parameters corresponds to a spatial root-mean-square radius for the nucleon interpolating operator of around 0.4 fm. The quality of our projection into the ground state can be evaluated from figs. 5a, 5b.
For the vector form factors studied here, we used the renormalization factor Z V determined using vector Ward identities in ref. [42] and summarized in tabs. 1, 2. The strategy we use in this paper is a perturbative expansion in powers ofθ (the expansion is fully performed in Euclidean space). This is justified by the small value ofθ estimated from experimental constraints. With this strategy every correlation function O θ , evaluated in aθ vacuum, is determined from a small-θ expansion where O is some multi-local operator,θ is the coefficient for the CP-violating term, and Q is the topological charge. The expectation values on the r.h.s of eq. (7) are computed on a standard QCD background. This allows us to use lattice QCD gauge configurations without generating new gauges for this specific calculation. We define the correlation functions used to determine the nucleon EDM induced by theθ term in the next section. The topological charge CP-violating operator that enters the correlation functions must in principle be normalized. We use the gradient flow [16] to define the topological charge which in this way has a finite continuum limit and does not need any additional normalization [16,17,18]. The reason is that the flowed fields are free from ultraviolet divergences [16,17] for all positive flow times, t f > 0. Additionally it can be shown that the topological charge defined with the gradient flow is flow-time independent [18] for all positive flow times, t f > 0, in the continuum limit. Details on how we numerically perform the flowing of the fields have been described in ref. [11].

Topological charge and the gradient flow
We define the topological charge at finite lattice spacing as where the topological charge density reads and G a µν (x, t f ) is a lattice discretization of the continuum field tensor defined with flowed gauge fields. As a lattice definition for the field tensor, we use the discretization suggested in ref. [43]. This definition suffers from small discretization effects and, in fact, the corresponding topological susceptibility is flow-time independent starting from a flow-time radius, 8t f , of about 1 fm for all lattice spacings we have investigated. This can be seen in fig. 1, where we show the flow-time dependence of the topological susceptibility computed for all M-(left) and A-ensembles (right). As expected, the region where the susceptibility is independent of the flow time extends towards smaller flowtime values for smaller lattice spacings. We used the topological charge to perform various checks on the quality of the ensembles. An important check for EDM calculations is to make sure that the ensembles sample the field space in such a way that no spurious CP-violation is induced. In other words we must check the expectation value Q(t f ) = 0 within statistical errors. In fig. 2 we show Q(t f ) evaluated on all our ensembles for various pion masses and lattice spacings. To properly estimate the statistical uncertainties we evaluate the autocorrelation function and the corresponding integrated autocorrelation time, τ int , as defined in ref. [44]. For all our ensembles the average topological charge vanishes within statistical errors. In fig. 3 we show the flow-time dependence of the integrated autocorrelation time for the topological charge. As expected the gradient flow, by smoothing out some of the short-distance fluctuations, allows a better determination of τ int that reaches a plateau for 8t f 0.2 fm for all our ensembles [45,46]. The integrated autocorrelation time τ int we obtain falls within the range 7 < τ int < 35 for the M 1 and M 2 ensembles and slightly smaller, 3 < τ int < 10, for our M 3 ensemble. We attribute this behavior with the rather short Markov Chain for the M 3 ensemble which most likely does not allow a more accurate determination of its τ int . We also observe from fig. 3 that τ int increases as we decrease the lattice spacing. This is an expected result [47,45] since the tunnelling between different topological sectors becomes increasingly difficult with decreasing lattice spacing, meaning that the sampling of different sectors, which would decrease τ int , is lessened.
For completeness in fig. 4 we show the difference of error determination if we were to use a standard resampling technique, such as bootstrap, instead of the error determination using the autocorrelation function. In this case, Q(t f ) = 0 within uncertainties. This demonstrates how a robust uncertainty determination for the topological charge requires both the estimate of the autocorrelation function and its corresponding integrated autocorrelation time.

Two-point correlation functions and the nucleon mixing angle
In this section we analyze the nucleon two-point correlation function, to extract the effective mass, as well as the nucleon mixing angle. The standard two-point correlation function with sink momentum p has the form where Π is some spin projector, and N is an interpolating field with the quantum numbers of a nucleon, inserted with a source-sink time separation of t. The spectral decomposition for this equation in the limit where T t 0, keeping implicit a sum over the polarizations, is where the lowest energy state β 0 for which Ω| N |β 0 = 0 and β 0 | N |Ω = 0 arises from the approximation 1 T t 0. The effective mass, which is shown in fig. 5, is given by the simple log ratio where Π + = (I + γ 4 )/2 is the positive parity projector, β + 0 is the lowest energy positive parity nucleon state, and again, T t 0. In fig. 5a, we compare our effective mass determinations for the M-ensembles to those computed in ref. [38] and find agreement within statistical errors. As shown in fig. 5b, we observe lattice-spacing dependence of the order of 10% between the finest and coarsest lattices.
The nucleon mixing angle [1], α N , can be extracted by utilizing the two-point correlator from eq. (11), and theθ-modified two-point correlator The mixing angle is defined using the small-θ expansion as in the region where t 0 and t f 0. In figs. 6a and 8a, we show the dependence of the nuclear mixing angle on the source-sink separation t (in fm) for the M-and A-ensembles, respectively. For the M-ensembles in figs. 6a, there is little to no excited-state contamination effects for t > 0.6 fm. The results suggest a non-trivial chiral behavior for α N . We will discuss this in detail in sec. 6 when we describe our EDM determination. For the lattice-spacing ensembles in fig. 8a, we require a minimum sourcesink separation of t ≈ a {5, : Left: m π dependence of the effective mass (in GeV) defined in eq. (13) plotted against source-sink separation time t. Bands correspond to values quoted in ref. [38]. Right: lattice-spacing dependence of the effective mass (in GeV) plotted against source-sink separation time t.  In figs. 6b we show the integrated autocorrelation time of α N for the M-ensemble results shown in fig. 6a. For the M 1 and M 2 ensembles a factor of 2 − 4 increase in autocorrelation as the source-sink separation approaches 0. Fortunately, a minimum source-sink separation of t ≈ 1 fm greatly decreases the autocorrelation correction that we apply in the determination of the nucleon mixing angle α N . Most importantly, in comparison to Q from fig. 3 (i.e. not in the presence of a nucleon), the autocorrelation effect is dramatically decreased by a factor of at least 4. We attribute this effect to the presence of a fermionic part, N N , in the correlation function. Numerical evidence suggests that the observables considered in this work containing fermion lines, such as α N and the EDM are less coupled to the slow modes contributing to the spectral decomposition of the autocorrelation function [45]. As this effect will be greater when analyzing the EDM (from three-point correlation functions), we resort to our standard bootstrap error propagation technique for the final EDM computation. We checked explicitly that error estimates from a bootstrap and an autocorrelation analysis give consistent results.

Improving the Nucleon Mixing Angle
In this section we describe a method previously explored in [48], that aims to reduce the statistical uncertainty of the determination of the nucleon mixing angle α N . The strategy can be described as an attempt to understand the space-time region where the overlap between the topological charge density and the fermionic part of the correlation function is maximal. To perform this investigation we define a spatially-summed topological charge density We then numerically study the dependence on τ Q of α N and corroborate our numerical findings with a spectral decomposition of the relevant correlators. The ratio α N and the modified two-point correlator G we therefore focus on the latter. Setting p = 0 and omitting it in our expressions, we define where the correlator in eq. (14) can be obtained by summing τ Q from 0 to the time extent of the lattice T G (Q) To focus on the region where the signal resides, we sum the spatially-summed topological charge density, Q(τ Q , t f ), symmetrically starting from the source location. That is we sum τ Q starting from 0 (and T ) up to a value t s (and T − t s ). The goal is to find a summation window t s small enough such that we capture all the signal and avoid the summation of unnecessary "noise". We define the partial summed correlator from which, using the periodicity of our lattice, the original correlator in eq. (14) is obtained as Although there are other choices for the starting point of our summation in τ Q , we only consider starting from τ Q = 0. In app. A we derive a spectral decomposition for the correlator in eq. (17).
We argue that in the limit t s t 0, the partially-summed correlator G (Q) 2 (t, Π, t f , t s ) is independent of t s and t, up to exponentially suppressed corrections. These corrections seem to be rather small, and in fact our numerical experiments indicate that we can safely stop the summation over τ Q at t s t. In this way we avoid to sum in the region between t and T /2 where numerically the correlators seem to vanish up to statistical fluctuations.
We first fix the source-sink separation t to a large enough value such that effects from excited states are suppressed. We then study the dependence of α N on the summation window t s . In fig. 10 we show the t s dependence of α N , for the M-ensembles (left), and A-ensembles (middle) and two different physical volumes corresponding to M 1 and A 2 ensembles (right). In all ensembles we observe that α N reaches a plateau when t s t, consistent with the expectation that contributions for t s > t are exponentially suppressed and below our statistical accuracy. We do observe a very small drift of α N for larger values of t s for the ensembles M 1 , and a smaller drift for the ensemble M 3 , for small values of t. We attribute this to statistical fluctuations that could arise from small local parity-violating effects induced by non-vanishing matrix elements of Q(τ Q , t f ) between two states of the same parity, β| Q(τ Q , t f ) |β = 0. These local fluctuations are averaged out when the charge density is summed over the whole space-time volume as shown in fig. 2. Nevertheless, all values of α N determined with the improved method are statistically compatible with the results obtained with the standard analysis.
To compare the improved extraction of α N to the standard determination described in sec. 5, we show in figs. 11, 12 the standard and improved determination of α N as a function of the  Euclidean source-sink separation t. The values of t s considered are summarized in tab. 4. We observe a signal-to-noise improvement in all our ensembles, up to a factor 2, with the most significant observed in the ensembles M 1 , M 2 and A 3 . We observe the largest discrepancy between the improved and unimproved methods, of the order of 2.3 σ, in the M 1 and M 3 ensembles. We attribute this discrepancy to standard statistical fluctuations of the gauge fields. A summary of the fit ranges and results for the improved nucleon mixing angle is given in tab. 4, where for companions we added the values of α N determined in the standard way.

Electric Dipole Moment Results
The neutron (n) and proton (p) EDMs, d p/n , can be extracted from the CP-odd electric dipole form factor 2 F p/n which requires a lattice QCD computation of F 3 (Q 2 ). The variable Q 2 in this case refers to the momentum transfer and should not be confused with the topological charge. The smallθ expansion provides us a way of accessing F 3 from three-point correlation functions without the need for generating new gauge configurations at finiteθ and without relying on a problematic analytical continuation to imaginaryθ. To access F 3 (Q 2 ), we calculate the following three-point   (11) correlation functions with and without the insertion of the topological charge, respectively, where the electromagnetic current in terms of the quark currents is given by and N denotes standard proton or neutron interpolating fields.
Once the three-point correlation functions are computed, we remove the leading Euclidean time dependence and nucleon-to-vacuum amplitude contributions via the ratios: where we have implicitly defined ratios for the proton and the neutron. We define the square-root factor as The spectral decomposition of the ratio function R in eq. (24), in the limit T t 0, reads where the vector form factor contains all terms allowed by the symmetries of the theory For completeness the expression of A(E p , E p ) reads The data that we computed coming from the fixed-sink method is such that p = 0 (which implies q = −p), simplifying eq. (26) to The analogous modified ratio function with the insertion of the topological charge R (Q) in eq. (24) has, retaining only the ground state contribution, the following spectral decomposition at leading order inθ Due to subtleties between lattice quantities and physical, the form factor decomposition in presence of a CP-violating operator insertion, is written in terms of modified form factors, F 2 (Q 2 ) and F 3 (Q 2 ), related to the physical form factors by [6] The rotated form factor F 3 (Q 2 ) corresponds to the actual electric dipole form factor as measured in experiments. From now on, we will focus on this quantity. The ratio functions R and R (Q) become constant, as long as the large-time approximation T t τ 0 is satisfied to ensure ground-state dominance. As the fixed-sink method is employed to compute the three-point correlation functions, a region in which this large time approximation is satisfied for τ can be found and we denote the results of the fits as   The technique for fitting these ratio functions over τ is described in app. B. With this construction, a system of equations can be solved for form factors F i (Q 2 ) , i = 1, 2, 3 of the form: where the collective index A denotes any combination of the indices A = {j, k, l}. In other words, we run over all possible combinations of projectors Π, all current momentum q within a given Q 2 , and operator gamma matrix γ µ . The index A of the matrix A Ai (Q 2 ) corresponds to the coefficients for each form factor F i for the corresponding ratio function R or R (Q) . These coefficients are found by analyzing the spectral decomposition of R or R (Q) , which needs to be done for every evaluated index A. Using eq. (21), we extrapolate to F p/n 3 (Q 2 → 0)/(2M N ) = d p/n . We use a linear plus constant fit function, giving the extrapolated value d p/n at Q 2 → 0 (as well as slope in Q 2 providing S p/n ).
The final extraction of the neutron (left) and proton (right) CP-odd form factor F 3 (Q 2 ) 2M N is shown for the M-ensembles in fig. 13 and for the A-ensembles in fig. 14. Fig. 13 shows that all M-ensembles are statistically consistent evaluated, and with zero. Fig. 14 shows that there are no major discretization effects, as all the extrapolated Q 2 → 0 results are consistent.
The following figs. 15, 17, 16 are all displayed to understand the systematic effects resulting from varying the flow time t f , and different methods of determining the nucleon mixing angle α N used in the form factor decomposition (A Aj (Q 2 ) in eq. (35)). In fig. 15, for example, we show how the form factors F 3 , determined at different flow-time radii 8t f = 0.60, 0.65, 0.70 fm (green, red and blue), are statistically consistent for all three M-ensembles (left to right). From both fig. 16, where the improved method (see sec. 5.1) of determining the nucleon mixing angle α N (in red) is compared to the standard method for α N (in blue), and fig. 17, where we vary the fit range for extracting α N , it is clear that a more precise determination of α N has a negligible impact on improving the precision of the results for the CP-odd form factor F 3 . A summary of the Q 2 → 0 extrapolations for different ensembles is given in tabs. 5, 6.

Improving the Modified Three-Point Correlation Function
In this section, we utilize a similar improvement technique used for α N , but now applied to the modified three-point correlation function G  dependence of the spatially integrated topological charge density where τ Q signifies the temporal location of the topological charge Q defined in eq. (16). The spectral decomposition for this correlator has the form: where α, β, γ and δ are labels for the states propagating, and the 0 subscript indicates the lowest energy state propagating with the appropriate quantum numbers. We stress that t f = 0 implies the absence of any contact terms. From fig. 18, a clear signal is observed at τ Q = 0 on all ensembles. This motivates summing τ Q symmetrically around τ Q = 0 to obtain the summed three-point correlator: The resulting fit function to the sum range t s , for the range t s > t, is: Where γ ± and β ± represent the positive and negative (±) parity nucleon states. A 0 , A γ ± , A β ± ,γ ∓ , and A β ± are combinations of nucleon matrix elements, E γ ± is the energy of the propagating state γ ± , and E 0 ± is the lowest energy of the positive and negative parity nucleon states 0 ± . We construct the improved ratio function R (Q) in the same way as in eq. (24), but with the 3 . The value of the correlation function G (Q) 3 , used to extract the CP-odd form factor, is obtained in the limit t s → T . If the summation over τ Q is performed up to a value t s < T the neglected terms will be exponentially small as one can deduce from eq. (39). Our numerical results seem to indicate that indeed the neglected contributions for intermediate values of t s are well below the statistical accuracy of our calculation.
In fig. 19, the results for the symmetrically summed topological charge ratio function R (Q) are shown as a function of the sum range t s . In all cases, a plateau can be observed at t s = τ . This indicates that all the exponential terms in eq. (39) are suppressed for t s > τ . Coupled with the large statistical noise inherent in the data, we fit the result with a constant value once the plateau has formed. These fit ranges are displayed in tables 7 and 8, and are used for the form factor analysis in sec. 6.2.
Finally, fig. 20 displays a standard modified ratio function R (Q) 3 plot over current insertion time τ , where the improved ratio function (blue) is compared with the standard method (red). The improved ratio function uses the "min" time from tables 7 and 8.
In fig. 21, a comparison between the improved ratio functions (blue) and the standard integrated topological charge (red) used in the extraction of the neutron CP-violating form factor 2M N is shown. In all cases, a two-to-three times increase in the signal-to-noise is observed and all results are statistically consistent 3 .

Continuum Extrapolated Results with Improved Ratio Functions
Armed with the improved results for the nucleon EDMs, the next step entails the extrapolation to the physical pion mass and the continuum limit. From χPT we learn that the leading      (13) dependence of the nucleon EDMs on the pion mass is given by [49] d p/n (m π ) = C 1 m 2 π + C 2 m 2 π log( where C 1 and C 2 are fit constants. To account for the finite lattice spacing, we include an additional fit parameter, C 3 , The additional term ensures that the EDM only vanishes in the chiral limit after taking the continuum limit. We have performed a global fit with eq. (41) taking into account our 6 data points from ensembles A 1 -A 3 and M 1 -M 3 . In the four plots in figs. 22, 23, we show the EDM results for the proton and neutron separately as function of the pion mass and lattice spacing. Specifically, in Fig. 22 we show the extraction of the neutron (left) and proton (right) EDM plotted against their m 2 π values (in MeV). The blue band shows the extrapolation using the fit function in eq. (41), evaluated at d p/n (a = 0, m π ). This function evaluated at the physical pion mass is what we are interested in. In red we show the same extrapolation, where the fit is evaluated instead at d p/n (a = 0.09 fm, m π ), to study the role of discretization errors. In particular, we observe an uncertainty of the EDMs at the physical pion mass that is roughly twice larger at a = 0.09 fm. It is perhaps surprising that the uncertainty at the physical point reduces in the continuum limit. But the reason is clear. By fitting the nucleon EDMs to the fit function in eq. (41), the uncertainty on the fit parameters C 1 and C 2 is increased by the presence of the C 3 term. Now that the a 2 dependence is taken into account, we can perform Table 9: Neutron and proton EDM fit parameters C 1 , C 2 , C 3 extracted from the combine fits to all 6 ensembles using eq. (41), as well as the resulting χ 2 P DF . We also estimateḡ 0 using C 2 and eq. (3). an interpolation between the EDM in the chiral limit and the pion masses of our ensembles. In the continuum limit, a = 0, the resulting nucleon EDM at the physical pion mass has now less uncertainty because d n,p (a = 0, m π = 0) while d n,p (a > 0, m π = 0) = 0 and unconstrained.
and we include the determination for the fit parameters C 1 , C 2 , C 3 of eq. (41), as well as the chi-squared per degree of freedom parameter, χ 2 P DF , in tab. 9. The error on χ 2 P DF is determined from the bootstrap samples distribution. Since the correlators for the proton and the neutron EDM are different, it is possible to obtain different relative uncertainties in the two cases. It is not clear to us though, why we observe a relative larger uncertainty for the proton than for the neutron.
In fig. 23 we show the dependence of our EDM results on the lattice spacing a for the neutron (left) and proton (right) EDM. Overlaid on top, we have the evaluation of the fit function eq. (41) at two different values of m π : m π = 700 MeV (purple band) and m π = m phys π (green band). The ensembles analyized in this work do not allow us to study mass-dependent discretization effects, but we can still observe the impact of the chiral interpolation on the continuum limit. The continuum extrapolation has less uncertainty, thanks to the constraint that the EDM vanishes in the chiral limit. Adding more ensembles to study mass dependence cutoff effects is certainly desirable, but it does not change the main conclusion of this analysis.
We can extract a value of the CP-odd pion-nucleon LEC,ḡ 0 , which plays an important role in the EDMs of nuclei and diamagnetic atoms, by identifying our result for fit parameter C 2 with the coefficient of the log term in eq. (3) for the neutron EDM. This gives the relation leading to the extractionḡ at the physical pion mass. This result is in good agreement with the chiral perturbation theory prediction in eq. (5) and confirms the applicability of the fit function in eq. (41). A consistent result, with larger uncertainties is obtained for the proton EDM (see tab. 9). , (green) and in the chiral limit m π = 700 MeV (purple).

Schiff Moment of the Proton and Neutron
Apart from the EDMs of the neutron and the proton, the nucleon electric dipole form factor (EDFF) contains additional information. The EDFF can be decomposed as where d p/n denotes the proton or neutron EDM, S p/n denotes the proton or neutron Schiff moments defined by S p/n = (2M N ) −1 (dF p/n 3 /dQ 2 )| Q 2 =0 , and H p/n are functions that capture the remaining Q 2 dependence. Chiral perturbation theory allows for a calculation of the Schiff Moments and the H p/n functions from the analogous isovector and isoscalar quantities [49]. At leading order in the chiral expansion the nucleon EDMs are given in eq. (3). The leading-order Schiff moments are isovector and given by where we have used eq. (5). NLO corrections have been calculated in ref. [22] and reduce the leading-order result by roughly 50% and provide a tiny contribution to S n +S p = O(10 −5θ e fm 3 ). χPT thus predicts that the neutron and proton Schiff moments are equal in magnitude but with opposite sign. The leading-order H p/n are also isovector and given by where h (0) In the limit Q 2 m 2 π the nucleon EDFFs become such that the Schiff moments provide the dominant Q 2 dependence of the EDFFs. The nucleon EDMs and the LECḡ 0 are induced by theθ term and scale as d p/n ∼ḡ 0 ∼m * θ ∼ m 2 πθ . As such, the Schiff moments scale as S p/n ∼ḡ 0 /m 2 π which is pion mass independent. This statement is potentially confusing as we infer from eq. (2) that theθ term decouples in the chiral limit and the whole nucleon EDFF should vanish. eq. (50), however, requires Q 2 m 2 π . In the opposite limit, we obtain and the EDFFs vanish in the chiral limit as expected.
The goal is to extract S p/n from our lattice data as this allows for a direct comparison to the χPT prediction in eq. (46) and the extraction in the previous section based on the pion mass dependence of the nucleon EDMs. To extract S p/n , we first extrapolate our results to small Q 2 by fitting the EDFF to the function The effects of the h (0) 1 function turns out to have minimal impact on the extraction of d p/n (m 2 π , a 2 ) and S p/n (m 2 π , a 2 ), and we obtain similar results if we use the fit function This shows that our results are not precise enough to isolate the more subtle Q 2 behavior. Once we have obtained S p/n (m 2 π , a 2 ) we can extrapolate to the continuum limit and the physical pion mass. LO χPT predicts no dependence on the pion mass, and, having an O(a) improved lattice action, we add a quadratic dependence on the lattice spacing a with C 4 and C 5 fit constants. The results for the Schiff moments along with the continuum extrapolation are shown in figs. 24, 25. In fig. 25 we show the fit results with the a 2 dependence. We observe minimal discretization effects over the range a = {0 → 0.12} fm. In fig. 24 we show the fit results as a function of the pion mass m π . At a = 0 we perform a constant fit in the pion mass, to obtain the proton and neutron Schiff moments at the physical point. We do not extrapolate to the chiral limit because the χPT prediction that S p/n are pion-mass independent will break down at some point as inferred from eq. (51). We obtain for the Schiff moments at the physical point as well as the fit parameters C 4 = S p/n and C 5 from performing this fit in tab. 10. The uncertainties are significant and the magnitudes are somewhat below the LO χPT predictions in eq. (46), but in better agreement once χPT NLO corrections are included. There is some evidence for a dominantly isovector Schiff moment as predicted from χPT, but the uncertainties are too large to make strong statements. We perform a sanity check of our result by comparing the ChPT predictions for the fit coefficient C 2 and C 4 . From eqs. (43) and (46), we infer the LO ChPT prediction Our fit values for this ratio are given in tab. 10, and agree with this prediction within (large) statistical errors.  Figure 24: Determination of the Schiff moment S p/n for the neutron (left) and proton (right) for all 6 of our ensembles, plotted against their respective m π values. The bands are the fits to all the ensembles using eq. (55), evaluated in the continuum limit a = 0 (blue) and the lattice spacing corresponding to the M-ensembles (red). Table 10: Neutron and proton Schiff fit parameters C 4 , C 5 extracted from the combine fits to all 6 ensembles using eq. (55), as well as the resulting χ 2 P DF . We additionally include the ratio C2 C4 , where C 2 is the second fit parameter result from tab. 9.

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

Comparison with other works
In this section, we compare the results obtained for the neutron EDM d n with few lattice QCD results from the literature. As noted in [6], it is sometimes problematic to compare different EDM calculations, as most results preceding this paper do not consider the rotation of the CPodd form factor F 3 with F 2 and α N computed on the lattice. 4 The EDM d n (rotated to F 3 ) is shown in Fig. 26. Good agreement is seen from our results to the others [1,3] at m π ≈ 475 MeV, but we see a slight tension between the results of [4,3] and our results at m π ≈ 350 MeV. It must be stressed that the rotation requires knowledge of the phase α N and the unrotated form factors F 2 and F 3 which are not always easy to extract. To rotate the "C. Alexandrou et al, 2016 " [4] results, an estimation of F 2 was determined from [50]. To rotate the "F.-K. Guo et al, 2015 " [3] results, F 2 was determined from [51] (atθ = 0) and α N and F 3 estimated via a linear+cubic fit inθ performed by [6].
In particular, the lattice results for F 3 not obtained in this work do not take into account correlations between F 2 , F 3 and α N . As such, Fig. 26 is mainly shown for illustrative purposes and the error estimates for results not obtained in this work should be taken with a grain of salt.

Impact on EDMs of light nuclei
Armed with a non-perturbative determination of the nucleon EDMs as a function ofθ we can revisit EDMs of systems with more than a single nucleon. EDM experiments so far have mainly focused on neutral systems, but EDMs of charged particles can be probed if the particles are trapped in electromagnetic storage rings [52]. This technique has lead to a direct limit on the EDM of the muon [53], and to plans to pursue EDM measurements of protons and light nuclei in dedicated storage rings. Such measurements are still far away but impressive progress has been reported in refs. [54,55]. EDMs of light nuclei have been calculated using chiral EFT [24,25] d2 d3 in terms of the EDMs of nucleons and the CPV pion-nucleon coupling constantsḡ 0 andḡ 1 associated to the interactions  The results of d n from this paper, improved in light blue and not-improved in dark blue, compared to other lattice QCD results [4,3,1]. The light blue bands correspond to a chiral extrapolation, using the improved data (light blue), after having performed the continuum limit as described in sec. 6.2. To perform the rotation of the CP odd form factor F 3 of other lattice calculations see the main text. We underline that our error determination for other works is purely illustrative since, not having at our disposal the raw data, we do not take into account correlations in the data.
Values ofḡ 0 andḡ 1 can be obtained from chiral-symmetry arguments [56,33] that link these LECs to the hadron spectrumḡ and the smallness ofḡ 1 /ḡ 0 is due to approximate isospin symmetry.
In absence of direct lattice calculations of d n and d p we could only predict values for the combinations But with our lattice determination of d n (θ) and d p (θ), we can now estimate the EDMs of light ions directly in terms ofθ Due the dependence on the isoscalar nucleon EDM, d n + d p , the deuteron EDM is still very uncertain. The tri-nucleon EDMs however are predicted more than three standard deviations from zero and with a fixed sign. The total uncertainty arises in roughly equal amounts from uncertainties in eq. (63) and in the determination of the nucleon EDMs in eq. (42). If nonzero EDMs are measured these relations can be used to differentiate between the SM theta term and BSM sources of CP violation. They can also provide indirect evidence for the existence of a Peccei-Quinn mechanism, by finding EDM patterns in disagreement with eq. (65) [57].

Conclusion
In this paper we computed the proton and neutron EDM from dynamical lattice QCD using various pion masses at different lattice spacings and volumes, as enumerated in tabs. 1, 2. We found our results have rather small (within our statistical uncertainties) discretization effects, which greatly simplified our continuum limit extrapolations. We found satisfactory agreement with existing results, as discussed in section 7. With our measurements at multiple pion masses, we performed a chiral interpolation to obtain, at the physical pion mass and in the continuum limit d n = −0.00152(71)θ e fm and d p = 0.0011(10)θ e fm. The nonzero result for the neutron EDM confirms the existence of the strong CP problem at two standard deviations and limits θ < 1.98 × 10 −10 at 90% C.L. The dependence of the nucleon EDMs on the pion mass allowed us to extract of the CP-odd pion-nucleon coupling,ḡ 0 . The resulting value is in good agreement with chiral perturbation theory. Important to our analysis was the implementation of the gradient flow on the topological charge. In fact we can perform the continuum limit at fixed flow time with no need to calculate the normalization of the topological charge. As we have discussed and documented in section 5 and ensuing subsections, the gradient flow also allows a more robust determination of the integrated autocorrelation time that must be taken into account when estimating statistical uncertainties.
Be that as it may, the extraction of a non-zero EDM is notoriously difficult due to its poor signal-to-noise ratio. To address this issue we have employed a novel technique, first documented in [12], for reducing the noise in our measured observables. Instead of summing all space-time points in the calculation of ratios relevant for the extraction of our 3-point function and α N term, this method focuses on the space-time region where the signal is strongest. We argued that the neglected space-time region gives exponentially suppressed contributions to the correlation functions and this expectation has been confirmed by our numerical results. On some ensembles, this method enabled us to increase the signal to noise by a factor of ≈ 2. This method was described in detail in subsection 6.1 and Appendix A.
We have also analysed the Q 2 dependence of our form factors and performed an extraction of the Schiff moment with dynamical fermions. Our results of S p = 0.50(59) × 10 −4θ e fm 3 for the proton and S n = −0.10(43) × 10 −4θ e fm 3 for the neutron that are in reasonable agreement with chiral perturbation theory predictions. Our estimates for this value can be improved upon with more statistics and calculations on larger lattices (and thus lower Q 2 points), which would allow for a more robust extraction.
Our calculation is a big step towards a precise determination of the nucleon EDM and Schiff moment. Improvement of these results will most definitely come from increased statistics, and more calculations at different pion masses at several lattice spacings. We comment here on the necessity to perform calculations at the physical pion mass. In the chiral limit the EDM induced by the θ term vanishes (i.e. d p/n = 0 at m π = 0). In our view, given the small value of the EDM induced by the θ term and the additional standard reduction in signal-tonoise as the pion mass is lowered for nucleon correlators, calculations of these quantities at the physical pion mass have possibly less to gain than those at higher pion masses. Because of this constraint, it could be advantageous to have results at slightly heavier-than-physical pion masses and then robustly interpolate to the physical pion mass using χPT. The subsequent errors of the interpolation are stable and easily quantified precisely because one is doing an interpolation and not an extrapolation. We remark though that, for a chirally-breaking lattice action, such as the non-perturbatively O(a) improved Wilson-clover fermion action we have adopted, the nucleon EDMs vanish in the chiral limit only after performing the continuum limit. This emphasize the importance of the continuum limit when using a chirally-breaking action. In this respect the gradient flow allows us to perform a safe study of discretization effects.
To summarize, the ideal scenario of a direct determination at the physical point with statistical uncertainties under control, can be circumvented by simply investing more time in lattice QCD calculations at slightly heavier pion masses (where the signal-to-noise is not as prohibitive). It goes without saying that calculations of the EDM at heavier-than-physical pion masses can potentially be more cost effective than the physical-pion-mass calculations only if one has a robust description of the lattice data with χPT.
cial support from the Deutsche Forschungsgemeinschaft (Sino-German CRC 110). The authors gratefully acknowledge the computing time granted through JARA-HPC on the supercomputer JURECA [58] at Forschungszentrum Jülich.

A Alpha Improvement Derivation
Staring with the general three-point correlation function ∆ (O) 2 (p , t; q, τ ; Π) = a 6 x,y e −ip ·(x−y) e iq·y Tr Π N (x, t)O(y, τ )N (0, 0) , we handle the time ordering in the next two sections, by performing the spectral decomposition for t > τ and τ < t. This is the general expression for an arbitrary operator O and spin projector Π. For the computation of the nucleon mixing angle α N , we have O = Q.

A.1 Case t > τ
Starting with the specific time ordering t > τ in eq. (66), we perform the standard spectral decomposition to produce the correlation function where the sum over states α, β, γ have been reduced to states that only contain momenta p γ = q, p β = p − q and p α = p. The two approximation one can apply to this equation are T t and t 0, which are related to the source-sink separation of the two-point correlation function where α 0 is the lowest lying energy state that gives a non-zero contribution to ∆ The two approximations T t and t 0 are again applied where this time, γ 0 is the lowest lying state that gives a non-zero contribution to ∆ Tr{Π α| O |β β| N |γ 0 γ 0 | N |α } = 0. (72)

A.3 Total Form
Over the total range τ ∈ [0, T ], the expression is

A.4 Symmetric Partially Summed Current t s > t
This region is where the fit will take place, so we omit the derivation for t s < t.
noting the second sum is shifted to τ ∈ [t + a, t s ] using the lattice spacing increment a. One thing to note here, is the terms τ = 0 and τ = t are contact terms, which need to be handled properly (operator product expansion, or gradient flow). Next we group τ terms in preparation for the t s sum The sums can be computed by using: substituting in this expression, and including the terms where α = β and γ = β separately As the final 2 terms are only exponentially dependent on t s , we write these terms as exponentials of single energy indices

A.5 Explicit form for O = Q
As Q is a parity violating operator, the nucleon states that propagate before and after this operator must be opposite in parity. This removes the first and second terms as β| Q |β = 0. As well as this, the terms with sums over two terms either require (α, β = α + , β − ) or (α, β = α − , β + ) where the subscript ± refers to the state having positive or negative parity. The projector is selected to be Π = γ 5 Π + = γ 5 The terms e −Eα + (T +t) and e −Eα + (T −ts) in the final sum are exponentially suppressed as T T /2 ≥ t s and T t G (Q) 2 (p , t; q, t s ; γ 5 Π + ) = a β =γ From this complicated expression, the t s dependence only appears exponentially in the final term. Therefore, we can fit the two-point correlation function with f it(t s ) = A + Be −Ets . (82) Due to the statistical noise of the data and high correlation in the data with respect to t s , we elected to neglect the excited state term by fitting a constant in the region where Be −Ets A.

B Ratio Function Fit Range Selection
In this appendix, we present the technique used for extracting the CP-odd form factor F 3 (Q 2 ) from the ratio function in eq. (34). Since only constant ("one-state") fits are implemented for the ratio functions, careful consideration to excited state effects is needed. The method employed to account for fit range dependence in our error estimates, is to include multiple fit ranges that satisfy some χ 2 per degree of freedom (χ 2 P DF ) criterium. For this study, we only select fits that satisfy χ 2 P DF ∈ [0.5, 1]. Using the multiple fit range deturminations of R and R Q , we extend eq. (35) to include different fit ranges: where the extra index f (A) refers to which fit range is used, which depends on the collective index A = {j, k, l}, A ∈ [1, . . . , N A ].
Since each ratio function selected by index A has f (A) different ways to extract the quantity, the system is solved for every combination of f (A) ∀ A ∈ [1, . . . , N A ]. This results in A F (A) independent system of equations to solve, where F (A) is the number of different fits accepted (using the χ 2 P DF criterium) for index f (A) = 1, 2, . . . , F (A). Once the form factors have been solved over different fit range combinations, the result we obtain is F i,f (Q 2 ), where the (A missing) index f refers to which combined set of fit ranges were used. Since the extrapolation to Q 2 → 0 must be performed to compute the nucleon EDM, this must be performed for every f (Q 2 ) combination (analogous to f (A) above). So in addition to above, we increase the number of fits to Q 2 F (Q 2 ), where F (Q 2 ) is the number of fits computed for index f (Q 2 ) = 1, 2, . . . , F (Q 2 ).
Combining both these studies together, the resulting nucleon EDM has been computed using different fit ranges, indexed by

B.1 Computational Viability
As one may notice, the above formulation is of order O(A!), assuming a fixed number of fit ranges selected. A stochastic estimation of the fit range variation is highly recommended, which can be employed when solving the form factor eq. (83), as well as when taking the form factor Q 2 → 0 extrapolation.
At the form factor solving stage, this is employed by randomly selecting N χ different fit range that satisfy the χ 2 P DF criterium. The resulting number of systems of equations to be solved are N N A χ .
For the form factor extrapolation in Q 2 → 0, a random selection of N F results of index f (Q 2 ) in F N,f (Q 2 ) (Q 2 ). The resulting number of fits to be performed using this estimation is N

B.2 Results computed in this paper
The results computed in this paper use the fit criterium χ 2 P DF ∈ [0.5, 1], and excluded fits of length 2. The cutoff for the number of fit ranges per ratio function is N χ = 4 and the cutoff for the form factor extrapolation is N Q = 4 as well.
As for the the number of equations to solve, results at lattice (q/a) 2 = 1, 4 has 1024 equations, (q/a) 2 = 2 has 4096 equations and (q/a) 2 = 3 has 16384 equations. Multiplying these numbers by 200 bootstrap samples, will give the individual number of system of equations solved. Once this is complete, we avoid computing ∼ 70 trillion equations by performing the stochastic estimate which only requires 256 equations to solve for. Although it may seem the values for N χ and N F are insufficient in size, the results shown in fig. 27 demonstrates minimal variation when analyzing each individual