Neutron Electric Dipole Moment on the Lattice

For the neutron to have an electric dipole moment (EDM), the theory of nature must have T, or equivalently CP, violation. Neutron EDM is a very good probe of novel CP violation in beyond the standard model physics. To leverage the connection between measured neutron EDM and novel mechanism of CP violation, one requires the calculation of matrix elements for CP violating operators, for which lattice QCD provides a first principle method. In this paper, we review the status of recent lattice QCD calculations of the contributions of the QCD $\Theta$-term, the quark EDM term, and the quark chromo-EDM term to the neutron EDM.


Introduction
Electric dipole moment (EDM) of the neutron measures the separation of positive and negative charge in the neutron and is necessarily aligned along the spin axis. For a neutron to have an EDM, the theory of elementary particles must violate parity (P) and time reversal (T) invariance, and charge conjugation and parity (CP) symmetry if CPT is conserved. Under a parity transformation, the EDM of any particle changes its sign while the direction of the spin remains the same. Under time-reversal transformation, the spin flips its direction while the EDM stays unchanged.
So far, a nonvanishing neutron EDM (nEDM) has not been observed, and the current experimental upper bound is 3.0×10 −26 e·cm [1]. In recent experiments, the nEDM is measured by the change in the spin precession frequency of ultracold-neutrons aligned in a magnetic field under a flip in the direction of a strong background electric field. Figure 1 shows the history of the experimental upper bound on the neutron EDM, as well as the target precision, ∼ 5 × 10 −28 e·cm [2][3][4][5][6], of proposed experiments.
In the standard model, the leading contribution to nEDM due to the CP violating phase in the CKM matrix comes from three-loop and higher order diagrams, and the expected size is more than 5 orders of magnitude below the current experimental bound [7]. In extensions of the standard models, however, nEDM can appear at one-loop due to novel CP violating interactions. In some of the most popular models, such as supersymmetric (SUSY) models, the expected size of the nEDM is between 10 −25 -10 −28 e·cm as shown in Figure 1. Therefore, a nEDM smaller that 10 −28 e·cm will put a serious constraint on those models [8][9][10].
There are two main outcomes of neutron EDM experiments for BSM physics. First, a non-zero measurement of nEDM would establish new sources of CP violation, which is one of Sakharov's three conditions for weak-scale baryogenesis [11]. The CP violation in the standard model is not sufficient to explain observed the baryon asymmetry of the universe, and it requires new sources of In the standard model, the leading contribution to nEDM from the CP violating phase in the CKM matrix arises at three-loops, and the expected size of nEDM is more than 5 orders of magnitude below the current experimental bound, as shown by the green shaded area in the plot [7]. In extensions of the standard models, however, nEDM can appear at one-loop with the new CP violating interactions. In some of the most popular models, such as SUSY, the expected size of nEDM covers the regions of the current and future experimental bound, as shown by the red shaded area [8][9][10].
CP violation from beyond the standard model (BSM) [12]. Neutron EDM is a good probe for such novel CP violation. Second, most BSM theories have additional sources of CP violation. There are many BSM scenarios predicting a nEDM between 10 −25 -10 −28 e·cm, and upcoming experiments will put constraints on them [2][3][4][5][6]. This requires both the measurement of, or a bound on, the neutron EDM, and the calculations of the matrix elements of novel CP violating operators within the neutron states. Current non-lattice estimates of the matrix elements have large uncertainties [10]. Lattice QCD can provide a first principle, non-perturbative method for calculating these matrix elements.
To analyze novel CP violation, we work within the framework of effective field theories and classify operators by their canonical dimension [13]. At hadronic scale, the effective Lagrangian for the CP violating interactions at dimension 4, 5 and 6 can be written as where g is the QCD coupling constant, and G µν,b = ε µναβ G b αβ /2. The first term on the right hand side (r.h.s) is the dimension-4 Θ-term already allowed in QCD, and the next two terms are the dimension-5 quark EDM (qEDM) and quark chromo-EDM (cEDM) terms. There are two types of dimension-6 terms in the second line in the above equation: the Weinberg three-gluon operator and various fourquark operators.
If the matrix elements of all the operators are O(1) and the anomalous dimensions are small, then the O(1/Λ D−4 scale ) suppression implies that the lower mass dimension operators are more important. However, it turns out that the current bound on the nEDM already makes the dimension-4 QCD Θterm unnaturally tiny;θ ≤ O(10 −9 − 10 −11 ) [13][14][15][16], whereθ is the total effective CP violating anglē θ = θ + arg det M CKM with the quark mixing matrix M CKM . The unexpected smallness of this number is known as the strong CP problem. One of the popular models explaining the strong CP problem is the Peccei-Quinn mechanism [17,18], which requires the existence of the, as yet unobserved, axion. The two dimension-5 operators typically arise at the TeV scale as dimension-6 operators involving a Higgs. At the hadronic scale, 1 GeV, the Higgs fields are replaced by its vacuum expectation value, v EM , and the operator becomes dimension-5 with the coefficients suppressed by v EM /M 2 BSM . How well these dimensional arguments work and which terms in Eq. (1) make the largest contribution to nEDM depends on the details of the BSM theory. But, while the couplings depend on the starting BSM model, the matrix elements of these operators within the neutron states are model independent results. Thus, lattice QCD calculation of the matrix elements of these operators will play an important role in connecting the measured neutron EDM and novel CP violation in BSM scenarios, i.e., knowing the matrix elements along with the bound on (or the value of) the nEDM, one can put bounds on the couplings and thus on the parameter space of allowed BSM theories.
In this paper, we review recent progress in Lattice QCD calculations of the matrix elements of these CP violating operators. We start by discussing the contribution of the QCD Θ-term to the nEDM and use it to illustrate the various methods used in the calculations.

