On Lattice Calculation of Electric Dipole Moments and Form Factors of the Nucleon
M. Abramczyk S. Aoki T. Blum T. Izubuchi H. Ohki S. Syritsyn





Physics Department, University of Connecticut, Storrs, CT 06269, USA


Center for Gravitational Physics, Yukawa Institute for Theoretical Physics,
  Kyoto University, Kyoto 606-8502, Japan


Physics Department, University of Connecticut, Storrs, CT 06269, USA
RIKEN/BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA


RIKEN/BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA


RIKEN/BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA


Jefferson Lab, 12000 Jefferson Ave, Newport News, VA 23606, USA
Kavli Institute for Theoretical Physics, UC Santa Barbara, CA 93106, USA
Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA

11.15.Ha, 12.38.Gc, 12.38.Aw, 21.60.De
RBRC-1226





We analyze commonly used expressions for computing the nucleon electric 
dipole form factors (EDFF) F_3 and moments (EDM) on a lattice
and find that they lead to spurious contributions from the Pauli form factor F_2 
due to inadequate definition of these form factors when parity mixing of lattice nucleon 
fields is involved.
Using chirally symmetric domain wall fermions, we calculate the proton and the neutron EDFF
induced by the CP-violating quark chromo-EDM interaction using the corrected expression.
In addition, we calculate the electric dipole moment of the neutron using background
electric field that respects time translation invariance and boundary conditions,
and find that it decidedly agrees with the new formula but not the old formula for F_3.
Finally, we analyze some selected lattice results for the nucleon EDM and observe that after 
the correction is applied, they either agree with zero or 
are substantially reduced in magnitude, thus reconciling their difference from 
phenomenological estimates of the nucleon EDM.


On Lattice Calculation of Electric Dipole Moments and Form Factors of the Nucleon
    S. Syritsyn
    18 november 2016
=================================================================================








§ INTRODUCTION



The origin of nuclear matter can be traced back to the excess of nucleons 
over antinucleons in the early Universe and it is one of the greatest puzzles
in Physics known as the baryonic asymmetry of the Universe (BAU).
One of the required conditions for the BAU is violation of the  symmetry ().
In the Standard Model (SM), the CKM matrix phases lead to  violations in weak interactions, 
but their magnitudes are not sufficient to explain the BAU, 
and signs of additional  are actively sought in experiments.
The most promising ways to look for  are measurements of
electric dipole moments (EDM) of atoms, nucleons and nuclei.
In particular, the Standard Model prediction for the neutron EDM is five orders 
below the current experimental bound, and represents a negligible background.
Near-future EDM experiments plan to improve this bound by 2 orders of magnitude,
and are capable of constraining various Beyond-the-Standard-Model (BSM) extensions 
of particle physics, purely from low-energy nuclear and atomic high-precision experiments.


Knowledge of nucleon structure and interactions is necessary to interpret these 
experiments in terms of quark and gluon effective operators and put constraints 
on proposed extensions of the Standard Model, in particular SUSY and GUT models 
as sources of additional . 
Connecting the quark- and gluon- to hadron-level effective  
interactions is an urgent task for lattice QCD (an extensive
review of EDM phenomenology can be found in Ref. <cit.>).



The proton and the neutron can have electric dipole moments only if the symmetry of 
the Standard Model Lagrangian is broken by additional P-,T-odd interactions.
The only such dimension-4 operator is the QCD θ̅-angle 
(θ̅ stands for the physically-relevant combination of the QCD θ 
angle and quark mass phases).
The θ̅-induced nucleon EDMs (nEDMs) have been calculated on a lattice 
from energy shifts in uniform background electric 
field <cit.>
or extracting the P-odd electric dipole form factor (EDFF) F_3(Q^2)
from nucleon matrix elements of the vector current in  vacuum <cit.>.
Nucleon EDMs 

have been studied using QCD 
sum rules, quark models, and chiral perturbation theory (see Refs. <cit.> 
to name a few).
On a lattice, quark EDM-induced nucleon EDMs have been recently
computed on a lattice in partially-quenched framework <cit.>.
Another important dimension-5(6)[
  These operators are sometimes referred to as “dimension-6” because in the SM they
  contain a factor of the Higgs field.] 
operator is the -odd quark-gluon interaction,
also known as the chromo-electric dipole moment (cEDM)

    _cEDM = i ∑_ψ=u,dδ̃_ψ/2ψ̅(T^a G^q_μν) σ^μνγ_5 ψ ,

and calculations of cEDM-induced nEDMs have recently started using Wilson
fermions <cit.>.





In this paper, we report several important achievements in studying nucleon EDMs 
on a lattice. 
First, we argue that the commonly accepted methodology for computing electric 
dipole form factors of spin-1/2 particles on a lattice has a problem to identify the electric
dipole moment form factor.
In particular, in the standard analysis of the nucleon-current correlators <cit.>, the electric dipole form factor 
F_3 receives large and likely dominant contribution from spurious mixing with 
the Pauli form factor F_2.
The energy shift methods <cit.> are not affected by
such mixing, but their precision has not been sufficient to detect the discrepancy.
This problem affects all the previous lattice calculations of the nucleon EDFFs and EDMs from
nucleon-current correlators, including those studying the θ̅-angle <cit.> 
as well as the more recent studying chromo-EDM <cit.>.
We demonstrate the problem formally in Sec <ref>
and also derive the correction for the results of Refs. <cit.> to subtract the spurious mixing with F_2.
In addition, in Sec. <ref> we study the energy shift of a neutral
particle on a Euclidean lattice in uniform background electric field.
We introduce the uniform electric field preserving translational invariance and periodic 
boundary conditions on a lattice <cit.>.
In order to satisfy these conditions, the electric field value has to be analytically continued to 
the imaginary axis upon Wick rotation from Minkowski to Euclid, and we demonstrate that 
the eigenstates of a fermion having an EDM are shifted by a purely imaginary value.
In Sec. <ref>, we apply this formalism to the analysis of neutron correlators 
computed in the presence of the quark chromo-EDM interaction (<ref>).


Calculation of the neutron EDM in the background field is independent from parity mixing
ambiguities, and it allows us to validate our new formula for the EDFF F_3 numerically.
The difference is evident only if the nucleon “parity-mixing” angle  is large, 
≳1.
The calculations with quark chromo-EDM generate very strong parity mixing compared to the 
θ̅-angle, which is beneficial for our numerical check.
In Section <ref> we calculate the proton and neutron EDFFs F_3p,n(Q^2) 
induced by the quark chromo-EDM interaction (<ref>), as well as the regular 
-even Dirac and Pauli form factors F_1,2.
In Sec. <ref> we compare the EDM results from the form factor and 
the energy-shift calculations, 
providing a numerical confirmation of the validity of our new EDFF analysis.

Finally, in Section <ref> we analyze some select results for nucleon EDM 
induced by θ̅-angle availiable in the literature <cit.> and attempt to correct them according to our findings. 










§ CP-ODD FORM FACTORS OF SPIN-1/2 PARTICLE






 §.§ Form factors and parity mixing



In this section we argue that the ubiquitously used expression for computing -odd electric
dipole form factor F_3 on a lattice does not correspond to the electric dipole moment
measured in experiments and leads to a finite and perhaps dominant contribution from the Pauli 
form factor F_2 to the reported values of EDFF F_3 and EDM of the proton and the neutron. 
First, we recall the lattice framework for calculation of the -violating form factor F_3
first introduced in Ref. <cit.>, and later used without substantial
changes in the subsequent papers <cit.> studying the QCD θ-term, as well as more recent
developments <cit.> to study the quark chromo-EDM.



