Neutron electric dipole moment using lattice QCD simulations at the physical point

We extract the neutron electric dipole moment $\vert \vec{d}_N\vert$ within the lattice QCD formalism. We analyse one ensemble of $N_f=2+1+1$ twisted mass clover-improved fermions with lattice spacing of $a \simeq 0.08 \ {\rm fm}$ and physical values of the quark masses corresponding to a pion mass $m_{\pi} \simeq 139 \ {\rm MeV}$. The neutron electric dipole moment is extracted by computing the $CP$-odd electromagnetic form factor $F_3(Q^2 \to 0)$ through small $\theta$-expansion of the action. This approach requires the calculation of the topological charge for which we employ a fermionic definition by means of spectral projectors while we also provide a comparison with the gluonic definition accompanied by the gradient flow. We show that using the topological charge from spectral projectors leads to absolute errors that are more than two times smaller than those provided when the field theoretic definition is employed. We find a value of $\vert \vec{d}_N\vert = 0.0009(24) \ \theta \ e \cdot {\rm fm}$ when using the fermionic definition, which is statistically consistent with zero.


Introduction
The discrete symmetries of parity P , charge conjugation C and time-reversal T play a significant role in the phenomenological structure of the Standard Model (SM). Although the neutron has an overall zero electric charge, a conundrum that physicists are trying to understand for decades is whether the distribution of the positive electric charge could possibly not coincide with the distribution of the negative electric charge. That is to say whether there is no invariance under CP parity, which can address the well known unsolved puzzle of the origin of the imbalance of matter and antimatter in the universe. Such an imbalance between positive and 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 current bounds are still several orders of magnitude below the SM prediction on CP violation induced by the weak interactions. A finite nEDM value would point towards physics beyond the standard model (BSM) [1] and it 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 d N provided by a number of experiments, such as those given in Refs. [2][3][4]. Until recently, the best measured upper bound was that given in Refs. [3,4] as | d N | < 2.9 × 10 −13 e · fm (90% CL) . ( 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 precession frequency. 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 Switzerland using Ramsey's method of separated oscillating magnetic fields with "ultra-cold" neutrons. The novelty of this experiment lies on using the Hg-199 co-magnetometer and an array of optically pumped cesium vapor magnetometers to cancel and correct for magnetic field changes. The result of this experiment yields an improved upper limit of | d N | < 1.8 × 10 −13 e · fm (90% CL) .
The nEDM value has also been extracted from different theoretical perspectives, such as model dependent studies [6][7][8][9][10][11][12][13][14], as well as effective field theory calculations [15][16][17][18][19][20][21][22][23], reporting values for d N in the range of | d N | ∼ θ · O 10 −2 − 10 −3 e · fm. Using either the result given in Eq. (1) or Eq. (2) one derives a bound of the order θ O 10 −10 − 10 −11 . 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 the lowest-dimensional (d = 5) effective quark-gluon CP -odd interaction [24,25]. In this approach, we include a θ-topological term in the QCD Lagrangian density: written in Euclidean time. The first two terms are CP -conserving while the θ-term is CP -violating that can give 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, is defined as By considering the electroweak (EW) sector of the SM, the action in Eq. (3) receives a contribution from the 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 and EW 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 bound on the nEDM value. This is referred to as the "strong CP problem". We assume that a near cancellation of θ and arg det M is extremely difficult to occur, thus, we take θ to be small and we perturbatively expand around it.
At low momentum transfer, the nucleon effective Lagrangian gives rise to the expression for the nEDM [1]: at leading order in θ, where m N denotes the mass of the neutron, Q 2 = − q 2 the four-momentum transfer in Euclidean space (q=p f − p i ) and F 3 (Q 2 ) is the CP -odd neutron form factor. We can, therefore, calculate the electric dipole moment by extracting the zero momentum transfer limit of the CP -odd form factor. This is the framework on which our work is based on. In practice, it is impossible to extract F 3 (Q 2 ) at Q = 0 due to the fact that the CP -violating nucleon matrix element, decomposes to Q k F 3 (Q 2 ) (k=1, 2, 3) and not just to F 3 (Q 2 ). Hence, a direct extraction of F 3 (0) is prohibited by the theory. We then need to parametrize the Q 2 -dependence of F 3 (Q 2 ) and then take the limit Q 2 = 0.
Since we expect that the θ-parameter is very small we can expand the exponential of the topological term in powers of θ. This enables us to integrate out (over space-time) the topological charge density, which gives the total topological charge Q: A consequence of this expansion is that the value of F 3 (Q 2 ) depends on the nucleon two-and threepoint functions correlated with the topological charge Q, namely with respectively. While the computation of nucleon two-and three-point functions within lattice QCD in the absence of the parity-violating term is well-known and rather precise [26][27][28], when correlating with Q, introduces large statistical fluctuations. The topological charge, if sampled adequately, has a gaussian distribution centered at zero with a width that increases with the volume of the lattice and would by itself average to zero. If there is no correlation, or mild correlation, between the correlation functions and the topological charge then no signal can be obtained for F 3 (Q 2 ). Different definitions of the topological charge may be more suitable to pick up the correlations with the nucleon two-and three-point functions [29,30].
Recently, a comparison of lattice QCD determinations of the topological susceptibility using a gluonic definition with gradient-flow, cooling and other equivalent smoothing schemes, as well as, a fermionic definition using spectral projectors [31][32][33][34] has been carried out [29,30]. Although both these gluonic and fermionic definitions lead to compatible results in the continuum limit, the gluonic definitions are much more affected by cut-off effects. For the spectral projectors one can choose a value of the spectral cutoff that would completely eliminate discretization effects in the topological susceptibility. Therefore, using spectral projectors is preferable since the discretization error on the topological susceptibility as well as the statistical error are suppressed. Given the statistical superiority of spectral projectors in calculation of quantities that depend on the topological charge, one needs to consider the fermionic definition as compared to the commonly used gluonic definition. One has to bear in mind that the computational cost of spectral projectors is orders of magnitude larger than that of the gluonic. However, for quantities calculated at a given lattice spacing without a continuum extrapolation, the fermionic definition can be the only option to avoid large cut-off effects. Since in this work we investigate the nEDM at physical pion mass we are, at the moment, restricted to use one N f = 2 + 1 + 1 ensemble [35], and therefore a continuum extrapolation is not possible. In the future, we plan to include two more ensembles with smaller lattice spacings in order to be able to take the continuum limit.
An alternative way to compute the nEDM is to use an expression from chiral perturbation theory to extrapolate the nEDM obtained at larger pion masses to the physical point. However, such an extrapolation in the pion mass carries its own uncontrolled systematic errors. Current values of F 3 (m π ) computed for m π > 135 MeV within lattice QCD, as can be seen in Fig. 10, are non-zero only within a few standard deviations or even consistent with zero within the given statistical accuracy. Furthermore, F 3 decreases as the pion mass decreases vanishing at the chiral limit. Thus, the smaller the pion mass the more difficult it is to compute F 3 (Q 2 ) with a statistical significance that would exclude a zero value. Therefore, performing a chiral extrapolation from such data can be difficult. In a recent study [36], the authors performed a chiral fit and obtained a statistically significant non-zero value for F 3 (0) at pion mass m π = 135 MeV, concluding that there is a signal for CP violation induced by the θ-term. However, large uncertainties on the few data points used in the fit, cast in doubt the reliability of the result. It is, therefore, very important to perform a first-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 3 (Q 2 ) from two and three-point correlation functions. Subsequently, in Section 3 we provide the details of the lattice QCD setup and the calculation, including the expressions of the correlation functions in terms of the topological charge and discussion on the topological properties (Section 4) of the ensemble used in this work. In Section 5 we present our results for the nucleon mixing angle and the CP -odd form factor F 3 (Q 2 ) as we take the Q 2 = 0 limit. We compare with other studies in Section 6 and finally, in Section 7 we conclude.

Method
The precise computation of the nEDM from first principles is one of the long-standing challenges and an active topic of research within lattice QCD. The lattice QCD formulation provides an ideal framework to access non-perturbatively the nEDM. The first pioneering attempt was reported nearly three decades ago [37], and was based on the introduction of an external electric field and the measurement of the shift in the associated energy. This approach of extracting the nEDM within lattice QCD, however, breaks unitary since it requires the introduction of an external electric field. Subsequently, two new approaches have been proposed. The most commonly used method involves the calculation of the CP -odd F 3 (Q 2 ) form factor by treating the θ-parameter as a small perturbation [36,[38][39][40][41][42]. Another approach is to extract the CP -odd F 3 (Q 2 ) form factor 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 by replacing the topological charge operator with the flavour-singlet pseudoscalar density employing the axial chiral Ward identities [38,39]. Although this approach provides a well-defined framework, it requires the production of new ensembles at various values of θ and analytic continuation. The cost of simulations of ensembles at the physical point for several values of θ and different values of lattice spacing is currently prohibitively high. Therefore, in this work we use the former approach, keeping θ real and expanding in powers of θ.
Allowing for the CP -violating θ term, the matrix element of the electromagnetic current can be decomposed in four form factors as follows The electromagnetic current is given by J µ e.m. = f e fψf γ µ ψ f , where e f is the electric charge of the quark field ψ f ,ū θ N (p , s ) is the nucleon spinor in the presence of the θ-term and q = p − p is the four-momentum transfer. In the above expression F 1 (Q 2 ) and F 2 (Q 2 ) are the Dirac and Pauli electromagnetic form factors, F 3 (Q 2 ) is the CP -odd form factor and F A (Q 2 ) is the anapole form factor, that vanishes for C-preserving actions. They are all expressed as function of the Euclidean four-momentum transfer squared Q 2 . The presence of θ in Eq. (7) indicates that it is the state with the action that includes the θ-term. Once the dependence of the form factors on the transferred momentum is known, one can extract the electric dipole moment from F 3 (Q 2 ) in the limit Q 2 → 0 using the relation given in Eq. (5).
As pointed out in Ref. [45], however, a spurious mixing between the F 2 (Q 2 ) and F 3 (Q 2 ) form factors occurs if these quantities are not carefully defined. In particular, the contribution to the J µ e.m. matrix elements coming from F 3 (Q 2 ) transforms as an axial 4-vector as expected, only if the spinor u θ N (p, s) appearing in Eq. (7) transforms as a regular Dirac spinor under parity, i.e. only if u θ N (p = (p 0 , − p)) = γ 4 u θ N (p). This holds if spinors satisfy the Dirac equations: where the real-valued m θ N is the mass of the nucleon in the θ = 0 vacuum. On the lattice, the above matrix elements can be extracted from the Euclidean three-point function given by G µ,(θ) where J N ( p f , t f ),J N ( p i , t i ) are the nucleon interpolating operators that respectively create a nucleon at time t i (source) with momentum p i and annihilate it at time t f (sink) and momentum p f . Inserting unity as a complete set of energy and momentum eigenstates in Eq.(10) leads, in the large Euclidean time limit, to the ground state contribution given as Excited states contribution are suppressed in the above expressions. The overlap between the interpolating operators and the nucleon state of a given momentum can be expressed as where the spinorsũ θ N ,ū θ N satisfy the following Dirac equations The imaginary phase in the mass terms that arises due to the CP -violation induced by the θ-term, is parametrized by the so called mixing angle α N (θ). These spinors are related to the u θ N ,ū θ N that are wellbehaved under parity by an axial rotation, namelỹ This can be verified by substituting Eq. (14) in Eq. (13) and noticing that one recovers Eq. (9). With this in mind, one can rewrite Eq. (11) as 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, as follows: while higher order contributions can be neglected, therefore in the following sections we will simply refer to the mixing angle and the CP -odd form factor as α N and F 3 (Q 2 ) correspondingly. The expectation values . . . θ can also be expanded by using Eq.(6). This gives The left hand side (l.h.s.) of Eq. (15) is now written in terms of expectation values of states with θ = 0 and can be computed using gauge configurations extracted with the standard CP -even QCD action. The three-point function reads: One can, thus, divide the right hand side (r.h.s.) of Eq. (15) into a CP -even and a CP -odd part and relate them respectively to the three-point functions of Eq. (19) and Eq. (20). The unknown normalization coefficients Z N can be canceled by taking an appropriate ratio with the two-point function, while the nucleon mixing angle α N can be measured by using the relation and expanding in powers of θ. The result is that the mixing angle α N can be determined from the computation of the following two-point functions

Lattice setup
We use one gauge ensemble of twisted mass fermions produced with 2 degenerate light flavours, a strange and a charm quark (N F = 2 + 1 + 1), with all quark masses tuned close to their physical values. We use the Iwasaki improved gauge action [46], given by where β = 6/g 2 , U 1×1 the plaquette and U 1×2 the rectangular Wilson loops. The Symanzik coefficients are set to c 0 = 3.648 and c 1 = (1−c 0 )/8. The fermionic sector is implemented using the twisted mass formulation of lattice QCD [47,48], which for the degenerate light quark double up and down has the form In the equation above χ (l) is the field representing the light quarks doublet, expressed in the twisted basis, m 0,l and µ l are respectively the untwisted and twisted mass parameters, τ 3 is the third Pauli matrix acting in flavor space and D W is the massless Wilson-Dirac operator. The clover term ∝ σ µν F µν is included in the action to suppress cut-off effects reducing the difference between the mass of the charged and neutral pions [35]. The strange and charm quarks are included as a non-degenerate twisted doublet where m 0,h is the bare untwisted quark mass for the heavy doublet, µ δ the bare twisted mass along the τ 1 direction and µ σ the mass splitting in the τ 3 direction. We tune the partial conserved axial current (PCAC) mass to zero in order to achieve maximal twist. This ensures automatic O(a) improvement for the expectation values of the observables of interest [50]. More details about the generation of this ensemble can be found in Ref. [35]. The lattice size is 64 3 × 128, with lattice spacing a = 0.0801(4)fm, as determined from the nucleon mass [35], pion mass m π = 139(1)MeV and Lm π = 3.62. We will refer to this ensemble as cB211.72.64. For the analysis, we used 750 gauge configurations, separated by 4 trajectories each. The parameters of the simulation are summarized in Table 1.

Correlation functions
For the computation of the nucleon two-and three-point functions we employ the standard proton interpolating field, namely where u(x) and d(x) are up and down quark fields in the physical base, and C = iγ 2 γ 4 is the charge conjugation matrix. 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 Gaussian smeared quark fields [51,52], with 125 smearing steps and parameter α G = 0.2. The smearing is tuned in order to approximately reproduce a mean square radius for the neutron of 0.5 fm. We apply 50 steps of APE-smearing [53] on the gauge links that enter the smearing operator with α APE = 0.5 in order to reduce gauge fluctuations. For the electromagnetic current J µ e.m. (x) we use the symmetrized lattice conserved vector current defined as which, in contrast to the local currentψ f (x)γ µ ψ(x), does not need renormalization. The two-and three-point functions are then given by where the projectors Γ 0 and Γ k are given by For the computation of the three-point functions, we consider only the connected contribution as shown in Fig. 1. We use the standard method of sequential inversions through the sink, taking the final momentum p f = 0. The time slice of the sink relative to the source (sink-source time separation) is kept fixed to t f − t i = 12a. From our previous investigation using an ensemble of N f = 2 + 1 + 1 twisted mass fermions with pion mass m π ≈ 370 MeV we showed that such a time separation is sufficient to suppress excited state contributions to the accuracy of the present study [54]. Larger values of the sink-source time separation, even if desirable for a better suppression of the excited states, lead to large statistical uncertainties, that for the current investigation would require a prohibitively high statistics. We compute the three-point functions using 54 randomly distributed source positions per configuration, that can be treated as statistically independent measures, leading to a total statistics of ∼ 40k data. The two-point functions entering the F 3 determination both explicitly, appearing in the ratio with the three-point function (see Eq. 40 in Section 5.2), and implicitly through the mixing angle α N . In the former case, we employ the same 54 source positions used for the threepoint functions, while for the computation of the mixing angle, we had higher statistics available, namely 200 source positions per configuration. A summary of the statistics used is given in Table 2.

Topological charge
As already mentioned, for the computation of the three-and two-point functions in Eqs. (30),(32) one needs the topological charge. In the continuum, the topological charge is defined as the integral over the four-dimensional volume of the topological charge density given in Eq. (4), namely The discrete counterpart of the above quantity can be obtained by replacing the gluonic field tensor with a lattice operator that reproduces the correct continuum limit. The choice is not unique, and operators with better finite-size effects can be obtained by using O(a)-improved discretizations of F µν . In particular, one of the definitions of Q we used in this work, is the symmetric or 'clover' definition, firstly introduced in Ref. [55]. It has the form and uses a discretization of the gauge strength tensor in terms of a "clover leaf" path C µν , made by the sum of the plaquettes P µν (x) centered in x and with all the possible orientations in the µν-plane, i.e.
This operator is even under parity transformations and exhibits O(a 2 ) discretization effects. We use the gradient flow [56] in order to suppress the UV fluctuations of the gauge field defining the topological charge. The smoothing action employed in the flow equation is the standard Wilson action. The elementary integration step is = 0.01 and the topological charge is computed on the smoothed fields at multiples of ∆τ flow = 0.1. The flow time must be chosen large enough such that discretization effects are canceled but to keep the topological properties of the gauge field unchanged. For this reason we study the dependence of our final quantities on τ flow searching for a plateau region. Accordingly to Ref. [56] we expect that this happens for a √ 8τ flow ∼ O(0.1fm).
In addition, we use a second definition of the topological charge based on spectral projectors as described in Refs. [32,33]. This definition allows one to extract the topological charge from the spectrum of the hermitian Wilson-Dirac operator D † W D W , by employing the relation where u i is the eigenvector related to the i-th eigenvalue λ i , Q 0 is the bare topological charge, and M 0 is the bare spectral threshold. It bounds the modes that enter into the sum in Eq. (37) by requiring λ i < M 2 0 . Accordingly to Ref. [33], the renormalization of these quantities are given by where Z P , Z S are the renormalization constants of the pseudoscalar and scalar densities, respectively. We calculate the lowest 200 eigenvalues of the squared twisted mass Dirac operator using the Implicitly Restarted Lanczos Method (IRLM) where polynomial acceleration is employed. The renormalization constants are computed in a massless renormalization scheme, employing the Rome-Southampton method or the so-called RI scheme [57]. We use five N f = 4 ensembles generated exclusively for the renormalization program at the same β value as that of the ensemble of interest. The five ensembles are generated at different pion masses in the range of 366 MeV to 519 MeV using a lattice volume of 24 3 × 48. This enables us to perform the chiral extrapolation to extract the Z-factors at the chiral limit [27]. For the non-perturbative calculation of the vertex functions we use momentum sources [58] that allow us to reach per mil statistical accuracy with O(10) configurations [59,60]. For the renormalization we need Z P and Z S , which are scheme and scale dependent. We use the MS-scheme, which is commonly used in experimental and phenomenological studies. The conversion procedure is applied on the Z-factors obtained on each initial RI scale (a µ 0 ), with a simultaneous evolution to an MS scale, chosen to be µ=2 GeV. For the conversion and evolution we employ the intermediate Renormalization Group Invariant (RGI) scheme, which is scale independent and connects the Z-factors between the two schemes. For more details about our renormalization program see Refs [26,27,61]. We find that Z P = 0.462(4) [27] and Z S = 0.620(4) 1 .
Using these values of the renormalization constants, and changing the upper limit in Eq. (37) up to the maximum 200 lowest eigenvalues corresponds to a threshold M thr that varies in the range 0 ÷ 65 MeV.
In the rest of the paper, we will refer to the definition of Eq. (35) as "gluonic" or "field theoretic" definition of the topological charge, while the one defined by Eq. (38) will be referred to as the "fermionic" or "spectral projectors" definition.

Topological susceptibility and scale setting
Before presenting results on the CP -odd form factors, we investigate the topological properties of the cB211.72.64 gauge ensemble. In the left panel of Fig. 2 we show the Monte Carlo (MC) history of the topological charge using the two definitions. One can observe that the autocorrelation time of Q is not critical, i.e. no topological freezing occurs at this lattice spacing, and the configurations explore several topological sectors. Moreover, there is a correlation of 72(2)% between results obtained using the gluonic and fermionic definitions. These results are computed at τ flow = 3.5, and M thr = 64.98 MeV but similar observations hold also for other values of the smoothing and cut-off scales.    In the right panel of Fig. 2, histograms of Q are compared. Both definitions lead to a Gaussian-like distribution 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 less fluctuations and, thus, exhibits a narrow distribution. This impacts the values of the topological susceptibility defined as Q 2 /V that differs between the two definitions, as can be seen from Fig. 3. This mismatch is actually expected, due to the fact that the spectral projectors definition shows milder cut-off effects with respect to Q 2 /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 field theoretical definition is essentially flat as a function of the time flow, the one computed using spectral projectors depends strongly on the threshold M thr . The absence of a plateau for the window of eigenvalues computed makes it harder to properly justify a good choice of M thr . In principle, there is no reason to expect to obtain a plateau for a large value of the renormalized cut-off. Actually, the results from Ref. [30] suggest that 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 particular choice of M thr becomes irrelevant once one takes the continuum limit, as long as the renormalized value of the threshold in physical units is kept fixed. Since in this work we use one gauge ensemble, the question that arises is what are the finite-size effects that result from taking a threshold that is not in the plateau region. To give a qualitative answer we use the results obtained in the work of Ref. [29]. In Fig. 4, we show the continuum 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 limit increases when we decrease the threshold, it is still milder than the slope of the susceptibility computed via the 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 on this observation, we use the maximum value of the threshold accessible with the current number of eigenvalues computed, i.e. M thr = 64.98 MeV. We can reasonably expect that similar considerations also hold for the quantities 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.

Nucleon mixing angle
For the determination of the CP -odd form factor F 3 (Q 2 ), one requires the knowledge of the mixing angle α N . We extract it from the following ratio of two-point functions at zero-momentum We take t i = 0 and thus suppress the dependence on t i . We seek for an interval for which the ratio becomes time-independent as a function of t f (plateau region). The time evolution of the ratio of Eq. (39) is illustrated in Fig. 5 as a function of t f /a. We use both the gluonic (left panel) and the fermionic (right panel) definitions for the topological charge that enters in the computation of α N . In particular, the one in the left panel is computed at timeflow τ flow = 3.5, while for the fermionic one, we used a cut-off of M thr = 64.98 MeV.  With the grey band we show the corresponding fits when using for the fit a constant plus an exponential term which takes into account the 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 fermionic definitions, respectively. Fitting to a constant, we find α N = 0.202(20)(4) for the gluonic topological charge and α N = 0.128(9)(3) for the fermionic one. Errors are respectively the statistical and the systematic coming from the choice of the plateau range. The latter is computed by varying the initial time slice in the range [8,12] and the final time slice in the interval [17,20], and taking the largest difference between mean values. For both definitions of Q, it is negligible if compared to the statistical uncertainty. To account for a residual time-dependence we consider the contribution of the first excited state too. As shown in Fig. 5, both Ansätze give 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 topological charge. 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.
This shows a stronger dependence on M thr . Such a behaviour is reminiscent of the one exhibited by the topological charge as discussed in Section 4.1.

Determination of the CP -odd form factor F 3 (Q 2 )
In order to extract the CP -odd form factor F 3 (Q 2 ), we construct an appropriate ratio of the relevant threepoint function given in Eq.(32) and a combination of two-point functions so that unknown overlaps and the exponential time dependence cancel in the large time limit of t f and t ins . This ratio is given by where p f = 0. The form factor F 3 (Q 2 ) is then extracted from where E N is the initial energy of the nucleon, C = (2m 2 N )/(E N (E N + m N )) is a kinematic factor and G E (Q 2 ) = F 1 (Q 2 ) + (q 2 /(2m 2 N ))F 2 (Q 2 ) is the electric Sachs form factor. F 1 (Q 2 ) and F 2 (Q 2 ) are the Pauli and Dirac form factors. We note that Eq. (41) derives from Eq. (55) of Ref. [54] with the correction given in Ref. [45]. The neutron electric form factor G E (Q 2 ) is extracted from Π 00 3pt (see Eq. (A4) of Ref. [28]) We note that Eq. (41) is not the only way to extract the F 3 (Q 2 ) form factor. Similar relations can be written that express the Π kk 3pt,Q and Π j =k 3pt,Q components as linear combinations of F 1 , F 2 and F 3 . We verified that they lead to worse signal-to-noise ratios, and thus we opt to use Eq. (41) in this work. In Fig. 7, we show the ratio defined in Eq. (40) as a function of the insertion time t ins at fixed sinksource time separation t f = 12a, for the smallest three non-zero values of the momentum transfer squared, i.e. Q 2 = 0.056 GeV 2 , Q 2 = 0.111 GeV 2 and Q 2 = 0.164 GeV 2 . The results shown for Π 0k 3pt,Q are aver-aged among momenta with non-zero k-component in all k-directions. Despite the large statistics employed (∼ 40000 samples), the errors are large and do not allow to increase further the sink-source time separation. We fit the ratio in symmetric intervals [−t fit , t fit ] and vary the fit ranges, taking t fit = 2, 3, 4. The resulting values are all compatible within our current statistical accuracy. It is worth to be noting that the results extracted using the fermionic definition of the topological charge show a significant error reduction with respect to the gluonic counterpart. Nevertheless, a zero value cannot be excluded for all three momentum transfers. This is made even more evident in Fig. 8, where the values of F 3 at different Q 2 are reported. In order to extract F 3 (0) and thus d θ N , we extrapolate to Q 2 = 0 by considering the weighted average of the values at the three smallest Q 2 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 final results are field theoretical or gluonic definition fermionic definition via spectral projectors |d θ N | = 0.0009(24) θ e · fm .  As can be seen, the fermionic definition based on spectral projectors of the topological charge provides a more accurate determination than the gluonic definition. Considering the absolute error as a bound for the magnitude of the nEDM, what we have is an improvement by a factor ∼ 2 when using the fermionic definition. Therefore, the additional cost coming from the computation of the eigenmodes for the fermionic definition of 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 in the computation of the results presented, is in some sense the most conservative choice. Indeed, by looking at the dependence of the extrapolated value of F 3 (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 increases with increasing M thr . This is in contrast with what is observed for the mixing angle α N , that instead changes significantly with M thr . This is because the changes in α N are small as compared to the uncertainty of Π 0k 3pt,Q and are thus not visible in F 3 at the current level of precision.  This work, a = 0.0801(4) fm  Table III of Ref. [45], where the spurious contribution coming from F 2 (Q 2 ) is subtracted. See Ref. [45] for further details.

Comparison with other determinations of nEDM
In Fig. 10 we provide a comparison of our result with those of other lattice QCD studies for a similar lattice spacing. We observe that our value has a statistical error that is comparable with the recent results in Ref. [36] that were, however, computed at much larger values of the pion mass. Since the errors grow with decreasing pion mass and so does the computational cost, achieving such an accuracy it is a major outcome of this work. We note that the authors of Ref. [36] using their results at these heavy pion masses perform a chiral extrapolation and find a non-zero value of |d θ N | = 0.00152(71)θ e · fm at the physical point. The accuracy of their determination is due to the chiral extrapolation where systematic errors from using chiral expressions cannot be determined. It is worth mentioning that their actual data have uncertainties similar to this work. Our result computed directly at the physical point does not exclude a zero value, but provide a significant bound to the value of nEDM.

Summary and Conclusions
We compute the neutron electric dipole moment using an ensemble of N f = 2 + 1 + 1 twisted mass fermions simulated at the physical point and with a lattice spacing of a 0.08 fm. The extraction of the CP -violating form factor F 3 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 definition by means of gradient flow and spectral projectors, respectively. This enables us to test what effect different definitions of the topological charge may have on the nEDM. F 3 (Q 2 ) cannot be extracted directly from the nucleon matrix element due to the presence of momentum appearing multiplicatively in front of the form factor. The usual approach is to compute F 3 (Q 2 ) for finite Q 2 and then extrapolate to Q 2 = 0 using some fitting form. We would like to highlight the following three crucial aspects of our investigation of the CP -odd form factor: • Determination of nEDM directly at the physical point: We perform the computation directly at the physical value of the pion avoiding uncontrolled errors that a chiral extrapolation may result in. We study F 3 (Q 2 ) as a function of the the momentum transfer Q 2 , and show that there is no noticeable dependence on low values of Q 2 within our current statistical accuracy. We, thus, perform a weighted average on the three lowest ones to extract the value of the form factors at Q 2 = 0.
• Definition of the topological charge: The approach is based on using spectral projectors to determine the topological susceptibility yielding a statistical error that is twice smaller than when using a gluonic definition whether using gradient flow, link smearing, or cooling are applied [29,30]. Furthermore, it has been shown that the topological susceptibility calculated using the spectral projectors approach shows milder lattice artifacts than when using a gluonic definition [29]. Both these observations have led us to investigate the approach based on spectral projectors for the definition of the topological charge entering the evaluation of the nEDM since we expected that these features would also hold by secondary quantities involving the topological charge. The necessity to reduce as much as possible the lattice cut-off effects stem from the fact that currently we do not have three ensembles with different lattice spacings at the physical point to enable us to take the continuum limit. Therefore, while F 3 (0) extracted using the fermionic definition is consistent with the value extracted using the gluonic definition, the statistical errors of the former are about half as compared to those when using the gluonic 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 the topological charge. We find |d θ N | = 0.0009(24) θ 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 observe a clear plateau for α N when using spectral projectors as a function of M thr . This is an artifact of carrying out the investigation at one lattice spacing. Only after a continuum extrapolation keeping the renormalised cut-off M thr in physical units fixed one can extract physical observables that are independent of the cut-off. The same holds for the field theoretic definition for which the gradient flow time in physical units must be fixed as one takes the continuum limit. This is illustrated in Fig. 3 for the susceptibility, where only after the continuum limit is taken there is an agreement among the approaches.
A question arises as of why when using the gradient flow with the Wilson smoothing action we obtain quantities 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. On 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 topological content which remains unaltered until the smoother starts affecting the large instantons contributing to the topological structure of the theory. On the other hand, the same quantities extracted using the spectral projectors reveal a behaviour that depends much stronger on the associated cut-off. This is because, for a theory that breaks chiral symmetry, as we sum up eigenvalues of the squared Dirac Matrix, the summation will keep increasing. This observation is fully supported by the results on the topological susceptibility demonstrated in Ref. [29]. One may argue that the existence of a plateau enables to reliably quote a value of a quantity at a given lattice spacing; this is of course true but has no physical importance since in the end a continuum extrapolation is needed in order to get rid of discretization effects. Nevertheless, we observe that, when using spectral projectors, F 3 appears to exhibit 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 difficulty of this particular method to deliver a statistically significant result. Practically, one needs to increase measurements of at least an order of magnitude to reduce the error significantly hopping to exclude a zero value. Employment of techniques, such as volume clustering [24,64] or spectral projectors as a definition 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 a dynamical θ-term one may improve the importance sampling for the nEDM signal inducing nonzero average topological charge. Exploratory investigations using this method with a field theoretic definition of the topological charge are under way.