QCD Θ-term
There are three strategies for the calculation of the contribution of CP violating operators to the neutron EDM on the lattice. We illustrate these three approaches for the case of the CP-violating Θ-term whose volume integral is the topological charge Q top . They are the external electric field method [19][20][21][22], expansion in θ [23][24][25][26][27], and simulation with imaginary θ [28-30].

External Electric Field Method
The first method uses an external electric field [19][20][21][22], in which the neutron EDM, d N , is extracted from the energy difference of the neutron states with spin S aligned [anti]parallel to the external electric field E: For a given value of θ, the effect of the Θ-term is included by reweighting the nucleon correlation function with the topological charge: where the expectation value on the r.h.s is evaluated on lattice configurations generated without the Θ-term. This method extracts the neutron EDM from a spectral quantity, not the form factor F 3 as discussed next, so the results are not affected by the mixing under parity violation problem discussed in Section 3.

Expansion in θ
Based on the current experimental bound of the neutron EDM and model studies, we know that the coupling θ is tiny. Hence it can be treated as a small expansion parameter. Any expectation value in presence of a small QCD Θ-term can then be written as follows [23][24][25][26][27]: In the second line of the above equation, the QCD Θ-term effect is included at the lowest order by measuring the correlation of the topological charge with the observable O, using configurations generated without the Θ-term.
To calculate the nEDM, the operator O is chosen as the matrix element of the vector current within the neutron states. Assuming PT invariance, the matrix element, calculated as a function of q µ , the momentum transfered by the vector current, can be decomposed in terms of Lorentz covariant form factors: Of these, the form factor of interest is F 3 at zero-momentum transfer. This is extracted by extrapolating F 3 (q 2 ), measured at finite q 2 , because the term containing F 3 only contributes at q ν 0 as can be seen from Eq. (5). Then the contribution to d N is given by This approach, which relies on calculating the form-factor F 3 from nucleon matrix elements, requires a careful handling of the neutron spinors in a theory with P violation. Most previous calculations did not correctly account for this subtlety [31]. The problem can be corrected retroactively, and in Figure 3 we show the resulting reduction in the value reported in recent calculations [26,27]).