To compute nucleon form factors on a lattice, one evaluates nucleon two- and three-point functions
(see Fig. <ref>)
in presence of -violating () interactions

    C_NN̅^(p⃗,t) 
         = ∑_x⃗ e^-ip⃗·x⃗⟨ N(x⃗,t) N̅(0)⟩_ ,
    
    C_NJN̅^(p⃗^',t_sep;q⃗,t_op)
         = ∑_y⃗,z⃗ e^-ip⃗^'·y⃗+iq⃗·z⃗⟨ N(y⃗,t_sep) J^μ(z⃗,t_op) N̅(0)⟩_ .

The subscript  indicates that these correlation functions are evaluated in 
vacuum, either with a finite value of the relevant  coupling or an infinitesimal value,
i.e. performing first-order Taylor expansion of the correlation functions.
As argued in Ref. <cit.>, as well as earlier model
calculations <cit.>, the  background leads to a 
 phase in the
nucleon mass in the Dirac equation that governs the on-shell nucleon fields N, N̅ 

    (i∂ - m_N^'e^-2iγ_5)N(x) = 0 .

where the real-valued m_N^'>0 is the nucleon ground state mass in the new vacuum.
The spinor wave functions ũ_p, ũ̅_p for the new ground states

    ⟨Ω|N|p,σ⟩_
      = Z_N^'ũ_p,σ  
      = Z_N^' e^iγ_5 u_p,σ ,

also satisfy the same Dirac equation

    (p - m_N^' e^-2iγ_5) ũ_p 
      = (p - m_N^' e^-2iγ_5) e^iγ_5 u_p = 0  ,

where the chirally-rotated spinors ũ_p, ũ̅_p have a 
Lorentz-invariant  phase similarly to the mass term,

    ũ=e^iγ_5  u , 
      ũ̅ = u̅  e^iγ_5 ,

where the regular spinors u_p, u̅_p satisfy the regular Dirac equation with a real-valued
nucleon mass,

    (p - m_N^') u_p = 0 ,
      u̅_p (p - m_N^') = 0

From the above equation (<ref>) it also follows that the spinors u_p,
u̅_p transform under spatial reflection (parity P) as the regular spinors,

    γ_4 u_p=(p⃗,E) = u_p̃=(-p⃗,E) .



Below we will discuss correlation functions on a Euclidean lattice, which depend on the
Wick-rotated 4-momentum and are more conveniently expressed using the Euclidean matrices
[γ^μ]_ (<ref>). 
Whenever a confusion may arise, we will explicitly specify the type 
 (Minknowski) or  (Euclidean) of γ-matrices and 4-vectors
(see App. <ref> for details). 
The Euclidean versions of the Dirac equations (<ref>) 
for the nucleon spinors are

    (ip_ + m_N^') u_p = 0 ,
      u̅_p (ip_ + m_N^') = 0 ,

where (-ip_) = (-i) p^μ_  [γ_μ]_
  = E [γ^4]_ - ip⃗· [γ]_,
in which the Euclidean on-shell 4-momentum p^μ_=(p⃗, iE) 
is contracted with Euclidean γ-matrices and
E=√(m_N^'^2 + p⃗^'2) is the real-valued on-shell energy of the nucleon.
Due to the chiral phase (<ref>), the nucleon propagator on a
lattice (<ref>) contains chiral phases e^iγ_5.
Keeping only the ground state and omitting the exponential time dependence for simplicity,
we get

    C_NN̅(p⃗,t)|_g.s.∼∑_σũ_p,σ̃̅u_p,σ/2E_N^'
      = -ip_ + m_N^'e^2iγ_5/2E_N^'
      = e^iγ_5[-ip_ + m_N^'/2E_N^']
        e^iγ_5 .


Analogously, the expression for the nucleon-current correlator (<ref>) 
contains the phases e^iγ_5:

    C_NJN̅(p⃗^',t_sep;q⃗, t_op)|_g.s.∼∑_σ^',σũ_p^',σ^'⟨ p^',σ^' |J^μ|p,σ⟩_ũ̅_p,σ
      = e^iγ_5[∑_σ^',σ u_p^',σ^'⟨ p^',σ^' |J^μ|p,σ⟩_u̅_p,σ] e^iγ_5 .



The problem with the commonly used expression for the three-point function comes from the fact
that the physical interpretation of a parity-mixed fermion state (<ref>)
on the lattice is not clear.
In Refs. <cit.>, it is assumed that the 
nucleon matrix elements of the vector current in  vacuum have the form 
(in Minkowski space, up to sign conventions for F_3 and F_A)

    ⟨ p^',σ^' |J^μ|p,σ⟩_?= ũ̅_p^',σ^'[
        F_1(Q^2) γ^μ 
        + F̃_2(Q^2) iσ^μνq_ν/2m_N^'
        - F̃_3(Q^2) γ_5σ^μνq_ν/2m_N^'
        + F_A(Q^2) (qq^μ - γ^μ q^2)γ_5/m_N^'^2] ũ_p,σ ,

where q=p^'-p, Q^2 = -[q^2]_2=-(q^4)^2+q⃗^2, 
F_1 and F̃_2 are the Dirac and Pauli form factors, 
F̃_3 is the electric dipole form factor (EDFF),
and F_A is the anapole form factor
(notations F̃_2,3 are introduced to avoid confusion with the true F_2,3 below).
The matrix element expression (<ref>), however, disagrees with the 
literature <cit.>, 

    ⟨ p^',σ^' |J^μ|p,σ⟩_
      = u̅_p^',σ^'[
        F_1(Q^2) γ^μ 
        + F_2(Q^2) iσ^μνq_ν/2m_N^'
        - F_3(Q^2) γ_5σ^μνq_ν/2m_N^'
        + F_A(Q^2) (qq^μ - γ^μ q^2)γ_5/m_N^'^2] u_p,σ ,

in which the vertex spin matrix 

    Γ^μ(p^',p)
      = F_1 γ^μ + (F_2 + i F_3γ_5) iσ^μνq_ν/2m_N^' 
      + F_A(qq^μ - γ^μ q^2)γ_5/m_N^'^2

is contracted with the spinors 
satisfying the regular parity transformations of spinors (<ref>).
Only in this case the contribution of the form factor F_3 to the matrix element 
⟨ p^',σ^' |J^μ|p,σ⟩ transforms as an axial 4-vector
so that F_3 is indeed the -odd coupling of the nucleon to the electromagnetic
potential <cit.>.
Let us show that this is not the case if the matrix element of the current has the
form (<ref>).
Upon spatial reflection, the true 4-vectors of momenta and current have to transform as

    (p⃗^('), p^4(')) → (-p⃗^('), p^4(')) ,
       (J⃗, J^4) → (-J⃗, J^4) .

while the axial vector current A^μ transforms with the sign opposite to J^μ:

    (A⃗, A^4) → (A⃗, -A^4) .

The chirally-rotated spinors transform as

    ũ_p⃗   →ũ_-p⃗ = e^2iγ_5γ_4  ũ_p⃗ , 
    ũ̅_p⃗   →ũ̅_-p⃗ = ũ̅_p⃗ γ_4 e^2iγ_5 ,

up to an irrelevant spinor-diagonal phase factor.
Finally, remembering that the spatial momentum q⃗ is also reflected and using the identities

    γ_4σ^iν (-q⃗, q^4)_νγ_4 = -σ^iν(q⃗, q^4)_ν ,
    
    γ_4σ^4ν (-q⃗, q^4)_νγ_4 = σ^4ν(q⃗, q^4)_ν ,

we observe that a combination of the F̃_2,3 form factors transforms as

    e^2i(F̃_2 + iF̃_3) →
    e^-2i(F̃_2 - iF̃_3) .

Therefore, we conclude that the axial-vector contribution of the matrix
element (<ref>) appears because of the parity-odd form factor combination

    Im[e^2i(F̃_2 + iF̃_3)]
      = sin(2) F̃_2 + cos(2) F̃_3 .

which is different from F_3 if π n.


Since the expression (<ref>) is used in lattice calculations so
ubiquitously, we present extensive arguments that it is not correct.
The form factor F_A is irrelevant for this discussion, and will be omitted[
  It is also worth noting that F_A is not affected by the parity mixing, unlike F_2,3.].
In Appendix <ref> we directly show that it is the 
expression (<ref>) that leads to the correct -odd EDM coupling 
∼E⃗·S⃗, 
and the forward limit Q^2→0 of form factors F_2(Q^2) and F_3(Q^2) yields the anomalous
magnetic κ and electric ζ dipole moments (in units e/(2m_N)), respectively.
In Section <ref> we calculate the mass shift of a particle governed by a
Dirac equation with chirally-rotated mass in background electric field.


In this section we offer several heuristic arguments why
expression (<ref>) is not correct.
First, revisiting the form factor expression (<ref>), we note that 
the only effect of the chiral phase is to mix form factors F_2 and F_3 into each other,

    ũ̅_p^'[
        F_1 γ^μ 
        + (F̃_2 + iF̃_3γ_5) iσ^μνq_ν/2m_N^'] ũ_p
    = u̅_p^'[
        F_1 γ^μ 
        + e^2iγ_5 (F̃_2 + i F̃_3 γ_5)
        iσ^μνq_ν/2m_N^'] u_p

while the form factor F_1, as well as the omitted F_A, are independent of .
Thus, the form factors F̃_2,3 computed in Refs. <cit.>
are linear combinations of the true form factors F_2,3 

    (F_2 + i F_3 γ_5) = e^2iγ_5 (F̃_2 + i F̃_3 γ_5) ,
       or {[                         F̃_2 = cos(2) F_2 + sin(2) F_3  ,;                         F̃_3 = -sin(2) F_2 + cos(2) F_3 ,;                              ].

which is also consistent with Eq. (<ref>).

It is easy to see that the effect of the phase e^iγ_5 can be completely
removed by a field redefinition N^'= e^-iγ_5 N.
After this transformation, the on-shell nucleon field N^' will satisfy a Dirac equation 
with the real-valued mass m_N^',

    (i∂ - m_N^')N^'(x) = 0 .

A similar transformation for the nucleon correlators 
C_N[J]N̅→ C^'_N[J]N̅
  = e^-iγ_5 C_N[J]N̅ e^-iγ_5 
will remove any dependence on  altogether.
Note, however, that this is the case only if Eq. (<ref>) is used for the 
nucleon matrix elements of the current.
Thus, this phase is purely conventional and similar to the operator normalization Z_N^',
in that physical quantities cannot depend on it.
In a lattice calculation, however, this phase is not known in advance and must 
be determined numerically to be removed from the two- and three-point
correlators (<ref>,<ref>).


To make this point evident, suppose one calculated nucleon form factors in -even QCD but
using unconventional nucleon interpolating fields e^iγ_5 N with some arbitrary 
.
If Eq. (<ref>) was used, the definition of F̃_2,3  would
depend on this arbitrarily chosen .
Consequently, because of the spurious mixing (<ref>), 
the electric dipole form factor would obtain the non-zero value F̃_3 = -F_2sin(2)
in absence of any  interactions.
Analogously, the apparent nucleon magnetic moment 
μ̃_N^'=G̃_M(0)=F_1(0)+F̃_2 
would have contributions from both F_2 and F_3.
In -even QCD vacuum, F_3=0 and the mixing (<ref>) would just reduce
the contribution of F_2 to μ̃_N^' by a factor of cos(2).
This would happen because the spin operator Σ^k=1/2ϵ^ijkσ^ij
was “sandwiched” between chirally-rotated 4-spinors and

    (2S⃗)^k = Σ^k 
      = ũ̅_p^'Σ^k/2m_N^'ũ_p
      = ξ^'†σ^kξ cos(2)

where the initial and final momenta p⃗,p⃗^'≈0
and ξ,ξ^' are the corresponding 2-spinors.


The resolution to this apparent paradox is hinted by the modified form of the Gordon identity
for the spinors ũ_p, ũ̅_p.
Since these spinors satisfy the Dirac equation with the chirally rotated 
mass (<ref>), the Gordon identity takes the form

    ũ̅_p^'γ^μũ_p 
      = ũ̅_p^'[
          (p^' + p)^μ + iσ^μν(p^' - p)_ν/2 m_N^' e^2iγ_5]ũ_p ,

which is obtained from the standard Gordon identity by replacing 
m_N^'→ m_N^' e^2iγ_5.
The form of the nucleon-current vertex must be compatible with the form of the Gordon identity,
which, among other things, relates form factors F_1,2 to G_M.
Therefore, to make the nucleon-current vertex compatible with the spinors
ũ_p, ũ̅_p, 
the nucleon mass in the F̃_2,3 terms in Eq. (<ref>) 
has to be adjusted similarly to Eq. (<ref>),
which leads back to the correct expression (<ref>).


Finally, we emphasize that Eqs. (<ref>) and
(<ref>) result in different prescriptions for analyzing the three-point
nucleon-current correlators:

    C_NJN̅(p⃗,t_sep;q⃗,t_op)|_g.s.?= 
        e^-E_N^'^'(t_sep-t_op)-E_N^'t_op
        e^iγ_5-ip^'_ + m_N^'/2E^'_N^'{e^iγ_5}^?
        Γ^μ_{e^iγ_5}^?
        -ip_ + m_N^'/2E_N^'
        e^iγ_5

where phase factors in curly braces {e^iγ_5}^? are present only if 
one uses the (incorrect) Eq.(<ref>).
In the above equation, we have introduced the Euclidean nucleon-current vertex

    Γ^μ_(p^', p) =[F_1γ^μ + (F_2 + iγ_5 F_3)
        σ^μνq_ν/2m_N^']_









 §.§ EDM energy shift from Dirac equation


We argued in the previous section that one has to use regular ”even” spinors satisfying
Eq. (<ref>) to evaluate the nucleon matrix elements even if the QCD vacuum
is -broken, contrary to the previous works <cit.>.
Most of the ambiguity must have resulted from the notion that in a -broken vacuum, particles
are no longer parity eigenstates, hence, the argument goes, the nucleon must be described by a
parity-mixed spinor.
This argument is rather confusing because parity transformations of fermion fields are fixed only
up to a phase factor, and only a fermion-antifermion pair may have definite parity.
To clarify this question, in this section we calculate the energy spectrum of a particle described
by the Dirac operator _N with the complex mass me^-2iγ_5 and with
magnetic and electric dipole interactions in the form (<ref>) in the
background of uniform magnetic and electric fields.
Such an operator is exactly the nucleon effective operator on a lattice.
The zero modes of this operator (i.e., the poles of its Green's function) must correspond to 
particle eigenstates, and their calculation avoids the spinor phase ambiguity completely.
The energy shifts linear in these fields are then identified with the magnetic κ and 
electric ζ dipole moments, respectively.



The effective action for the Euclidean lattice nucleon field in the  vacuum and
point-like electromagnetic interaction introduced via “long derivative” is

    _int = N̅ [i∂ 
          - Q γ^μ A_μ - m e^-2iγ_5] N ,

where we neglect the momentum transfer dependence of the nucleon form factors
for simplicity, setting F_1 to a “point-like” value Q=F_1(0)=const.
In the absence of electromagnetic potential A_μ, the nucleon
propagator ⟨ NN̅⟩ takes the form (<ref>).
We add effective point-like anomalous magnetic κ̃=F̃_2(0)
and electric ζ̃=F̃_3(0) dipole interactions 
to the interaction vertex

    Q γ^μ→ 
      Q γ^μ + κ̃iσ^μνq_ν/2m 
      - ζ̃γ_5σ^μνq_ν/2m

Using conventions (<ref>-<ref>) as well as
(<ref>,<ref>), the Dirac equation for N becomes

    [p - Qγ^μ A_μ 
      - (κ̃+ iζ̃γ_5) 1/2 F_μνσ^μν/2m
          - m e^-2iγ_5] N
      = 0

We are going to find the energy levels of the particle in presence of constant field strength 
F_μν. 
To avoid irrelevant complications, we consider only a neutral particle with Q=0. 
Using the identity (<ref>) to trade γ_5 for
F_μν→F̃_μν, we obtain the Dirac operator in the block-diagonal form 
in the chiral basis (<ref>):

    p - 1/2(κ̃F_μν 
          - ζ̃F̃_μν)σ^μν/2m 
      - m e^-2iγ_5
      =([      -M E-p⃗·σ⃗; E+p⃗·σ⃗    -M^† ]) ,

where 
M = me^2i - 1/2m(κ̃- i ζ̃)(H⃗ + iE⃗)·σ⃗.
In the rest frame, p⃗=(0⃗,E_0), and the operator (<ref>) has 
solutions if

    det([   -M  E_0;  E_0 -M^† ]) = 0
      ⇔  det(E_0^2 - M^† M) = det(E_0^2 - M M^†) = 0

Up to terms linear in κ̃, ζ̃, the normal operator M^† M is

    M^† M    = m^2 - 1/2[ e^ 2i(κ̃+ iζ̃)(H⃗ - iE⃗) 
                  + e^-2i(κ̃- iζ̃)(H⃗ + iE⃗)] ·σ⃗+ O(κ̃^2,ζ̃^2)
    
       = m^2 - [κH⃗ + ζE⃗]·σ⃗+ O(κ^2,ζ^2)

where in the last line we redefined 
e^2i(κ̃+ iζ̃)=(κ + i ζ), 
which is the same transformation as Eq. (<ref>).
Finally, the energy of the particle's interaction with the E&M background is

    E_0 - m = -κ/2mH⃗·Σ̂- ζ/2mE⃗·Σ̂+ O(κ^2,ζ^2)

where Σ̂ is the unit vector of the particle's spin.
From the interaction energy, we conclude that indeed 

    κ=F_2(0) ,   ζ=F_3(0)
 are the
particle's magnetic and electric dipole moments. 
For a neutral particle such as the neutron, the form factor F_2(0) is indeed the full magnetic 
moment.

Thus, we have shown that if particle's field is governed by the Dirac equation with a complex
mass (<ref>), electric and magnetic dipole moments have to be properly
adjusted (κ̃, ζ̃)→(κ, ζ).
This adjustment follows from redefining the field and the operator

    N → N^' = e^-iγ_5 N ,
      D̃_N→D_N
      = e^iγ_5D̃_N e^iγ_5

to remove the complex (chiral) phase from the mass,
where D̃_N(D_N) contains κ̃(κ) and ζ̃(ζ).




 §.§ EDM energy shift in Euclidean space

In order to verify our findings, in this paper we calculate the EDM of the neutron on a lattice
using two methods: from the energy shift in the background electric field method and using the
new formula for the -odd form factor F_3. 
The electric field is introduced following Ref. <cit.> and preserving the 
(anti)periodic boundary condition in time.
Such electric field <cit.> is analytically continued to an imaginary value.
If the particle's electric dipole moment is finite and real-valued, the energy shift will be
imaginary, which might present a problem in the analysis of corresponding lattice correlators.
However, in our methodology, the -odd interaction is always infinitesimal, and so are 
the electric dipole moments and the corresponding energy shifts, which are extracted from the 
first-order Taylor expansion of the nucleon correlation functions in the -odd interaction.
In this paper, we study only neutral particles, because analysis of charged particle propagators
is more complicated <cit.>.

In this section, we repeat the calculation of Sec. <ref> for a neutral
particle on a Euclidean lattice, which has on-shell Euclidean 4-momentum 
p_=(p⃗,iE), with energy E=√(E_0^2 + p⃗^2) up to discretization errors.
The energy at rest E_0 is modified from the mass m due to electric and magnetic dipole
interactions.
To avoid any confusion, we imply no relation between the Minkowski E⃗,H⃗ and Euclid
, electric and magnetic fields. 
Instead, we introduce ad hoc uniform Abelian fields on a lattice 
(see Fig. <ref>) preserving boundary conditions in both space and
time <cit.> that probe the MDM and EDM: 
magnetic ϵ^ijk^k = (∂_i _j - ∂_j _i) = n_ijΦ_ij
(no summation over i,j)

    {[                  _x,j      = n_ij Φ_ij  x_i;       _x,i|_x_i=L_i-1 = -n_ij Φ_ij  L_i x_j ].

and electric ^k = (∂_k _4 - ∂_4 _k) = n_k4Φ_k4 

    {[                  _x,4      = n_k4 Φ_k4  x_k;       _x,k|_x_k=L_k-1 = -n_k4 Φ_k4  L_k x_4 ]. ,

where Φ_μν = 6π/L_μ L_ν is the quantum of flux through a plaquette
(μν) and n_μν is the corresponding number of quanta.
The fractional quark charges Q_u=2/3, Q_d=-1/3 and periodic boundary conditions require 
that the flux through the edge of the lattice is L_μ L_ν·Φ_μν = 3·2π.
From potentials (<ref>,<ref>), the Euclidean field strength 
tensor _μν = ∂_μ_ν - ∂_ν_μ

    _μν 
    = ([       1   2   3   4;   1   0  ^3 -^2  ^1;   2 -^3   0  ^1  ^2;   3  ^2 -^1   0  ^3;   4 -^1 -^2 -^3   0 ])

with
=(n_23Φ_23, n_31Φ_31, n_12Φ_12) and
=(n_14Φ_14, n_24Φ_24, n_34Φ_34).




We start from the effective EDM and MDM coupling in the nucleon-current vertex.
The Dirac operator for the nucleon on a lattice is 
D+m=γ^μ(∂_μ + i Q _μ) + m, which we extend to include the
point-like effective interactions from Eq. (<ref>)

    [ip +m +iQ_q 
      -(1/2σ^μν_μν)  κ +iζγ_5/2m]_

with κ=F_2(0) and ζ=F_3(0).
We use Euclidean matrices γ^μ (<ref>) 
and [γ_5]_=(γ^1γ^2γ^3γ^4)_ (<ref>)
[
  The results are manifestly independent from the basis of γ-matrices used, 
  if the relation between γ_5 and γ^μ is kept unchanged.]
and the plain-wave fields ψ_p(x) and _q,μ(x) depending on 
the Euclidean 4-momenta p, q as

    ψ_p(x)∼ e^i p x ,
      ∂_μψ_p(x) ↔ i p_μ u_p ;
      _q,μ(x) ∼ e^i(p^'-p) x = e^i q x , 
      ∂_ν_q,μ(x) ↔ i q_ν_μ .

The mass m in the above equation (<ref>) is chosen real and positive, 
since any chiral phase factor may be removed with a field redefinition (<ref>), 
which at the same time rotates the dipole couplings (κ, ζ) into the physical magnetic 
and electric dipole moments, as has been shown in Sec. <ref>.
After setting the charge Q=0 and the momentum p⃗=0, we use

    σ_^ij 
      = ϵ^ijk([ -σ^k     ;      -σ^k ]) ,
      σ_^i4 
      = ([  σ^i     ;      -σ^i ]) ,

and transform the operator (<ref>) into the block-diagonal form,
and find the condition for on-shell fermion energies

    ([   _- -E_0; -E_0   _+ ])
    =0 
      ⇔  (E_0^2 - _+ _-) = (E_0^2 - _- _+) = 0

where

    _± = m + 1/2m(κ∓ i ζ)  (±)·σ⃗ .

The on-shell energies are then determined by the eigenvalues of the spin-dependent operator

    _∓_±
      =m^2 +κσ⃗·- ζσ⃗· i +O(k^2, ζ^2) ,
    
    
    E_0-m=κ/2m ·Σ̂- ζ/2m i·Σ̂+O(k^2, ζ^2) ,

where Σ̂ is the direction of the particle's spin.
Note that the electric field enters Eq. (<ref>) as i,
with an imaginary factor emphasizing that its value has been analytically continued 
to the imaginary axis, and the corresponding energy shift is purely imaginary.
Equation (<ref>) provides a prescription for extracting the EDM and MDM from
energy shifts of a neutral particle on a lattice in uniform background fields.







§ CEDM-INDUCED EDM AND EDFF ON A LATTICE


In our initial calculation of cEDM-induced nucleon EDMs,
we use two lattice ensembles with Iwasaki gauge action and N_f=2+1 dynamical domain wall
fermions: 16^3×32 with m_π≈420 MeV <cit.> 
and 24^3×64 with m_π≈340 MeV <cit.>.
The ensemble parameters are summarized in Tab. <ref>.
We use identical ensembles, statistics, and spatial sampling per gauge configuration 
in both calculation methods discussed in further sections.




We use all-mode-averaging <cit.> framework to optimize sampling, 
in which we  approximate quark propagators with truncated-CG solutions 
to a Möbius operator <cit.>.
We use the Möbius operator with short 5th dimension L_5s and complex s-dependent
coefficients b_s + c_s = ω_s^-1 (later referred to as “zMobius”) that approximates 
the same 4d effective operator as the Shamir operator with the full L_5f=32 (DSDR) 
or L_5f=16 (Iwasaki).
The approximation is based on the domain wall-overlap equivalence

    [^DWF]_4d = 1+m_q/2 - 1-m_q/2γ_5
    ϵ_L_5(H_T) ,
      
    H_T = γ_5_W/2+_W ,
    ϵ^Möbius_L_5s(x) 
      = ∏_s^L_5s(1+ω_s^-1x) - ∏_s^L_5s(1-ω_s^-1x)/∏_s^L_5s(1+ω_s^-1x) + ∏_s^L_5s(1-ω_s^-1x)≈ϵ^Shamir_L_5f(x) 
       .

where the coefficients ω_s are chosen so that the function 
ϵ^Möbius_L_5s(x) is the minmax approximation to the 
ϵ^Shamir_L_5f(x).
We find that L_5s=10 is enough for efficient 4d operator approximation.
Shortened 5th dimension reduces the CPU and memory requirements:
for example, L_5f=16 is reduced to L_5s=10 saving 38% of the cost.
We deflate the low-lying eigenmodes of the internal even-odd preconditioned operator,
to make the truncated-CG AMA efficient. 
The numbers of deflation eigenvectors N_ev and truncated CG iterations N_CG are given in
Tab. <ref>.
We compute 32 sloppy samples per configuration. 
To correct any potential bias due to the approximate  operator and the truncated CG
inversion, in addition we compute one exact sample per configuration using the Shamir operator.
The latter is computed iteratively by refining the solution of the “zMobius” 
to approach the solution of the Shamir operator, again taking advantage of the short L_5s
and deflation.





 §.§ Parity-even and -odd nucleon correlators
















The EDFF F_3 is a parity-odd quantity induced by  interactions.
To compute the effect of -odd interactions, we modify the lattice action 

    S→ S + i δ^ S = S + i∑_i,x c_i [^_i]_x

where c_i are the -odd couplings such as the QCD θ-angle, quark (chromo-)EDMs, etc.
We Taylor-expand QCD+ vacuum averages in the couplings c_i.
For example, for the three-point function, we get[
  In this section, all conventions for correlators, form factors, and momenta are Euclidean.]

    ⟨ N [q̅γ^μ q] N̅⟩_
      = 1/Z∫  U ψ̅ψ e^-S - iδ^ S  N [q̅γ^μ q] N̅ 
      = C_NJN̅ - i∑_i c_i  δ^_i C_NJN̅
        + O(c_ψ^2) ,

where C_… and δ^ C_… stand for the -even and -odd correlators.
Similarly, we also analyze the effect of  interaction on the nucleon-only correlators.
In total, we calculate the following two- and three-point -even correlators as well as 
three- and four-point -odd correlators,

    C_NN̅   = ⟨ N N̅⟩ ,
         δ^_i C_NN̅   = ⟨ N N̅ ·∑_x [^_i]_x ⟩ ,
    
    
    C_NJN̅   = ⟨ N [q̅γ^μ q] N̅⟩ ,
         δ_i^ C_NJN̅   = ⟨ N [q̅γ^μ q] N̅·∑_x [^_i]_x ⟩ ,

where ⟨·⟩ stand for vacuum averages computed with -even QCD action S.
In Sec. <ref>, we also modify the action S to include the uniform background electric
field as the probe of the electric dipole moment.
In this work, we study only the quark chromo-EDM as the source of  violation,

    _ψ G^  
      = 1/2ψ̅[G_μν]^clovσ^μνγ_5ψ
      = 1/2ψ̅(g_S G^a,cont_μν T^a)σ^μνγ_5ψ ,

where G^a,cont is the continuum color field strength tensor and 
the “clover” [G_μν]^clov gauge field strength tensor 
on a lattice is (see Fig. <ref>)

    [G_μν]^clov   = 1/8i[
      ( U^P_x,+μ̂,+ν̂ + U^P_x,+ν̂,-μ̂ 
      + U^P_x,-μ̂,-ν̂ + U^P_x,-ν̂,+μ̂) 
      - h.c.]

Insertions of the quark-bilinear cEDM density (<ref>) can generate both connected 
and disconnected contractions, similarly to the quark current.
In this work, we calculate only the fully connected contributions to these correlation functions 
shown in Fig. <ref>. 
The disconnected contributions (see Fig. <ref>) are typically much more
challenging to calculate, and we will address them in future work. 
Neglecting the disconnected diagrams will not affect the comparison of the form factor and the
energy shift methods, because they are omitted in both calculations.










To compute the connected diagrams, we insert the quark-bilinear cEDM density (<ref>)
once in every ψ-quark line of C_N JN̅ diagrams, generating the four-point functions 
shown in Fig.<ref>.
We evaluate all the connected three- and four-point contractions using the forward 
and the set of sequential propagators shown in Fig. <ref>.
In addition to the usual one forward  and two backward (sink-sequential)  propagators,
we compute one cEDM-sequential  and four doubly-sequential ({cEDM, sink}-sequential)
(+) propagators per sample.
For every additional value of the source-sink separation t_sep and sink momentum p⃗^', additional backward  and doubly-sequential (+) propagators must be
computed, i.e.

    N_ = N_ = 1 , 
       N_=N_q N_sep N_mom ,
       N_+ = N_q N_ψ N_sep N_mom

where N_q is the number of separate flavors in the quark current and N_ψ is the number of
separate flavors in the  operator.
The connected -even two- and three-point correlators do not require any additional inversions.
In this scheme, we perform only the minimal number of inversions required for computing 
all the diagrams for the neutron and proton EDM induced by 
a connected flavor-dependent quark-bilinear  interaction 
with the two degenerate flavors u and d.
Compared to Ref. <cit.>, in which a finite small O(ϵ)
-odd perturbation term is added to the quark action that results in modified quark
propagators

    _m^-1→(_m + iϵσ^μνG̃_μν)^-1 ,

our four-point contractions correspond to directly computing the first derivative 
(∂ C_2,3^ϵ / ∂ϵ)_ϵ=0, thus avoiding 
any higher-order dependence on ϵ and obviating the ϵ-extrapolation.
As a cross-check, we have verified our contraction code on a small test lattice by replacing 
propagators _m^-1η with

    _m^-1η→ 
      [_m^-1 - _m^-1 (iϵΓ) _m^-1]η

to approximate [_m + iϵΓ]^-1η,  
where Γ=1/2 G_μνσ^μνγ_5.
Using these “-perturbed” propagators, each of which needed two inversions, 
we have computed the nucleon C^ϵ_NN̅ and nucleon-current 
C^ϵ_NJN̅ correlators, and compared their finite-difference 
ϵ-derivatives to δ C_NN̅ and δ C_NJN̅.

We use only one value of the sink momentum p⃗^'=0.
We compute nucleon-current three- and four-point correlators with two source-sink separation 
values t_sep={8,10}a={0.91,1.15} fm for the 16^3×32 ensemble, 
and three t_sep={8,10,12}a={0.88,1.11,1.33} fm for the 24^3×64 ensemble.
For the 24^3×64 lattice, we use the Gaussian-smeared sources with APE-smeared gauge links
using parameters optimized for overlap with the ground state <cit.>,
while for the 16^3×32 ensemble we used the smearing parameters from
Ref. <cit.>.
The effective nucleon mass plots for the ensembles are shown in Fig. <ref>.
Correlators C_NJN̅ and δ^C_NJN̅ are computed with
the polarization projector

    T^+_S_z+=1+γ_4/2(1+Σ_3)=1/2(1+γ_4)(1-iγ_1γ_2) ,
,
while correlators C_NN̅ and δ^ C_NN̅ are computed with all 16
polarizations and saved to be used later for disconnected contractions.
We reduce the cost of computing backward propagators with the widely-used “coherent”
trick combining 2 backward sources from samples separated by L_t/2 into one inversion.
Combining 4 samples resulted in a large increase in the statistical uncertainty negating the
cost-saving advantages.







 §.§ Nucleon form factors
 

Following the discussion in Sec.<ref>, we use the form factor decomposition that
is different from the previous works <cit.>,

    ⟨ N_p^' | q̅γ^μ q | N_p⟩
      = u̅_p^'[ F_1(Q^2) γ^μ
          + F_2(Q^2) σ^μν q_ν/2m_N
          + F_3(Q^2) iγ_5σ^μν q_ν/2m_N] u_p .

where the spinors u_p,u̅_p^' have positive parity.
Details of evaluating kinematic coefficients for form factors F_1,2,3
are given in Appendix <ref>.
We use the standard plateau method to evaluate both -even and -odd matrix elements of
the nucleon 

    [δ^] _NJN̅ (t_sep,t_op)
      = [δ^] C_NJN̅(t_sep,t_op)/c_2^'(t_sep)√(c_2^'(t_sep)/c_2(t_sep)c_2^'(t_op)/c_2(t_op)c_2(t_sep-t_op)/c_2^'(t_sep-t_op))

where the two-point functions are projected with the positive-parity polarization matrix
T^+=1/2(1+γ_4),

    c^(')_2(t) = [ T^+· C_NN̅(p⃗^('), t) ]  .

The three central points on the ratio plateaus are taken as the estimate of the ground state 
matrix elements.
This is a crude estimate and improved analysis of excited states is necessary for better
control of systematic uncertainties.
However, we find that our results change insignificantly with increasing source-sink separation 
(see Figs. <ref>, <ref>), therefore we conclude
that excited state effects cannot influence the main conclusions of the paper.


We calculate the Dirac and Pauli form factors F_1,2 using a correlated χ^2 fit to the
matrix elements of the quark vector current (“overdetermined analysis”). 
The system of equations for form factors is reduced by combining equivalent equations to
reduce the system dimension and make estimation of the covariance matrix more stable 
(see, e.g., Ref. <cit.> for details).
The quark current operator is renormalized using renormalization constants 
Z_V=0.71408 for 24^3×64 <cit.> and 
the chiral-limit value Z_V=Z_A=0.7162 for 16^3×32 <cit.> ensembles.
We show the momentum dependence of the resulting Sachs electric and magnetic form factors 

    G_E(Q^2) = F_1(Q^2) - Q^2/4m_N^2 F_2(Q^2) ,
      
    G_M(Q^2) = F_1(Q^2) + F_2(Q^2) ,

for the proton and the neutron (connected-only) for both ensembles in Fig. <ref>.
Our data for form factors G_E,M show no significant systematic variation with 
increasing the source-sink separation t_sep.








In order to compute the form factor F_3, we first need to calculate the parity mixing angle
 in order to subtract the F_1,2 mixing terms.
Expanding the nucleon two-point function C_NN̅^(t) to the first order in 
∼ c_ψ G and assuming that the ground state dominates for sufficiently large t,

    C_NN̅(t) - i c_ψ G δ^ C_NN̅(t) +O(c_ψ G^2) 
      t→∞= |Z_N|^2[1+γ^4/2 + iγ_5 + O(^2)] e^-m_N t

we use the projectors T^+=1+γ_4/2 and T^+γ_5 to calculate 
the “effective” mixing angle (t) normalized to c_ψ G=1

    ^eff(t) 
      = -[ T^+γ_5 ·δ^ C_NN̅ (t)]/[ T^+ · C_NN̅(t)]t→∞ = /c_ψ G

The time dependence of the ratios (<ref>) for both ensembles is shown in 
Fig. <ref>.
The quark flavors in the cEDM interaction are shown respective to the proton, and 
for the neutron must be switched u↔ d due to the isospin symmetry.
The plateau is reached for time t≥8, and we extract the  values from a constant fit 
(weighted average) to points t=8…11.
An interesting observation is that the mixing angle depends very strongly on the flavor 
involved in the  interaction.
Thus, for the proton P_δ = u_δ(u^T C γ_5 d), in which the d-quark enters
together with u as a scalar diquark, the d-cEDM does not lead to any parity mixing.


Finally, the electric dipole form factor F_3 is calculated from the -odd 
four-point correlator δ^ C_NJN̅.
Similarly to the extraction of  above, we can expand the  
three-point function in the -odd interaction.
We extract the matrix elements using the ratios (<ref>) of 
polarization-projected three-point functions [T·^_NJN̅]
to -even two-point functions (<ref>).
Expanding the ratio in ∼ c_ψ G, we get

    [T(_NJN̅
        - i c_ψ G δ^_NJN̅ +O(c_ψ G^2))]
    
    t→∞= 
      ∑_i=1,2[ _ i^(T) + i_ i^({T,γ_5})]F_i
        + _ 3^(T) F_3 + O(^2)

where ^(T)_ 1,2,3 are the kinematic 
coefficients (<ref>-<ref>) for form factors F_1,2,3 
computed with the polarization matrix T and with →_ (<ref>).
Matching the O(c_ψ G^1) terms in the above expansion and neglecting excited states, 
we obtain

    i _ 3^(T)F̂_3 
      = [T·δ^_NJN̅]_g.s. 
        + ∑_i=1,2_ i^({T,γ_5}) F_i

The second term in the RHS of the above equation is the mixing subtraction. 
Its form indicates that the mixing between form factors F_1,2 and F_3 happens only
because of the mixing of the polarization of the nucleon interpolating fields on a lattice.
This is substantially different from expressions used in Refs. <cit.>,
which also include additional subtraction term (-2 F_3) because of spurious mixing of
F_2 and F_3 in the vector current vertex (<ref>).

Although both timelike and spacelike components of the current can be used to calculate F_3, 
in practice we find that the time component J^4 yields much better precision than 
the spacelike component J^3. 
Due to the larger uncertainty of the J^3 signal, combining both components did not result in
improved precision of the  F_3 form factor.
If only the J^4 component is used, the overdetermined fit to matrix elements is not required,
and for T=T^+_S_z+ = 1+γ_4/2(1-iγ_1γ_2) 
from Eqs. (<ref>,<ref>)

    (1+τ) F_3(Q^2)
      = m_N/q_3 _[T^+_S_z+·δ^_N J^4N̅]
        -  G_E(Q^2)

where τ is the kinematic variable (<ref>).
It is remarkable that for the neutron the subtraction term ∼ G_E is zero in the
forward limit.
In fact, if one uses the traditional formula for extracting the neutron EDM 
d_N = F_3(0)/(2m_N), a large contribution (-2 F_2(0))/(2m_N) comes from 
the spurious mixing if  is not zero.
In Section <ref> we will discuss the currently available lattice results for the 
neutron and proton EDM induced by the QCD θ-term.

To compute form factors from data with each source-sink separation t_sep, 
we use the value = ^eff(t_sep) in Eq. (<ref>)
to subtract the mixing.
The results for the EDFF F_3 are shown in Fig. <ref>.
Despite relatively high statistics, the signal for the cEDM-induced form factor is noisy.
There is no significant dependence on the source-sink separation t_sep.
Since the cEDM operator is not renormalized, it can include contributions from other operators of
dimension 5, as well as operators from lower dimension 3 <cit.>.
One peculiar feature of these results is that, similarly to , the contribution to the
proton EDM comes mostly from the u-cEDM, while the contribution to the neutron comes mostly from
the d-cEDM.
However, a substantial increase in statistics, as well as more elaborate analysis of excited
states, are required to confirm these observations.

The electric dipole moment is determined by the value of the form factor F_3(Q^2) at zero.
This value is not directly calculable, and one has to extrapolate the Q^2>0 data points to
Q^2=0. 
In Figure <ref> we show linear extrapolation of these form factors
using the three smallest Q^2>0 points.
Other fit models are not warranted until the statistical precision is substantially improved.











 §.§ Neutron electric dipole moments from energy shifts


Calculation of the dipole moment using uniform background field has an advantage that no form
factor extrapolation in momentum is required, because the energy shift depends on the forward
matrix element of the nucleon.
This calculation is easier for the neutron than for the proton, because in case of a charged 
particle there are additional complications due to its constant acceleration,
which makes analysis of its correlation functions more complicated <cit.>.

On the other hand, since the uniform background field is quantized on a lattice, these fields
cannot be made arbitrarily small.
In fact, the field quanta are very large and their magnitudes are comparable to the QCD scale, 
especially on the smaller 16^3×32 lattice.
Because of the fractional charges of quarks, there is additional factor of 3 in the minimal value 
of the electric field, which is quantized in multiples of _0 = 6π/a^2 L_x L_t. 
The _0 values are shown in Tab. <ref>, and for the smaller 16^3×32 
lattice the minimal electric field is _0=0.110 GeV^2 = 560 MV/fm.
Such electric field pulls the u quark in the neutron with tension 
≈(270  MeV)^2, or approximately 40% of the QCD string tension,
and may deform the neutron too far away from the ground state.

We introduce the uniform electric field on a lattice as described in
Sec. <ref> along the z direction.
Using modified QCD+U(1) gauge links, we calculate the regular nucleon correlator 
C_NN̅,, as well as the correlator with the insertion of 
-odd interaction, in full analogy with Sec. <ref>, e.g.

    δ^ C_NN̅,(p⃗, t) 
      = ∑_y⃗ e^-ip⃗(y⃗ - x⃗)⟨ N(y⃗, t) 
        N̅(x⃗, 0)·^_ψ G⟩_

The modified gauge links are used in both computing the propagators and construction of the
smeared sources and sinks.
In fact, since the individual quarks are charged, smearing their distributions with only 
the QCD gauge links breaks the covariance and makes the calculation dependent on the choice 
of the gauge of the electromagnetic potential.
The QCD links used in Gaussian smearing are first APE-smeared, and then the electromagnetic potential 
is applied to them.

From Eq. (<ref>), the energy of a particle on a lattice with the spin polarized 
along the electric field =ẑ is shifted by the imaginary value 
δ E = -(ζ / 2m) i.
The nucleon correlator at rest (p⃗ = 0) thus must take the form

    C^_NN̅,(p⃗=0, t) 
      = |Z_N|^2 e^iγ_51+γ_4/2[
          1+Σ_z/2 e^-(m+δ E)t
         +1-Σ_z/2 e^-(m-δ E)t]e^iγ_5

As with the -odd form factor F_3, expanding the correlator up to the first order in
c_ψ G∼∼δ E∼ζ, we get

    C_NN̅, - i c_ψ G δ^ C_NN̅,t→∞= |Z_N|^2 e^-m_N t [1 + γ_4/2 + iγ_5 
        - Σ_z δ E  t ]

for the electric dipole moment we obtain the following estimator for the 
effective energy shift:

    ζ^eff(t) = 2 m_N d_N^eff(t) = -2m_N/_z[R_z(t+1) - R_z(t)] , 
       R_z(t) = Tr[ T^+  Σ_z δ^ C_NN̅,_z(t)]/Tr[ T^+   C_NN̅,_z(t)] .

We have computed the neutron correlation functions with two values of the electric field 
=_0 and 2_0.
The results for both ensembles are shown in Fig. <ref>. 
We choose t=6…9 as the common plateau to estimate the value of ζ on both ensembles and 
both flavors in the cEDM operator.
In the case of d-cEDM, we observe non-zero values of the energy shift.
Also the EDM values computed with =_0 and 2_0 agree well with each other,
indicating that the energy shift is linear in  and our EDM result does not depend on the
polarizing effect of the electric field.














 §.§ Numerical comparison of the form factor and energy shift methods


The normalization and the sign convention of the dimensionless EDM ζ in Sec. <ref>
are identical to those of F_3(0) in Sec. <ref>, and we plot
them for comparison in Fig. <ref>.
We observe satisfactory agreement between the values of ζ computed in the uniform background
method and the values obtained from the Q^2→0 extrapolation of form factors F_3n(Q^2).

In order to check how the spurious mixing affects the results, 
in Fig. <ref> we also plot the values of form factors 
computed with the old formula used in Refs. <cit.>

    F̃_3 = F_3 - 2 F_2 .

This formula obviously gives a value for F̃_3 different from F_3 only if 
 is large.
In the case of u-cEDM, the value  for the neutron is small, and there is no observable
difference between F_3 and F̃_3.
However, in the case of d-cEDM, the difference is remarkable.
Neither of the three sources of uncertainty: excited state bias in the energy shift calculation, 
excited state bias in the form factor calculation, nor the Q^2→0 extrapolation of the form
factors, can plausibly change the outcome of this comparison, due to the large value of .
The agreement between the new form factor extraction formula and the energy shift method is one of
the main results of this paper, and serves as a numerical cross-check of the analytic derivation.





We collect the values of , extrapolated F_3(0), and ζ_n from the background field
method in Tab. <ref>.









§ CORRECTIONS TO EXISTING Θ-INDUCED NEDM LATTICE RESULTS





In Section <ref> it has been shown that the commonly used formula for extracting the
form factor F_3 from  nucleon matrix elements on a lattice is incorrect.
This formula has been used in all of the papers that compute QCD θ-induced nucleon 
EDM <cit.>.
Fortunately, the correction has a very simple form (<ref>), in which
F̃_2,3 refer to the old results and F_2,3 refer to corrected results.
Unfortunately, Refs. <cit.> offer a broad spectrum of conventions 
for F̃_3 and  differing in sign and scale factors. 
However, by comparing expressions for polarized -odd 
matrix elements of the timelike component of the vector current J_4
we can deduce the appropriate correction 
using that reference's conventions.
For example, using Eq.(55) from Ref. <cit.>,

    Π^0_3pt,Q(Γ_k=i/4(1+γ_0)γ_5γ_k)
         ∼i Q_k/2m_N[α^1 (F_1 + E_N + 3m_N/2m_N F_2) 
        + E_N + m_N/2m_NF̃_3]
       = i Q_k/2m_N[ α^1 G_E 
        + (1 + τ)
    (F̃_3 + 2α^1 F_2)
    _F_3] ,

where τ=E_N-m_N/2m_N introduced in Eq.(<ref>) and G_E=F_1-τ F_2 is the
Sachs electric form factor.
Comparing the above equation to the expected form (<ref>), 
for the corrected value of F_3  we obtain[
  Note that this correction is the opposite compared to Eq. (<ref>), which results
  from a difference in used conventions.]

    F_3(Q^2) = F̃_3(Q^2) + 2α^1 F_2(Q^2)  ,

which should hold for any value of Q^2.





Although it is more suitable that the original authors of Refs. <cit.> 
reanalyze their data with these new formulas,
it is interesting to examine whether the presently available lattice calculations 
necessarily yield non-zero values θ̅-induced nucleon EDM 
after corrections similar to Eq. (<ref>) have been applied.
The most precise result for F_3n(0) that also allows us to perform the correction
unambiguously is Ref.<cit.>, 
which reports an 8σ non-zero value for F_3(0)=-0.56(7) from calculations with dynamical 
twisted-mass fermions at m_π=373 MeV. 
However, when we apply the corresponding correction (<ref>), 
the value becomes 0.09(7) and essentially compatible with zero.


Calculations with finite imaginary θ-angle <cit.> yield the most precise 
values of the neutron EDM to date. 
However, they do not contain sufficient details to deduce the proper correction for F_3.
It must also be noted that it is not clear if the sign of the -odd interaction ∼G̃G 
is consistent in all of the Refs. <cit.>.
On the other hand, all the reported non-zero results for the proton and neutron EDM agree in sign 
with F_3n(0) < 0 and F_3p(0) > 0, and it is reasonable to assume that any 
differences in the conventions are compensated in each final reported EDM value.
Furthermore, because the θ-angle is equivalent to a chiral rotation of quark fields, 
it is then reasonable to assume that upon conversion to some common set of conventions,
e.g., those of Ref. <cit.>, the sign of the chiral rotation angle α 
agrees between different calculations.
Based on these plausible assumptions, we deduce that the results in <cit.> 
must be corrected as F^θ_3 = F̃^θ_3 + 2α(θ) F_2[
  Strictly speaking, for finite values of θ̅ and α̅(θ̅), one has to 
  use the hyperbolic “rotation” formula cosh(2α)F_3 = F̃_3 + sinh(2α) F_2, 
  in which we neglect O(α^2) terms because |α|≲0.15, while the precision is 
  only ≈10%.],
where α<0, in analogy with Ref. <cit.>.
The data for α̅^θ and F̃_3^θ(0) at finite θ̅ values are 
extracted from figures in Ref. <cit.>.
The original F̃_3^θ(0) and the corrected F_3^θ(0) values are shown in
Fig. <ref>.
Following Ref. <cit.>, the corrected F_3^θ(0) values are interpolated to 
θ̅→0 using a linear+cubic fit F_3(0)θ̅+ Cθ̅^3
and the resulting normalized values F_3(0)=dF_3^θ/dθ̅|_θ̅=0 
are given in Tab. <ref>.
We observe that the corrected values at both the finite and zero θ̅ values agree with
zero at ≲2σ level.

Corrections to other results <cit.>,

 may be done on the similar basis[
  Correction to results in Ref. <cit.> require the corresponding values for F_2,
  which we could not locate in published works.].
The resulting values are also collected in Tab. <ref>, and in most cases
they are compatible with zero, deviating at most 2σ. 
We emphasize that, apart from Ref. <cit.>, these corrections are made using the
sign assumptions discussed above.
If our assumptions are wrong, the corrected central values will be approximately twice as large
compared to the originally reported values.
Although we find our assumptions plausible, and thus the corrected values in
Tab. <ref> most likely valid, it is up to the authors of
Refs. <cit.>
to reanalyze their data and confirm or deny our findings.
It is possible that the difference between the lattice values of the neutron EDM and
phenomenological estimates d_n∼ O(10^-3…10^-2) θ̅ e·fm <cit.>, 
which has been ascribed to chiral symmetry breaking of lattice fermions and the 
heavy quark masses used in simulations, can disappear when the proper corrections are applied.





§ SUMMARY AND CONCLUSIONS


Among our most important findings in this paper is the new formula for analysis of nucleon-current
correlators computed in  vacuum and extraction of the electric dipole form factor F_3.
We have demonstrated, both analytically and numerically, that the analysis of the 
θ̅-induced nucleon EDM in previous calculations <cit.> 
received contribution (-2κ) from spurious mixing with the anomalous magnetic 
moment κ of the nucleon.
Fortunately, the correction is very simple and requires only the values of the nucleon 
anomalous magnetic moments from calculations on the same lattice ensembles.
Applying this correction properly is somewhat complicated due to differences
in the conventions used in these works.
Under some plausible assumptions we have demonstrated that, after the correction, 
even the most precise current lattice results for θ̅-nEDM may be compatible with zero.
If this finding is confirmed in detailed reanalysis of Refs. <cit.>, 
the precision of the current lattice QCD determination of θ̅-nEDM may be completely 
inadequate to constrain the QCD θ̅ angle from experimental data.
The entire modern Physics program to search for fundamental symmetry violations as signatures
of new physics relies on our understanding of the effects of quark and gluon  
interactions on nucleon structure.
The importance and urgency of first-principles calculations of these effects
hardly needs more emphasis, and we have to conclude, that they will likely be 
even more difficult that thought before.







In this paper, we have performed calculations of nucleon electric dipole moments 
induced by -odd quark-gluon interactions using two different methods. 
In the first method, we have successfully calculated the nucleon-current correlators 
with and without the -odd interaction, evaluating up to four-point connected nucleon 
correlation functions.
We have demonstrated that this novel technique works well and we argue that it is both cheaper 
and has fewer uncertainties than the technique used in
 <cit.> to compute the same observables with modified
Wilson action.
One of the obstacles to applying the technique of
Refs. <cit.> is that low-eigenmode deflation used to
accelerate calculations will be more expensive, because the eigenvectors have to be computed for 
every modification of the fermion action.
This may also be partially true for recently introduced multi-grid methods, in which
operator-dependent subspace null vectors have to be computed in the multi-grid setup phase, 
which has considerable cost.


In the second method, we computed the neutron EDM using its energy shift in uniform background
electric field and in the presence of the same -odd interaction.
The energy-shift method to compute nucleon EDM has been used before <cit.>,
but our calculation is the first one that uses the uniform background electric
field that respects boundary conditions <cit.>.
We perform calculations with identical statistics in both methods and can directly compare
the central values and the uncertainties of the results.
We find that the EDM results agree if the new formula for extraction of the EDFF F_3 is used. 
Also, both methods yield comparable uncertainty, and the energy shift method may be preferable in
the future because it does not require forward-limit extrapolation and the excited states may be
easier to control <cit.>.





Our calculations on a lattice are far from perfect and require improvement of the treatment of
excited states and forward-limit extrapolation of the form factors.
However, the associated systematic uncertainties are too small to cast doubt on the numerical
comparison of the energy shift and the form factor methods.
Although our calculations lack evaluation of the disconnected diagrams and 
renormalization and mixing subtractions of the quark chromo-EDM operator, 
these drawbacks apply equally to both methods, therefore do not affect said validation.


Future calculation of disconnected contributions to the F_3 form factors will 
be an extension to the present work, in which the quark-disconnected loops with 
insertions of the quark current, chromo-EDM, and both, will be evaluated and used 
together with the existing nucleon correlators.
The disconnected contractions do not require four-point correlators and are simpler 
to construct, although the stochastic noise will likely be a much bigger problem 
than for the connected contractions.
We expect that with advances in numerical evaluation of the disconnected
diagrams <cit.>, this problem will be tractable.






T.B. is supported by US DOE grant DE-FG02-92ER40716.
T.I. is supported in part by US DOE Contract AC-02-98CH10886(BNL). 
T.I. is also supported in part by the Japanese Ministry of Education 
Grant-in-Aid, No. 26400261.
H.O. is supported by the RIKEN Special Postdoctoral Researcher program
S.N.S. was supported by the Nathan Isgur fellowship program at JLab
and by RIKEN BNL Research Center under its joint tenure track
fellowship with Stony Brook University.
This material is based upon work supported by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics under contract DE-AC05-06OR23177.
The U.S. Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish
or reproduce this manuscript for U.S. Government purposes.
S.A. and S.N.S are also grateful for the hospitality of Kavli Institute for Theoretical Physics 
(UC Santa Barbara) during the “Nuclear16” workshop.
This research was supported in part by the National Science Foundation 
under Grant No. NSF PHY11-25915. 





Gauge configurations with dynamical domain wall fermions used in this work
were generated by the RBC/UKQCD collaboration.
The computation was performed using the Hokusai supercomputer of the RIKEN ACCC facility
and Jlab cluster as part of the USQCD collaboration.
The calculations were performed with the “Qlua” software suite<cit.>.







§ CONVENTIONS


In this appendix section, we collect conventions for γ-matrices implicitly or explicitly
used throughout the text.
In Table <ref>, we also provide notes the transformation between 
Minkowski () and Euclid () notations to avoid any ambiguities in matching
Minkowski and Euclidean form factor expressions for matrix elements and vertices.

In Minkowski space with metric {-1,-1,-1,+1}, we use the chiral γ-matrix basis

    [γ^i]_ = ([       σ^i; -σ^i      ]) ,   
    [γ^4]_ = ([   1; 1   ]) ,

and with ϵ^4123=+1 we define the chiral γ_5 matrix 

    [γ_5]_ 
      = -i/4![ϵ^μνρσγ_μγ_νγ_ργ_σ]_
      = i[γ^4γ^1γ^2γ^3]_
      = ([ -1   ;     1 ]) ,

For the spin matrix σ^μν=i/2[γ^μ,γ^ν] we 
will also need the relation 

    [σ^μνγ_5]_
      = i/2[ϵ^μνρσσ_ρσ]_

In accordance with Tab. <ref>, the γ-matrices in Euclidean space are

    [γ^i]_ = ([       -iσ^i; +iσ^i       ]) ,   
    [γ^4]_ = ([   1; 1   ]) ,

in which γ^1,3 have the opposite sign compared to the deGrand-Rossi basis used in most 
of the lattice QCD software.
This difference is inconsequential because all results are manifestly covariant with respect 
to unitary basis transformations.
Finally, we use the γ_5 definition that agrees with the lattice software,

    [γ_5]_ = [γ^1γ^2γ^3γ^4] = 
      ([  1   ;    -1 ]) ,

and note that the kinematic coefficients for vector form factors derived in
Sec. <ref> depend on a particular γ_5 definition in terms of
γ^μ, but the numerical lattice results are invariant as long as the same
[γ_5]_ is used in both Eqs. (<ref>)
and (<ref>).








§ ELECTRIC AND MAGNETIC DIPOLE MOMENTS AND FORM FACTORS


In this Appendix section, we recall the connection between the form factors F_2,3 and
the magnetic and electric dipole moments of a spin-1/2 particle.
Although this is discussed in many textbooks, we find it useful to perform a rigorous derivation
expanding the matrix element (<ref>) in the momentum transfer q=p^'-p
and taking the limit q→0.
For completeness and to avoid any ambiguities, in addition to the γ-matrices 
in Sec. <ref>, we collect all relevant conventions for E&M fields, 4-spinors, 
and their interaction.
The discussion in this Section assumes Minkowski conventions  with 
g_μν=diag{-1,-1,-1,+1}.

The fermion-photon interaction is determined by the form of the “long” derivative,

    D_μ = ∂_μ + ieA_μ,
       = ψ̅(iD_μγ^μ - m) ψ 
      = ψ̅(i∂ - m)ψ - e A_μ J^μ ,

which leads to the interaction Hamiltonian

    H_int = ∫ d^3x (-_int)
      = e∫ d^3x  A_μ J^μ = e ∫ d^3x  (ρϕ - J⃗·A⃗)

where the EM potential A^μ=(A⃗, ϕ), EM current J^μ=(J⃗, ρ),
and the electric coupling (charge) e=|e|. 

To evaluate the matrix element (<ref>) in the
interaction (<ref>), we use the chiral γ-matrix representation
summarized in Appendix <ref>.
The on-shell spinors satisfying the regular Dirac equation with a real-valued mass m>0 
and energy E^(')=√(m^2+p⃗^(')2)
take the form

    u_p 
         = ([ √(E-p⃗σ⃗)ξ; √(E+p⃗σ⃗)ξ ]) 
      = √(m)[1 + p⃗Σ⃗/2mγ_5 + O(p⃗^2)]
      ([ ξ; ξ ]) ,
    
    u̅_p^'   = ([ √(E^'-p⃗^'σ⃗)ξ^'; √(E^'+p⃗^'σ⃗)ξ^' ])^†γ^4
      = √(m)([ ξ^'; ξ^' ])^†[1 - p⃗^'Σ⃗/2mγ_5 + O(p⃗^2)]

where
Σ^k = 1/2ϵ^ijkσ^jk 
  = ([ σ^k    ;     σ^k ]) .
We will use these spinors to evaluate matrix elements of the Hamiltonian (<ref>),
treating the E&M field as classical background. 
Note that in order to treat these matrix elements as the interaction energy, 
they must be normalized as non-relativistic,

    E_int 
      = ⟨p⃗^',σ^'|H_int|p⃗,σ⟩_NR
      = eA_μ1/√(2E^'·2E)u̅_p^'Γ^μ u_p
      ≐ eA_μΓ^μ ,

where we introduced the notation 
X=1/√(2E^'·2E)u̅_p^' X u_p for convenience.
In the limit of small spatial momenta |p⃗|,|p⃗^'|→0, only the spatial components
σ^ij give non-vanishing contributions when contracted with
spinors (<ref>):

    σ^ij 
      = 1/√(2E 2E^')u̅_p^'σ^ij u_p 
      = ϵ^ijk ξ^'†σ^kξ + O(|p⃗|,|p⃗^'|) ,
      σ^4k 
      = u̅_p^'σ^4k u_p 
      = O(|p⃗|,|p⃗^'|) ,

 

Recalling the conventions <cit.> for the EM potential A^μ,

    (E⃗)^i 
         = -∂/∂ x^i A^4 - ∂/∂ t (A⃗)^i ,
    
    
    (H⃗)^i    = (curlA⃗)^i
      = ϵ^ijk∂/∂ x^j (A⃗)^k ,

which result in the following field strength tensor F_μν and its dual
F̃_μν=1/2ϵ_μνρσ F^ρσ,
ϵ_1234=+1

    F_μν   = ([         1    2    3    4;    1    0 -H^3  H^2 -E^1;    2  H^3    0 -H^1 -E^2;    3 -H^2  H^1    0 -E^3;    4  E^1  E^2  E^3    0 ]) ,
    
    F̃_μν   = ([         1    2    3    4;    1    0  E^3 -E^2 -H^1;    2 -E^3    0  E^1 -H^2;    3  E^2 -E^1    0 -H^3;    4  H^1  H^2  H^3    0 ]) ,

where the rows and the columns are enumerated by μ and ν, respectively.
With the following conventions for the fermion and photon fields with definite momenta
p^(') and q, respectively,

    ψ_p(x)∼ e^-ipx ,
      ψ̅_p^'(x) ∼ e^ip^' x ,
       A_q,μ(x)∼ e^-i(p^'-p)x = e^-iqx ,

the derivatives acting on these fields are translated into momentum factors,

    ∂ψ = γ^μ∂_μψ   →γ^μ (-ip_μ) ψ = (-i)pψ ,
    
    
    F_μν(x) = ∂_μ A_ν - ∂_ν A_μ   → (-i) ( q_μ A_ν - q_ν A_μ) .

Applying the Gordon identity to Eq. (<ref>) and omitting the F_A form
factor, we get 

    ⟨ p^',σ^' |J^μ|p,σ⟩_
      = u̅_p^'[
        F_1 (p^'+p)^μ/2m 
        + (G_M + iγ_5 F_3) iσ^μνq_ν/2m] u_p,σ ,

where G_M=F_1+F_2 is the magnetic Sachs form factor determining the full magnetic moment 
μ=Q+κ=G_M(0).
The first term is independent of the spin and is equal to the electromagnetic interaction 
of a scalar particle, which we omit as irrelevant.
With the use of (<ref>) and (<ref>), the spin-dependent 
interaction energy takes the form 

    E_int,spin 
      = i q_ν A_μ [
            eG_M σ^μν/2m
           - eF_3 1/2ϵ^μνρσσ_ρσ/2m]
      = 1/2(e G_M/2m  F_μν - e F_3/2m F̃_μν) σ^μν .

Neglecting all but the leading order in O(|p⃗|,|p⃗^'|), 
we only have to keep the spatial components σ^ij:

    E_int,spin 
      = -e G_M/2m H⃗·Σ̂- e F_3/2mE⃗·Σ̂ ,

where the unit spin vector Σ̂= ξ^'†σ⃗ξ, |Σ̂|=1.
The coupling coefficients to the magnetic and electric fields in the above equation have to be
identified with the magnetic and electric dipole moments, respectively,

    μ_N = G_M(0) ,
       d_N = F_3(0) .

which both are expressed here in the particle magneton units e/(2m).

Note that the above derivation could be repeated for the chirally-rotated spinors and the
nucleon-current vertex (<ref>). 
It can be easily shown that the only change compared to Eq. (<ref>) would be that the
magnetic and electric fields would couple to some orthogonal linear combinations of
F̃_2,3, and that these combinations would reproduce F_2 and F_3 exactly in
agreement with Eq. (<ref>).

Finally, we note that if one uses the chirally-rotated spinors to calculate the spatial 
matrix elements σ^ij, they are reduced by a factor of cos(2)
while the timelike matrix elements σ^4k become non-zero:

    e^2iγ_5σ^ij   = cos(2) σ^ij + sin(2) ϵ^ijkσ^4k ,
    
    e^2iγ_5ϵ^ijkσ^4k   = -sin(2) σ^ij + cos(2) σ^4k .

As we noted above, u̅_p^'σ^iju_p couples to the magnetic field, 
while u̅_p^'σ^4ku_p couples to the electric field.
This “mixing” of electric and magnetic fields compensates exactly the mixing in
Eq. (<ref>) induced by using the chirally-rotated spinors
ũ̅, ũ_p instead of the regular spinors u̅_p^', u_p.






§ KINEMATIC COEFFICIENTS


In this section, we present expressions for the kinematic coefficients for form factors
F_1,2,3 on a Euclidean lattice.
We use two types of the polarization projectors, (1) spin-average T^+ and (2) polarized
T^+_S_z. 
Both projectors also select the upper (positive-parity) part of the nucleon spinors

    T^+=[1+γ^4/2]_ ,  
    T^+_S_z=[1+γ^4/2(-iγ^1γ^2)]_ .


Form factor F_3 can be extracted from μ=3,4 components of the vector current matrix
elements between S_z-polarized nucleon states.
Using handy notations for the positive-parity nucleon spinor matrices,

    = -ip_ + m
     ,  ^' = -ip_^' + m ,

the form factor expression for the  nucleon-current correlation function on a lattice 
C^_3pt=C^_NJ^μN̅ can be written as

    Tr[T_pol C^_N J^μN̅(p⃗^',t; q⃗,t_op)]
         =e^-E^' (t-t_op) - Et_op/2E^'·2E[ e^iγ_5 T e^iγ_5·^'·Γ_^μ(p^',p)·]
         =e^-E^' (t-t_op) - Et_op/2E^'·2E[ (T + i{γ_5,T} + O(^2)) ·^'·Γ_^μ(p^',p)·] ,

where, assuming that the -odd interaction is small, we have expanded in
the -odd mixing angle .

Below we quote formulas for contributions to the last line of Eq. (<ref>)
computed for zero sink momentum p⃗^'=0, 

    source p⃗   = - q⃗ ,     E   =√(m^2+q⃗^2) ,
    sink p⃗^'   = p⃗ + q⃗ = 0 ,     E^'   =m ,

with the nucleon spin projectors T^+ and T^+_S_z.
The -independent contribution is 

    [T^+^'Γ_^μ] 
       = 4m^2([   i q_1 / m -iτ q_1 / m           0;   i q_2 / m -iτ q_2 / m           0;   i q_3 / m -iτ q_3 / m           0;      2(1+τ)    -2τ(1+τ)           0 ]) ,
    
    [T^+_S_z^'Γ_^μ] 
       = 4m^2([         -q_2 / m         -q_2 / m q_1 q_3 / (2m^2);          q_1 / m          q_1 / m q_2 q_3 / (2m^2);                0                0   q_3^2 / (2m^2);                0                0  -i(1+τ) q_3 / m ]) ,

where the rows correspond to the Lorentz components μ=1,2,3,4 
and the columns correspond to the form factors F_1,2,3.
We have also introduced the frequently used kinematic variable τ

    τ≐Q^2/4m^2 p⃗^'=0≡ E-m/2m .

The coefficients of the contributions ∼ are

    [{γ_5,T^+}^'Γ_^μ] 
       = 4m^2([          0          0 -τ q_1 / m;          0          0 -τ q_2 / m;          0          0 -τ q_3 / m;          0          0   2iτ(1+τ);            ]) ,
    
    [{γ_5,T^+_S_z}^'Γ_^μ] 
       = 4m^2([                       0      i q_1 q_3 / (2m^2)                       0;                       0      i q_2 q_3 / (2m^2)                       0;                    -2iτ -2iτ + i q_3^2 / (2m^2)                       0;                -q_3 / m               τ q_3 / m                       0;                         ]) .

Up to order O(), the  nucleon correlation functions are[
  Note that both  and F_3 are proportional to the -odd perturbation, therefore we
  consider F_3=O() and drop terms F_3 and higher.
] 

    [T^+ C^_N J^3 N̅]
         = [
          iq_3/m G_E
        + O(^2) ] ,
    
    [T^+ C^_N J^4 N̅]
         = [
          2(1+τ)G_E
        + O(^2)] ,
    
    [T^+_S_z C^_N J^3 N̅]
         = [
          2τ G_M
        - q_3^2/2m^2F_2
        + q_3^2/2m^2 F_3 
        + O(^2)] ,
    
    [T^+_S_z C^_N J^4 N̅]
         = [
        - iq_3/m G_E
        - i(1+τ)q_3/m F_3 
        + O(^2)] ,

where G_E=F_1-τ F_2 is the electric and G_M=F_1+F_2 is the magnetic Sachs form
factor, and 

    =m/Ee^-E^' (t_sep-t_op) - E t_op

is the time dependence combined with kinematic
factors.
In the analysis of the C_NJN̅/C_NN̅ ratios (<ref>), the exponential
time dependence is canceled, and the kinematic coefficients have to be modified 
to take into account the traces of the nucleon two-point functions:

    _=m/√(2E(m+E)) .



In addition, we evaluate the extra contributions to the kinematic coefficients
∼{γ_5,Γ_^μ} that comes from spurious mixing of F_2,3 

    [T^+^'{γ_5,Γ_^μ}] 
       = 4m^2([          0          0 2τ q_1 / m;          0          0 2τ q_2 / m;          0          0 2τ q_3 / m;          0          0 -4iτ (1+τ);            ]) ,
    
    [T^+_S_z^'{γ_5,Γ_^μ}] 
       = 4m^2([                0 -i q_1 q_3 / m^2     -2 i q_2 / m;                0 -i q_2 q_3 / m^2      2 i q_1 / m;                0   -i q_3^2 / m^2                0;                0 -2 (1+τ) q_3 / m                0;                  ]) ,

which in Refs. <cit.>
contributes to the polarized nucleon-current correlators as

    δ[T^+_S_z C^_N J^3 N̅]
         ?= [
          q_3^2/m^2F_2
        + O(^2)] ,
    
    δ[T^+_S_z C^_N J^4 N̅]
         ?= [
        - 2i(1+τ)q_3/m F_2
        + O(^2)] ,

If the terms (<ref>,<ref>) 
are erroneously added to the kinematic 
coefficients (<ref>,<ref>), analysis of the same 
lattice correlation functions will result in incorrect values of EDFF 
F̃_3 = F_3 - 2 F_2, in full agreement with Eq. (<ref>).




aip