Imaginary θ Simulation
The Θ-term in Euclidean space-time is purely imaginary, and Monte Carlo simulations including the Θ-term face the sign problem. One way to circumvent this issue is to make the action real by taking θ to be purely imaginary [28][29][30]. Assuming that the theory is analytic around θ = 0, the results can then be continued to real θ for small θ. Furthermore, this calculation can be carried out using the axial anomaly whereby the QCD Θ-term is chirally rotated to the fermionic term S q θ as Lattice simulation can then be performed by adding S q θ to the original QCD action. In this approach, the EDM is again extracted from the form factor analysis, so the results are affected by the mixing induced by parity violation discussed in Section 3. Again, in Figure 3 we show one of the recent results [30] with and without the correction due to this mixing problem.  Figure 2. CP violating phase α (denoted as α 1 in the figure) as a function of the cutoff radius R presented in Ref. [32]. The central value of the α appears to saturate after R is larger than about 15 lattice spacings, which corresponds to 1.71 fm in physical unit, while the statistical error grows as R is increased.

Variance Reduction using Cluster Decomposition
In the calculation of neutron EDM induced by the QCD Θ-term using θ-expansion or external electric field method, one needs to calculate correlators reweighted by the topological charge. Recently, the authors of Ref. [32] reported an error reduction technique using cluster decomposition. The main idea is to calculate the topological charge only in the vicinity of the sink in the momentum projection summation. Let us consider a neutron two-point correlator reweighted with the topological charge Q: where G is the source grid. In conventional calculations, Q is calculated by summing the topological charge density q(y) from all lattice points. In the new formula, however, the topological charge is calculated only up to a certain radius R from the sink point under the assumption that correlations due to long-distance points are negligible. Figure 2 presents the CP violating phase, labeled α 1 and calculated using the Eq. (8), as a function of the cutoff radius R. Data show that the mean value of α 1 saturates at some reasonably small value of R. Including the contributions of the topological density from larger R only increases the statistical noise. Considering that the mean value still exhibits significant fluctuations, a more detailed costbenefit analysis of the approximation is required.

Extracting F in a Theory with Parity Violation
In a theory with P and CP symmetry, the spinor u of the neutron state satisfies the following Dirac equation: with γ 4 the parity operator: u p → γ 4 u − p under P. When CP, but not PT, is violated, the Dirac equation is only modified by a phase factor to where the spinor that solves the modified equation is defined to beũ. Furthermore, for this equation, γ 4 is no longer the parity operator, rather under parityũ p → e −2iαγ 5 γ 4ũ− p . The two solutions u andũ are related asũ When calculating the form factors F 2 and F 3 , it is important to enforce that F 2 is the parity-even magnetic form factor while F 3 is the P and CP odd electric form factor. There are two ways to enforce these symmetry requirements. The first is to properly include the phase defined in Eq. (11) in the definition of the matrix element. In that case the decomposition into the form factors is the same as in a theory with CP symmetry with unmodified spinors u. Alternatively, one can calculate the standard matrix element, and then undo the mixing between F 2 and F 3 due to the non-standard parity transformation with phase α. Then, the two form factors F 2 and F 3 are given in terms ofF 2 andF 3 as in Eq. (14).
In both cases, the phase α has to be extracted from the nucleon 2-point function. In this extraction, it is important to note that the phase α is state dependent. In Section 5, we show that the α corresponding to the ground state is given by the long-time behavior of the nucleon 2-point function, i.e., when the ground state dominates.
In the first approach, one includes the phase in the calculation of the n-point functions: The resulting F 3 (0)/2M N is the desired contribution to the nEDM. In the second approach, one calculates and one has to extract F 3 from the following mixing structure: The two approaches are equivalent.
Unfortunately, as pointed out in Ref. [31], all calculations based on the F 3 extraction prior to their work, starting with the 2005 work in Ref. [23], used the second approach but did not carry out the rotation defined in Eq. (14) to get the F 3 . Their quoted results are, therefore, forF 3 . The authors of Ref. [31] show that correcting for this omission reduces the value of an already hard to measure F 3 by about a factor of ten. Thus, all previous estimates of the contribution of the Θ-term and the cEDM based on evaluating F 3 need to be revised. Figure 3 shows some of the recent lattice results of F 3 induced by the QCD Θ-term, before and after the parity mixing correction. The interesting point is that all corrected numbers are close to zero. Previous phenomenological estimates were an order of magnitude smaller than lattice QCD results. This tension may disappear after a proper analysis of the phase α introduced in the neutron states due to the breaking of parity. Further discussion of the mixing between F 2 and F 3 when parity is violated can be found in Ref. [31].  [31,33]. As noted in Ref. [31], some assumptions are made to calculate the corrected F 3 numbers as some of the original papers do not contain full information needed for the parity induced rotation Eq. (14), so the F 3 results presented in this plot may not be precise.

quark EDM
The electromagnetic current, V µ ≡ δL/δA µ , gets an additional contribution when the quark EDM operator is added to L SM . Thus, while the standard conserved vector current interacts with a quark with charge q f , this additional contribution, which is is just the flavor diagonal tensor operator for each quark flavor, couples with strength d γ f g f T , where d γ f are the BSM model dependent CP violating couplings and g u T , g s T , g s T , . . . are the tensor charges, i.e., matrix elements of the quark bilinear tensor operator within the nucleon states. Furthermore, the effects analogous to the mixing discussed in Section 3 are suppressed by powers of the electromagnetic coupling α EM ∼ 1/137. Thus, the leading contribution of the EDM of the quarks to the nEDM is given by the flavor diagonal tensor charges: In many models, such as the supersymmetric models, the quark EDMs are proportional to the corresponding quark masses (d q ∝ m q ), because in these theories all connections between left and righthanded quarks are mediated by a common set of Yukawa interactions. Since the strange quark is much heavier than the u and the d quarks, the strange quark contribution gets enhanced, by m s /m d ≈ 20. On the other hand g s T g l T . Since the contribution is the product of the two, it is, therefore, important to determine the strange quark, and perhaps even the charm quark, tensor charge precisely.
There are two classes of diagrams in the calculation of flavor diagonal tensor charges: quarkline connected and the quark-line disconnected diagrams as illustrated in Figure 4. In case of the strange (charm) quark tensor charge, only the disconnected diagram with strange (charm) quark loop contributes. This is the basis of the expectation that g s T g l T . In Refs. [34,35], the tensor charges were calculated using clover fermions on the HISQ lattices, and extrapolated to continuum and physical pion mass limit, as illustrated in Figure 5 (left). The dis-  Figure 5. (Left) Extrapolation of g s T to the continuum limit. The fit is performed simultaneously in a and M π but only the projection to the a-plane with data extrapolated to the physical M π , is presented. (Right) Bounds on d u and d d are obtained using the formula given in Eq. (15), the tensor charges given in Eq. (16) and the experimental bounds |d N | < 2.9 × 10 −26 e·cm. Since the strange quark tensor charge is consistent with zero within statistics, it imposes no constraint on d s . Both figures are taken from Ref. [34]. connected diagrams were estimated using a stochastic method, and it turned out that the disconnected diagram contribution to the tensor charges is very small. Their results are (7), g s T = 0.008 (9) .
There are BSM models in which the quark EDMs are the dominant sources of CP violation at low energy. In these scenarios, constraints can be placed on the couplings of the EDM of each quark flavor using the experimental bound on the neutron EDM and the lattice results for the tensor charges. Neglecting all other sources of CP violation, figure 5 (right) shows the allowed region of d u and d d using the lattice results, Eq. (16), described in Ref. [34,35]. Of all the CP violating operators listed in Eq. (1), the calculation of the quark EDM, including renormalization, is theoretically wellestablished and the results in Eq. (16) for g u T and g d T indicate that these are already available at the 15% level. Further improvement in precision will occur as higher precision data are generated in future calculations.

Quark Chromo-EDM
The last dimension-5 term is the quark chromo-electric dipole moment (cEDM) operator. The modification to the action in presence of the CP violating cEDM term is: whered q are again BSM couplings evolved to the hadronic scale for each quark flavor. The lattice calculation consists of evaluating the matrix element of the insertion of the product of the electromagnetic current V µ and theq(σ · G)γ 5 q operator within the neutron states. Then, the contribution of cEDM operators to the nEDM, ignoring the complicated operator mixing problem under renormalization discussed in Section 5.4 below, is obtained from qdq N|V µq (σ · G)γ 5 q|N . Lattice QCD studies of the cEDM operator have started only recently, and three methods are being explored: an expansion ind q [31], external electric field method [31], and the Schwinger source method [36,37].

Expansion ind q
The expansion ind q method proceeds analogously to the expansion in θ in the QCD Θ-term calculation. For small BSM cEDM couplings,d q , the neutron matrix element of the electromagnetic current can be written in terms of the expectation values evaluated on the standard CP even lattices, Their contribution to the neutron EDM is extracted from the CP violating form factor, F 3 as defined in Eq. (5). The calculation is challenging because it requires calculating the expectation value of fourpoint functions because one has to insert both the vector current and the cEDM operator between the neutron source and sink, with the cEDM operator defined as a sum over all space-time points. Figure 6 (top) shows all the required quark-line connected diagrams. These four-point correlators are constructed using five types of propagators used as building blocks that are shown in the bottom panel of Figure 6 (bottom). F and B are the usual forward and backward propagators. C is the cEDMsequential propagator. It is constructed by starting with the regular quark propagator and inserting the chromo-EDM operator at all lattice sinks points. This is then used as the source for calculating the chromo-EDM inserted propagator.
Results from Ref. [31] for the signal in F 3 and the source-sink separation dependence of the signal are shown in Figure 7 (left). F 3 has larger value when the cEDM is inserted in the d-quarks of the neutron. In Figure 7 (right), the results are extrapolated to zero-momentum to obtain the contribution F 3 (0) as needed for the neutron EDM.

External Electric Field Method
The chromo-EDM contribution to the neutron EDM also can be calculated using the external electric field. Similarly to the QCD Θ-term, the neutron EDM is extracted from the energy difference of nucleon states under spin flip in presence of a uniform external electric field iE: where the neutron correlators in the presence of the cEDM term are obtained by reweighting, This method needs only two-and three-point functions of the neutron, the latter with the insertion of the cEDM bilinear operator. In Ref. [31], F 3 is extracted using the external electric field method on the same ensemble of lattices used for the calculation using thed q -expansion method, and the results are presented in Figure 7 (right). The signal is poorer in the external electric field method compared to those in expansion ind q method, but the results of the two methods are consistent. Note that the external electric field method is not affected by the parity mixing problem described in Section 3, while thed q -expansion method is. Therefore, the consistency of the results from the two methods, after taking care of the parity mixing rotation given in (14), suggests that both methods are yielding a reliable signal.

Schwinger Source Method
Another method to calculate the contribution of the cEDM operator to the neutron EDM is the Schwinger source method [36,37]. Noting that the cEDM operator is a quark bilinear, iq(σ · G)γ 5 q, one can add it to the QCD fermion action: where D clov is the Dirac operator for the Wilson-clover action, and ε is a small control parameter the dependence on which will be removed by taking a derivative with respect to it at the end of the calculation. Because cEDM is a quark bilinear, the integration over the fermion degrees of freedom in the path integral can still be carried out as before. In the case of the clover fermions, the addition of the cEDM operator is equivalent to shifting the coefficient of the dimension-5 clover term by iεγ 5 : In this approach, the insertion of the cEDM term is carried out by calculating valence quark propagators including this modified clover term, and using this propagator in all diagrams requiring a cEDM insertion. Having taken care of the cEDM term, the calculation of the original 4-point functions reduces to three-point functions, with an insertion of the vector current within the neutron state, made up of propagators with and without cEDM insertion.
To perform this calculation on ensembles generated with just the standard QCD action, however, requires taking into account the change in the Boltzmann factor under the addition of the cEDM term to the action as shown in Eq. (22). This, a priori, non-unitary formulation between sea and valence quarks can be accounted for at leading order in ε by reweighting each configuration by ratio of the fermion determinant with and without the cEDM term: Diagrammatically, all the terms that need to be calculated in this Schwinger source method are shown in Figure 8. There are three classes of diagrams at leading order in ε: the reweighting factor for each configuration, and the quark-line connected and disconnected diagrams on that configuration. Current calculations are at the stage of demonstrating a signal in the full evaluation of Figure 8 with both the cEDM and the γ 5 operator with which it mixes under renormalization as discussed in Section 5.4. First, in Figure 9, we show the success at extraction of the phase α that is induced in the ground state nucleon spinor by the CP violating cEDM and γ 5 operators. In the right panel of Figure 9, we also show that this α is linear in ε, which allows us to tune the value of ε. We want to use a large ε to increase the signal but remain in the linear response regime.
First examples of the quality of the signal in the contribution of the connected diagrams to the form factor F 3 induced by the cEDM and γ 5 operators are shown in Figures 10 for two different source and sink separations. At present, the signal is consistent with zero for both operators. Red crossed box is the cEDM operator insertion, P is the regular quark propagator without the cEDM operator insertion, and P ε is the quark propagator with the cEDM operator insertion. An identical calculation has to be done with the cEDM operator replaced by the γ 5 operator in order to define a finite renormalized quantity.

Renormalization of cEDM Operator
The renormalization of the cEDM operator is studied both in one-loop perturbation theory with Twisted-mass fermions [38], and in nonperturbative RI-SMOM scheme [39]. The most challenging outcome is the divergent mixing with lower-dimensional operators: in particular the 1/a 2 mixing with the pseudoscalar quark bilinear operator: Because the mixing is 1/a 2 divergent, it needs to be calculated precisely so that a finite operator can be constructed, ensemble by ensemble, by subtracting the two terms. In Ref. [38], the authors presented the results of the 1/a 2 mixing coefficients for the Twisted-mass fermions calculated nonperturbatively and using the one-loop perturbation theory. In Ref. [39], authors defined a momentum subtraction scheme, RI-SMOM, for the non-perturbative renormalization of the cEDM operator on the lattice, and provided one-loop matching coefficients to the MS scheme in continuum limit.

Summary
In this talk, we reviewed the recent lattice QCD calculations of the contribution of dimension four and five CP violating operators to the EDM of the neutron.
We reviewed the three approaches, external electric field method, expansion in θ, and simulation with imaginary θ, used for the calculation of the QCD Θ-term in Section 2. We also summarized the observation made in Ref. [31] concerning the subtlety of properly defining the neutron spinor in the presence of parity violating interactions and its consequences vis-a-vis mixing between the form factors F 2 and F 3 in Section 3. All calculations based on extracting F 3 , previous to Ref. [31], had missed this mixing. After correction, current lattice QCD results for the QCD Θ-term, are reduced by about a factor of ten and do not show a statistically significant non-zero signal.
The quark EDM contribution to the neutron EDM can be written in terms of the neutron tensor charges g T q . Currently, the g T u,d are determined to within 15% uncertainty including the quark-line disconnected diagrams, and g T s is known to be very small. The results are summarized in Section 4. Lattice QCD calculation for the quark chromo-EDM operator have just started. Three methods for the calculation of the contribution of cEDM and the construction of a finite renormalized operator are described in Section 5. The methodology for calculation of the 3-point (or 4-point) correlation functions is under control but a non-zero signal in all the diagrams has yet to be demonstrated. Once a signal is achieved, we will still be a long way away from a full calculation due to the divergent mixing between the cEDM and γ 5 operators.
A few brief words on other lattice QCD calculations related to the neutron EDM that have not been covered in this review. First, preliminary results for the Weinberg three-gluon operator, one of the dimension-six operators, were presented at the Lattice 2017 conference [40]. Calculations of the Weinberg operator are at the stage of demonstrating a signal and have not even begun to address the issue of renormalization that is potentially even more complicated than that for the cEDM operator. Calculations exploring the four-fermion dimension six operators are yet to begin. Second, there are efforts to study the neutron EDM using lattice spectroscopy, and other matrix elements combined with chiral perturbation theory [41], in addition to the direct calculation of the matrix elements with the CP violating operators.
To summarize, the theoretical calculations of matrix elements needed to extend the importance and reach of nEDM experiments to constrain BSM physics have begun in earnest but are very challenging and will require many new innovations in both theory and computations.