On the fermionic couplings of axionic dark matter

In the non-relativistic limit, two types of dark matter axion interactions with fermions are thought to dominate: one is induced by the spatial gradient of the axion field and called the axion wind, and the other by the time-derivative of the axion field, generating axioelectric effects. By generalizing Schiff theorem, it is demonstrated that this latter operator is actually strongly screened. For a neutral fermion, it can be entirely rotated away and is unobservable. For charged fermions, the only effect that can peek through the screening is an axion-induced electric dipole moment (EDM). These EDMs are not related to the axion coupling to gluons, represent a prediction of the Dirac theory analogous to the g = 2 magnetic moments, are not further screened by the original Schiff theorem, and are ultimately responsible for inducing the usual axioelectric ionization. The two main phenomenological consequences are then that first the axion-induced neutron EDM could be several orders of magnitude larger than expected from the axion gluonic coupling, and second, that the electron EDM could also become available, and could actually be highly sensitive to relic axions.


Introduction
The axion mechanism [1,2] is currently our best solution to the strong CP puzzle.The nonobservation of a neutron electric dipole moment (EDM) constrains the QCD theta term, θG µν Gµν to be tiny, θ ≲ 10 −10 [3].As this coupling receives contributions from two unrelated sectors of the Standard Model (SM), a topological QCD contribution and an electroweak contribution from the quark Yukawa couplings, both a priori of O(1), such a tiny value requires an unacceptable fine-tuning.The axion solution [4,5] relies on the axion a being the Goldstone boson associated to the spontaneous breaking of an anomalous U (1) symmetry, the PQ symmetry [1,2].This ensures a coupling of the axion to gluons, aG µν Gµν , which develops into a potential for the axion in the low-energy limit.At the minimum of this potential, the axion field absorbs the θ term, making it unobservable.
The initial implementation of the axion mechanism relied on the axion emerging at the electroweak scale, and was quickly ruled out as this would imply way to large couplings to matter particles.Invisible axion scenarios were then developed, most notably the DFSZ [6,7] and KSVZ [8,9] models, in which the axion is very light and very weakly coupled.In most cases, being in addition very long-lived, the axion emerged as a viable dark matter (DM) candidate (for a review of the axion in a cosmological context, see e.g.Ref. [10]).Interestingly, experimental strategies could then take advantage of the rather high flux of such dark matter axions (see e.g.Ref. [11,12]).In practice, dark matter axion production mechanisms ensure the axion is rather cold, and being in addition very light, it can be represented by a classical coherent pseudoscalar field, typically a(r, t) = a 0 cos(Et − p • r), E 2 = p 2 + m 2 a , m a the axion (or axion-like particle) mass, and a 0 set by the local DM density, m a a 0 = √ 2ρ DM with ρ DM ≈ 0.4 GeV/cm 3 [13].
The goal of the present paper is to analyze the couplings to SM fermions of such a dark matter axion background, in the non-relativistic limit.The usual starting point is the axion Lagrangian (to simplify the notation, a coupling constant g = m/Λ with Λ the PQ breaking scale is understood to be absorbed into a throughout this paper) Such a derivative interaction to the axion is reminiscent of its Goldstone boson nature.The corresponding Dirac equation is i∂ t |ψ⟩ = H D |ψ⟩ with where ȧ = ∂ t a.In the Dirac representation, where γ 0 is diagonal, γ 5 directly couples the fermion and antifermion degrees of freedom of |ψ⟩.Consequently, in the non-relativistic limit, the ȧ term receives a dependence on p = −i∇: These two leading interactions have been extensively studied in the literature [14,15].The so-called axion wind term, γ 5 γ • ∇a, leads to a coupling of the gradient of the axion field to the spin of the fermion.It can be searched for experimentally e.g. using NMR techniques [16][17][18][19] or magnons [20].
The second term is dubbed the axioelectric effect [21][22][23][24].It translates as a coupling of ȧ to the combination p • S of the momentum p and spin S of the fermion.As a result, sufficiently energetic axions could kick bound electrons out, in analogy with the photoelectric effect.The sun could produce a consequent flux of such axions, whose possible detection via these ionization processes, or more generally electron recoil effects, gave rise to a rather intense experimental activity [25][26][27][28][29][30][31][32].The corresponding constraints on the axion are reviewed e.g. in Ref. [33], as well as more recently in Refs.[34,35] in the context of the excess events observed at XENON1T [36].Note, though, that these experiments also probe different mechanisms and/or the coupling of the axion to photons.
A peculiar feature of Goldstone bosons is that there are different ways to parametrize them.For the axion, an equally valid Lagrangian uses the so-called polar or exponential parametrization: The derivative interaction is replaced by an infinite tower of interactions, starting by the pseudoscalar coupling a ψiγ 5 ψ.The corresponding Hamiltonian is then and its non-relativistic limit can be worked out to be The same axion wind and axioelectric interactions emerge, but the coefficient of the latter differs by a factor two. Historically, this fact was well known in the context of nucleon-pion interactions.
The equivalence of the pseudoscalar and axial interaction was first discussed by Dyson in 1948 [37] (see also Ref. [38,39]), on the basis of the axion wind term being the same.Later, this ambiguity in the time-dependent term, as well as in some higher order terms in the non-relativistic expansion, generated a lot of attention [40][41][42][43][44][45].As we will see, part of the issue was related to the truncation of the exponential parametrization.After all, many of these works date back to before Goldstone theorem was formulated, let alone the pion identified as a pseudo-Goldstone boson of the chiral symmetry breaking.Nowadays, the equivalence between the derivative and exponential representation is an established fact, but surprisingly, a non-relativistic expansion truly reflecting this has not been worked out yet.This is the purpose of the present paper.
In particular, adopting a modern language, we will see that the γ 5 {γ • p, ȧ} coupling can be systematically rotated away for a neutral fermion.The demonstration is actually quite simple and can readily be given.First, remember that a non-relativistic expansion is not unique 1 .As customary in quantum mechanics, unitary transformations cannot change the physics.So, performing such a transformation, and provided the block-diagonal nature of the Hamiltonian is maintained, an equally valid non-relativistic expansion is found.Now, as proposed a long time ago in Ref. [40,42], consider |ψ⟩ → ψ ′ = exp(iS) |ψ⟩ , S = µ 4m 2 γ 5 {γ • p, a} . If ).Thus, acting on H NR D with µ = 2, or on H NR E with µ = 1, the γ 5 {γ • p, ȧ} coupling is replaced by O(1/m 3 ) and higher terms.
At the time, this was interpreted as an ambiguity that should cancel out in physical observables.Here, we will go one step further and argue that γ 5 {γ•p, ȧ} does not encode any true physical effects.In other words, for neutral fermions, the operator γ 5 {γ • p, ȧ} is totally screened at O(1/m 2 ) in the non-relativistic expansion.One reason for this interpretation has to do with Schiff's theorem [47], which states that charged fermion EDMs are screened.The transformation S of Eq. ( 7) is closely related to Schiff's transformation, and even becomes the Schiff's transformation for a charged fermion.The consequence in that case is that the covariant γ 5 {γ •P, ȧ} coupling becomes equivalent to an axion-induced EDM operator, aσ • E at O(1/m 2 ).Phenomenologically, these two operators are indistinguishable, ensuring that physics is independent of the choice of parametrization.Yet, this equivalence makes it manifest that unexplored but promising avenues do exist to search for dark matter axions.
The paper is organized as follow.To set the stage, we start in the next section by a brief overview of the construction of the non-relativistic expansions via the Foldy-Wouthuysen method [46].This also gives us the opportunity to introduce Schiff's theorem [47] and its generalizations.Then, in Sec. 3, we enter the core of the subject, and perform the non-relativistic expansion of the axionic Hamiltonian up to and including O(1/m −3 ) terms, firstly in the absence of electromagnetic (EM) fields, secondly for a charged fermion minimally coupled to EM fields, and thirdly for a neutral fermion having electric and magnetic dipole interactions with the EM fields.These results are then put to use in Sec. 4 to analyze axion-induced lepton and nucleon EDMs, showing how and when some new effects could be expected.Finally, in Sec. 5, our results are summarized along with their phenomenological consequences.

Brief overview of the non-relativistic expansion
The techniques used in the present paper are covered in most textbooks on relativistic quantum mechanics.In particular, recovering the Pauli equation by a non-relativistic expansion of the Dirac Hamiltonian for a spin 1/2 field minimally coupled to EM fields, where P = p − eA, p = −i∇, and the EM potential is A µ = (ϕ, A), is a standard exercise.Though well-known, we think it is nevertheless useful to briefly review this so as to fix our notations, and because it forms the backbone on which we will add axions later on.Further, once the magnetic moment and electric moment operators σ µν F µν and σ µν F µν are added, it permits to introduce the Schiff's theorem that will be central to the axion discussion.

Foldy-Wouthuysen transformation
The Dirac equation involves four-dimensional spinors, and thus includes both particles and antiparticles simultaneously.In the non-relativistic limit though, the energy is not sufficient for pair creation, and the antiparticle degrees of freedom are not dynamical.In practice, the Dirac equation must reduce to a decoupled pair of two-dimensional Pauli equations, describing the dynamics of spin 1/2 particles only.Several procedures exist to perform this reduction, starting historically by Pauli's elimination method [48].To set the stage, let us briefly describe the main idea.We first adopt the Dirac representation for the gamma matrices, that is, where σ are the usual Pauli matrices.Note also the identities γ i γ j = (−δ ij − iε ijk σ k )1 and σ ⊗ 1 = − γ 0 γ 5 γ, as well as the fact that γ † = γ 0 γγ 0 = −γ, but γ 0 † = γ 0 and γ 5 † = γ 5 .The diagonal form of γ 0 is instrumental for performing the non-relativistic expansion.Indeed, if the Dirac spinor |ψ⟩ is split into a pair of two-component spinors, the Dirac equation Eq. ( 8) takes the matrix form (after χ → −χ) Factoring out the large time evolution due to the rest mass and defining This is the essence of Pauli elimination method that can be generalized to the presence of other interactions and to higher orders.In those cases though, the method becomes very cumbersome because hermiticity of the Hamiltonian is not guaranteed, and additional renormalizations of the φ field are in general required [49].The Foldy-Wouthuysen (FW) procedure is designed to systematize the block-diagonalization of the Dirac Hamiltonian [46].Starting from i∂ t |ψ⟩ = H |ψ⟩, the idea is to construct a unitary rotation ψ → ψ ′ = e iS ψ such that i∂ t |ψ ′ ⟩ = H ′ |ψ ′ ⟩ with with H ′ now block-diagonal.This decouples the large two-component spinor from the small one, and should be valid as long as the energy involved does not allow for pair creation.Since we started by performing a unitary transformation, there is no hermiticity issue with H ′ .However, in the presence of interactions, an exact solution for S cannot be found in general, and one relies instead on a perturbative expansion in c −1 .That is, instead of a single unitary transformation S, a sequence of unitary transformations is performed to bring H ′ to a block diagonal form, up to some order c −n .For dimensional reasons, an expansion in 1/c is essentially identical to an expansion in 1/m, so we will rather concentrate on the latter and keep c = 1.Details of this construction are in Appendix A. In summary, one first uses the diagonal γ 0 to write the Hamiltonian as where O stand for odd terms, Oγ 0 = −γ 0 O, and E for even terms, Eγ 0 = γ 0 E. In general, O and E are differential operators that do not commute.The term O is the offending one that couples small and large components.So, in the first step, we must remove it by some unitary transformation S.
Since to leading order, Performing that transformation cancels the O term in H, but brings back odd terms at higher orders (proportional to [O, E], O 3 , etc), so the procedure must be iterated up to some given order in 1/m.After three steps, the Hamiltonian becomes where When applied on a four-component spinor, the upper two and lower two components are decoupled.Given the choice of γ 0 , only the large upper component needs to be kept, as the lower small component dynamics is dampened by the rest mass, i.e., by a P/m factor.
The FW transformation will be the first step in all our developments.Yet, it is important to stress that it is not the end of the story.As was realized comparing various block-diagonalization methods, including the elimination method, there are some ambiguities in the final form of H NR .This simply reflects the fact that additional unitary transformations ψ → ψ ′ = e iS ψ are still allowed as long as S is even (for a review, see e.g.Ref. [49]).This feature, at the root of Schiff's theorem, will be used extensively in the following.

Application to electromagnetic interactions
Taking O = γ • (p − eA) ≡ γ • P and E = eϕ, and keeping only terms up to O(1/m 3 ), the standard result is recovered: By convention, ∇ acts on the quantity immediately to its right, but P acts on everything.Note that the σ matrices occurring are to be interpreted as 1 ⊗ σ, since this Hamiltonian still acts on four-dimensional spinors.Yet, being diagonal, the reduction to the Pauli equation is now trivial.
As is well-known, one can identify the Zeeman magnetic coupling σ • B = 2S • B with S the spin operator and g = 2 the magnetic moment, the spin orbit coupling σ • (E × P), and the Darwin term ∇ • E.
A more interesting application starts by including the higher order magnetic moment and electric moment operators where electromagnetic fields satisfy 15), keeping in mind that δ µ and d are O(m −1 ), and discarding terms of O(m −4 ) and higher, the block-diagonal Hamiltonian is now The magnetic operator thus describes the deviation of the magnetic moment from its Dirac value, a = (g − 2)/2.The σ • E term describes the electric dipole interaction, with d the EDM.If we remember the identities the Darwin and spin-orbit couplings are identified inside [γ • P, γ • E], now modified by a magnetic moment contribution and accompanied by magnetic interactions induced by d.

Schiff 's theorem and beyond
As stated before, the FW transformed Hamiltonian can still be unitarily rotated without breaking its block-diagonal character.The simplest such transformation is The transformation exp(iS 1 ) is unitary, and importantly, it commutes with the mass term γ 0 m.One should not be put off by the fact that this transformation involves the external fields via P.
Actually, we already did many such transformations to block-diagonalize the Hamiltonian, since the first FW transformation is exp(iS) with iS = O/(2m) and O = γ • P + iγ 0 γ • (δ µ E + dB).All that differs here is the γ 5 factor, making S 1 even with respect to γ 0 .The new Hamiltonian H ′ = e iS 1 (H − i∂ t ) e −iS 1 can be expanded as before, and since iS 1 ∼ O(m −1 ), we need to compute: Now, the key in Schiff's theorem [47] is to note that the O(1/m) terms miraculously combine as Thus, with α = md/e, the EDM term in H NR EM is rotated away!More accurately, we should say that it is transformed into higher order corrections coming from the rest of Eq. ( 21).After some algebra, the transformed Hamiltonian is found to be, keeping only terms at most linear in either a or d, since these quantities are experimentally small, where the only non-trivial reduction is [γ•P, P 2 ] = −eγ 0 γ 5 [γ•P, γ•B], the rest being straightforward algebraic manipulations.

Higher order Schiff transformations and operator redundancies
Central to Schiff's theorem is the presence of the − Ṡ1 piece that directly enters in the transformed Hamiltonian in Eq. ( 21), and can thus directly interfere with the other terms.When applied on γ 5 γ • P, it generates a γ 5 γ • Ȧ term out of which γ 5 γ • E emerges without an additional m −1 factor.This same trick can be used for any term that involves time derivatives of external fields.For instance, consider now Since it is already of O(m −3 ), only the leading commutator with eϕ needs to be computed.Again, [iS 2 , eϕ] combine with the Ṗ in − Ṡ2 to give a ∇ϕ + Ȧ = −E factor: This time though, we find a redundancy among O(1/m 3 ) operators, up to higher order corrections.Our preferred choice is to take β = 1 to get rid of the Ė and Ḃ operators, but one could equally well decide to keep the Ė operator and eliminate the E 2 term, or keep the Ḃ term and eliminate the E • B couplings.A third possible transformation is which also introduces a redundancy among O(1/m 3 ) operators, up to higher order corrections, These redundancies can be used to reduce the number of relevant operators.In Appendix B, we present one possible choice of S 2 and S 3 that bring H NR EM to a somewhat optimal form.It should be stressed though that the final coefficients for the higher order operators in H NR EM should not be taken too literally.Indeed, once adopting an effective description with the O(1/m) couplings σ µν F µν and σ µν Fµν , one could in principle also include O(1/m 2 ) or O(1/m 3 ) operators.For example, if one adds the F µν F µν or F µν F µν operators in H EM , their coefficients will directly correct those of Finally, it is worth to stress that this list certainly does not exhaust possible unitary transformations, and that not all such transformations encode useful information.For example, consider which is even and hermitian for η real.The change in the Hamiltonian is This transformation just adds the {γ • P, γ • E} operator to the Hamiltonian, up to higher order terms.To understand why this has no impact on the physics, let us first expand it using Eq.(19), If we could take η imaginary, this operator would interfere with the Darwin and spin-orbit operator i[γ • P, γ • E], but this would make exp(iS 4 ) non-unitary.Actually, this operator has no impact because ∇ • E and 2E • ∇ compensate each other when acting on wavefunctions (both are standard forms for the Darwin operator), while the σ • (∇ × E) = −σ • Ḃ term drops out for static fields (and could be rotated away by a dedicated unitary transformation with S 5 ∼ σ • B anyway).

Schiff theorem and charged fermion EDMs
Schiff's theorem shows that the energy of a charged particle cannot be influenced by its EDM at leading order.The naive interpretation of this result is that a charged particle plunged in an electric field would feel the Lorentz force and fly away.The Schiff's transformation is then viewed as a translation that moves us in some sort of rest frame for the charged fermion in which there is no electric field anymore, hence where the EDM operator vanishes and cannot contribute to Stark energy shifts.Thus, for charged fermions, [γ • P, γ • B] encodes the leading impact the EDM has on the particle energies in the non-relativistic limit.Using Eq. ( 19), one can recognize in this term the spin-dependent σ • (p × B) coupling discussed originally by Schiff [47].To feel the EDM with electric fields, one has to go fetch the O(1/m 3 ) operators γ 5 {γ • P, {γ • P, γ • E}} ] which, thanks to Eq. ( 27), can both be reduced to γ 0 {P 2 , σ • E}.This operator thus encodes the leading relativistic corrections.For the case of an electron in a heavy atom, significant enhancements of this operator have been found that guarantee an experimental sensitivity to the electron EDM [50].Finally, it should be mentioned that another way to evade the shielding of the EDM is to account for finite-size effects, that clearly go beyond the current formalism (for a review, see Ref. [51,52]).Schiff's theorem is a statement about the Hamiltonian, and thus applies to the energy levels of bound charged fermions.It does not mean that the EDM operator cannot be felt using other observables.In particular, the electric field does exert a torque on the spin of charged fermions, leading to its precession in adequate experimental settings.Specifically, to leading order, the spin operator S = σ/2 evolves according to which is nothing but one term of the generalized Bargmann-Michel-Telegdi equation [53,54].After the Schiff rotation, the dσ • E operator is removed from H, and the spin operator in that basis satisfies The crucial point that makes this equation compatible with that of S is that the spin operator does not commute with the Schiff transformation, so that S ′ = e iS 1 Se −iS 1 ̸ = S. Plugging this in the above equation implies which obviously holds by construction, since [iS 1 , eϕ] − Ṡ1 = −γ 0 dσ • E, see Eq. (22).This exercise provides another interpretation of Schiff's theorem.In a gauge in which ϕ = 0, S ′ appears constant in time, so the Schiff rotation of Eq. ( 20) is actually that to the rotating frame in which the spin appears static.This could have been guessed from the start if one notes that S 1 actually involves the helicity operator, γ 5 γ • P = −γ 0 σ • P.

Schiff theorem and neutral fermion EDMs
Schiff's theorem cannot apply to neutral fermions.Indeed, one can simply send e → 0 to decouple EM fields in Eq. ( 17) while keeping an explicit EDM term, but the parameter of the Schiff's transformation in Eq. (20) has to be set to α = md/e, which is undefined in that limit.Said differently, it is only through a delicate interplay with the couplings to the external EM fields that the Schiff's transformation can interfere destructively with the EDM term.Thus, for the neutron, all one can do is to eliminate the Ė and Ḃ couplings, and starting from Eq. ( 17) in the e → 0 limit, one ends up with It is not possible to rotate away the EDM.At the fundamental level, d is induced by all the CPviolating operators involving gluons and/or quarks (for a review, see e.g.Ref. [55]).The most important contribution is that of the θ term, at the root of the strong CP puzzle, and which is estimated as [56,57] In the SM, the CKM contribution is negligible, but in principle, some New Physics may also induce fundamental EDMs for the quarks (see e.g.Ref. [58] and references cited there).As those are certainly far from non-relativistic inside a neutron, Schiff's theorem should be largely evaded.In a SU (6) model, the neutron EDM receives then the additional contribution [55] Note that the same rather naive model gives , in fairly good agreement with the measured δ µ = −1.913e/(2mN ).Though this hardly suffices to justify Eq. ( 36) as there is no analog of Schiff's screening for the magnetic moment, it is in fairly good agreement with recent lattice calculations [59] This shows that Schiff's screening theorem does not apply to quarks, as could have been expected since those are bound not by the electromagnetic interactions but by the strong interactions.
3 Axion interactions in the non-relativistic limit Nowadays, the equivalence between the pseudoscalar and derivative axial interactions is understood as particular application of the general reparametrization theorem to Goldstone bosons [60].Let us recall the essence of the argument (see Ref. [61] for more details).For a typical axion model, one starts with a spontaneously broken chiral symmetry, U (1) P Q .Then, the statement that Goldstone boson a must interact derivatively leads to the unique interaction term (remember that g = m/Λ is absorbed into a): Obviously, this interaction is invariant under constant shifts of the Goldstone field.This is the standard form for most axion analyses, but one should emphasize that the Goldstone field is actually parametrized non-linearly in this representation, since its simple shifts a → a + mθ must span the vacuum manifold as θ varies between 0 and 2kπ.Another peculiarity of this representation is that the fermion field ψ ends up neutral under the original U (1) P Q symmetry, even though it must initially be charged otherwise it would not occur in the Noether current of the U (1) P Q symmetry.This interpretation is made clearer by adopting a different parametrization for the fields, that in which the fermion field keeps its original charge.This is achieved via a chiral reparametrization of the fermion field When plugged in L D , the derivative coupling is replaced by a tower of pseudoscalar interactions A trivial mass term m ψψ necessarily breaks the axial symmetry U (1) P Q , so the fermion mass must arise through the symmetry breaking itself, like in the SM.Such a mass term is still invariant under the original chiral symmetry because the phase the fermion field acquires under U (1) P Q , ψ → exp(iθγ 5 )ψ, is compensated by the shift of the Goldstone field a → a + mθ.Now, to leading order in a, this interaction produces a pseudoscalar coupling of the fermion to the axion: Truncating the theory in this way, one should remember that O(a 2 ) terms and above are neglected.This approximation is only valid for on-shell fermions, since by integration by part, ψ(γ µ γ 5 ∂ µ a)ψ = − ψ(2imγ 5 a)ψ upon enforcing the free equation of motion (i As mentioned in the Introduction, part of the historic controversy on the equivalence between the axial and pseudoscalar descriptions of nucleon-pion interactions has to do with this truncation.Nowadays, the equivalence between both representations is built in chiral effective theories.For axion models, it is not always fully embedded yet, as we will see.Further, additional care is needed because the U (1) P Q symmetry being anomalous, so is the chiral reparametrization Eq. (38).As analyzed in details in Ref. [61] (see also Refs.[62,63]), the two representations are then equivalent only up to the presence of specific anomalous contact interactions of the axion to gauge bosons.In the present section, these effects are not relevant and will not be discussed further, but we will come back to them when analyzing the passage from the quark to the nucleon level in Sec. 4.
Our goal is to construct and analyze the axion-fermion interactions in the non-relativistic limit.To treat both representations simultaneously, we adopt the trick proposed by Friar a long time ago and described in Ref. [45].Specifically, let us start from L D .The Euler-Lagrange equation gives Then, we partially perform the fermion reparametrization Eq. ( 38), which is nothing but a unitary transformation Calculating µ) with the help of exp(iαγ 5 ) = cos α + iγ 5 sin α, we find So, this form permits to interpolate between the exponential and derivative representations, with H(0) = H D and H(1) = H E .Let us now perform the non-relativistic expansion of this expression, first as it stands, and then adding electromagnetic interactions.

In the absence of EM fields
In the Hamiltonian Eq. ( 43), the terms γ•p and γ 5 γ•∇a = γ 5 [γ•p, a] are diagonal since γ 5 γ = −γ 0 σ, but the γ 5 piece coming from the exponential and γ 5 ȧ are not.Splitting the exponential using exp(iαγ 5 ) = cos α + iγ 5 sin α, the elements to be used for the FW transformation are with S a ≡ m sin(2µa/m) = 2µa+... and C a ≡ m 2 (cos(2µa/m)−1) = −2µ 2 a 2 +....The calculation, though cumbersome, does not present any particular difficulty and we find with The non-derivative sine and cosine terms are singled out in the last line of Eq. ( 45) because they can be dropped.The specific combination S 2 a + 2C a already gives a term of O(a 4 ) and, when combined with the O(1/m 3 ) terms, gives the totally negligible contribution So, even though the polar representation initially involves non-derivative operators in a n , n > 1, none of them survive in the non-relativistic limit.This fact was not realized in Ref. [45], where only terms linear in the pseudoscalar field were kept.At this stage, we recover the expression in Eq. ( 3) and ( 6) by setting µ = 0 and µ = 1, respectively.As stated there, the axion wind term is independent of the parametrization, and On the other hand, the time-derivative term is not [45] This coupling even disappear for the specific choice µ = 2. Since there are no other O(1/m 2 ) terms, for this to make sense, this operator must not embody any real physical effects.

Schiff 's transformations
Let us first concentrate on the O(1/m 2 ) terms.In analogy with the transformation done in Sec.2.3 to eliminate the EDM operator, we can perform the unitary transform ψ → e iS 1 ψ with [40,42] This transformation is unitary and commutes with the mass term, [exp(±iS 1 ), γ 0 m] = 0.This means that, with iS 1 ∼ O(m −2 ) and H NR (µ) − γ 0 m ∼ O(m −1 ), only the first term of the expansion needs to be kept Explicitly, plugging in the expression of S 1 , The Ṡ1 term cancels precisely the O(1/m2 ) terms, by construction.That is Schiff's theorem trick in action.What it means is that this operator is actually a higher order effect, now embodied in the O(m −3 ) operators.In other words, we have succeeded at replacing the O(1/m 2 ) terms involving time-derivatives by O(1/m 3 ) terms involving only space derivatives, that is, axion wind operators.At this stage, it is clear that Schiff's trick can be used to remove or simplify the terms in H 3 involving time derivatives 2 .Specifically, we can perform to remove the term γ 0 γ 5 [γ • p, ä], up to some O(1/m 4 ) contributions.The final transformation we consider presents us with an alternative.Let us now rotate H NR (µ) with This kills the äS a coupling and corrects the ȧ Ṡ(a) in precisely the right way to make it µ independent: Now, we could have done the opposite, that is, make the äS a coupling µ independent by removing entirely the ȧ2 coupling.This time, the Schiff's transformation is not removing an operator, but telling us that two of them are redundant, up to higher order corrections.

Final Hamiltonian in the non-relativistic limit
All in all, after the unitary transformations S 1 in Eq. ( 50), S 2 in Eq. ( 53), and S 3 in Eq. ( 54), and after expanding S a and C a and keeping only terms up to O(1/m 3 ), the Hamiltonian becomes with where the last term with the 16/3 coefficient comes from the expansion of the O(m −1 ) term involving ∇S a , see Eq. (48).At this stage, algebraic manipulations of H 3 using commutator and anticommutator identities, e.g., permit to show that its O(a 3 ) terms cancel out completely, and its O(a 2 ) and O(a) terms become independent of µ.The final Hamiltonian is very simple and contains only five non-trivial operators: with the further information that ȧ2 can be freely traded for äa.Remember that the axion coupling constant has to be put back by a → ga with g = m/Λ and Λ the PQ breaking scale.Three comments are in order.
• It is remarkable that all µ dependences have cancelled out, and this involved highly nontrivial cancellations.In our opinion, it shows that the essential physical content is correctly identified, and redundancies kept at a minimum.Interestingly, this Hamiltonian cannot be obtained by setting µ to some value in H NR (µ) of Eq. ( 45).This is evident since S 1 , S 2 , and S 3 do not all vanish for the same value of µ.Said differently, the sequence of Schiff transformations S 1 , S 2 , and S 3 does not trivially undo the original Dyson rotation of Eq. ( 42).Note though that in practice, setting µ = 2 in Eq. ( 45) already goes a long way since S 1 has the most impact but vanishes for that value, at least for operators up to O(m −3 ).
What is different though is that we do not expect significant violations of this screening.First, finite size effects were relevant for the EDM as the electric charge density is far from constant in atomic systems.By contrast, the axion background should be relatively homogenous, even on macroscopic scales.Second, relativistic corrections were found significant for the EDM.But, as discussed in Sec.2.3, the relativistic corrections to γ • E were embodied in the very similar {P 2 , γ • E} operator.Here, the relativistic corrections replacing γ 5 {γ • p, ȧ} are totally different in nature: they all involve the axion wind and even vanish if ∇a = 0 (note that p 2 , {γ • p, a} = {γ • p, [p 2 , a]}).In that ∇a = 0 scenario, the relativistic corrections replacing γ 5 {γ • p, ȧ} would at best arise at O(m −4 ).For these reasons, we expect the screening of γ 5 {γ • p, ȧ} to be particularly effective.
• The leading fermionic coupling in a ∇a = 0 scenario is ȧ2 /(2m3 ), which is not a relativistic correction to γ 5 {γ • p, ȧ} but a genuine independent coupling.In this case though, being quadratic in the axion field, it is presumably totally negligible, and better windows could exist.In particular, in most scenarios, the axion also couples to photons.Classically, the aF µν F µν coupling can generates a ȧB term that act as a current density [64].
• On a technical note, let us stress that it is crucial to use the full exponential parametrization to correctly identify the final operators.Had we truncated the polar representation to its leading term by setting S a = 2µa and C a = 0, not only would there still be O(a 3 ) operators in the final Hamiltonian, but the µ dependence would not have cancelled completely [45].This explains why historically, the µ dependence was interpreted as an ambiguity.Now, we see that requiring reparametrization invariance actually points to a preferred basis of operators for H NR .
The fact that the axioelectric operator is screened can be demonstrated in an alternative way, shedding a different light on the mechanism at play behind the Schiff transformation.Let us assume for now that [γ • p, a] = γ • ∇a = 0, and define In the interaction picture, ψ I (t) = exp(iH NR 0 t) |ψ(t)⟩, the time-evolution of ψ I (t) can be encoded into the evolution operator with T the time-ordered product, and such that ψ I (t) = U (t, t 0 ) ψ I (t 0 ) .The interaction picture perturbation is the same as the Schrödinger one at leading order since [γ 0 m, V(t)] = 0. Now, we see that whenever V(t) = ∂ t X (t), the evolution operator collapses to the universal U (t, t 0 ) = exp(i(X (t) − X (t 0 ))), and drops out when the initial and final times are set by the experimental conditions.Thus, perturbations that are total time derivatives do not change energy levels 3 .Note also that U (t, t 0 ) is precisely the Schiff transformation done in Eq. ( 50) to get rid of the perturbation in the first place.It corresponds to ψ I (t) → ψ ′I (t) = exp(iX (t)) ψ I (t) , with then i∂ t ψ ′I (t) = 0 since V ′I (t) = V(t) − ∂ t X (t) = 0.In this picture (in the quantum mechanical sense), since X (t) commutes with H NR 0 up to terms of O(m −3 ), ψ ′I (t) stays fixed to some linear combination of eigenstates of the free Hamiltonian H NR 0 .In this sense, performing the Schiff transformation to get rid of V(t) produces a non-relativistic Hamiltonian that better reflects the physics of the system.That is the same idea as the original Schiff transformation for EDMs: H NR EM in Eq. ( 23) better reflects the energy level of the system than that in Eq. ( 18).

For charged fermions in an external EM field
The situation described in the previous section changes in a crucial way in the presence of minimally coupled electromagnetic fields.To show this, let us repeat all the steps of the previous section, but starting from Note that [γ • P, a] = [γ • p, a] = −iγ • ∇a since a is electrically neutral.This Hamiltonian can be block-diagonalized by plugging in Eq. ( 15).This produces where H NR EM is the electromagnetic Hamiltonian, Eq. ( 16), and H 3 is obtained from the neutral one in Eq. ( 46) by replacing γ • p → γ • P and p 2 → P 2 + eγ 0 γ 5 γ • B (which is nothing but (γ • p) 2 → (γ • P) 2 ).Compared to the neutral case, the only unexpected new addition is the EDM coupling S a γ 5 γ • E = 2µaγ 5 γ • E + ....Because it does not arise starting from the axion derivative interaction, it does not appear in the literature (though it is present in Ref. [45]).
As in the free case, to get a better handle on the physical couplings, let us perform the sequence of Schiff transformations: After this, the Hamiltonian becomes with Let us now expand S a and C a and keep only terms up to O(m −3 ).This calculation is simpler than it seems because most of the algebra done in the neutral case relied on the use of commutator and anticommutator identities, see Eq. ( 58), which remain essentially valid.One only has to pay attention to the extra γ Putting all together, the µ dependence again cancels out precisely, H NR (µ) = H NR , and H 3 greatly simplifies to only a few operators: where Apart from the new EDM coupling, this expression is identical to the neutral case, but for γ • p → γ • P.This is not our final form for the Hamiltonian.Because of their importance, we think it is crucial to keep track of the redundancies when they involve operators of the same order.So, let us reintroduce two free parameters explicitly and perform a final unitary transformation Then, we obtain: with the understanding that the choice of α and β is totally free.Remember that the axion scale has to be put back in these operators by writing a → ga with g = m/Λ and Λ the PQ breaking scale.To be more explicit, the leading operators read where we used γ 5 γ = −γ 0 ⊗ σ to put the operator in the standard form.

On the axioelectric -axionic EDM equivalence
At this stage, we have two operators at O(1/m 2 ) whose relative weight can be freely tuned, but whose overall impact must be identitical.In other words, starting from H NR (α, β), α and β must drop out of physical observables.Clearly, this means that the covariant axioelectric operator and the axionic EDM operators must be equivalent: At first sight, these operators appear to encode different physics.One depends on ȧ but not on E, while the other depends on E but not on ȧ.Yet, as we now discuss, there are several ways to interpret this equivalence, and to understand that at the level of observables, both operators always end up being strictly equivalent.
Table 1: Schematic representation of the equivalence of Eq. ( 74) at the level of observables.In the top line, time-dependent perturbation theory is understood, while for the bottom line, the operators are understood to be sandwiched between bound fermionic states.
• Double screening: The interplay between the γ 5 {γ • P, ȧ} and aγ 5 γ • E operators should have been expected.We know from the previous section that in the absence of electromagnetic fields, γ 5 {γ • P, ȧ} → γ 5 {γ • p, ȧ} can be eliminated.And, Schiff's theorem is telling us that if the axion field is constant, aγ 5 γ • E becomes a fixed EDM coupling that can be rotated away.So, we see that both a time-varying axion field and minimal couplings to the external electromagnetic fields are required to get a physical effect.Each form of the operator makes one of these screening manifest, but only their equivalence embodies their true impact on the physics.
• Duality in the Dirac equation: Before the Schiff transformation, the derivative and polar representations do not match and rather produce, from Eq. ( 66), Derivative : Polar : Yet, upon the equivalence of Eq. ( 74), one can choose to put the whole O(1/m 2 ) part of the Hamiltonian into the form of an EDM coupling −eaγ 5 γ • E/m 2 , which then has a very straightforward interpretation.For a charged particle, it is well-known that the Dirac equation predicts a magnetic moment g = 2 via the Zeeman term, The axion coupling to fermion can be viewed as a pseudoscalar mass term, dual to the scalar mass term.As a result, the Dirac equation then predicts an electric moment, since duality interchanges B and E. In this sense, the prediction d = ea/m 2 for the axionic EDM is the exact analogue of g = 2 for the magnetic moment.It represents an inescapable consequence of the Dirac equation whenever the charged fermion has a pseudoscalar coupling to the axion4 .The observability of this oscillating EDM is another question though, because one must fight the various screening effects, and will be discussed in Sec. 4.
Basically, this contribution is kinematically suppressed, and encoded into axion wind operators of O(1/m 3 ).Then, at leading order, what the Schiff transformation Eq. ( 67a) is telling us is that a coupling ȧσ • A is equivalent 5 to a coupling aσ • Ȧ, that is, aσ • E, exactly like the transformation Eq. ( 67c) is telling us that ȧ2 encodes the same physics as aä.These pairs of operators must give the same result when acting on fermion wavefunctions within physical observables (see Table 1).This is clearly in accordance with the interaction picture evolution of Eq. ( 61), where a perturbation like aσ • E or aä is to be integrated over time.This also shows how the original Schiff screening comes back if a(t) becomes constant, as the time-integral of aσ • E = ∂ t (aσ • A) then sums up to an unobservable constant rephasing of the fermion wavefunction.
• The axioelectric effects as EDM-induced: If one encodes entirely the O(1/m 2 ) terms into an EDM operator, the usual axioelectric effect is nevertheless still there.Because the two forms in Eq. ( 74) are equivalent, the same matrix elements as in Ref. [23] must be recovered.
The equivalence in the ϕ = 0 gauge was discussed above, so let us now concentrate instead on that with A = 0, so that E = −∇ϕ.In this case, the covariant axioelectric operator reduces as γ 5 {γ • P, ȧ} = γ 5 {γ • p, ȧ}.Though identical in form with the neutral fermion axioelectric operator, it cannot be rotated away when ϕ ̸ = 0 (there is an extra term in Eq. ( 52) from [iS 1 , ϕ] ̸ = 0 with S 1 in Eq. ( 50)).As explicitly calculated in Ref. [23], the γ 5 {γ • p, ȧ} can then induce transitions between energy levels for an electron bound in the potential ϕ.Now, starting instead from the EDM operator, we can write for an electrostatically bound electron 6 , eaσ • E = −eaσ • ∇ϕ = −iaσ • [p, H NR ] + ....With this, the transition matrix element of Ref. [23] is trivially recovered as ⟩ and integrating by part over time (see Table 1).Note that these mathematical steps essentially undo the Schiff transformation of Eq. (67a).This shows that whatever the operator, the same matrix element for the axioelectric effect is obtained, as it should since physics must not depend on the representation.Yet, we think it sheds new light to interpret the axioelectric effects rather as a manifestation of an axion-induced EDM, especially in view of the other points discussed previously.

For neutral fermions having an EDM interaction
The final application is the non-relativistic limit of the Hamiltonian for a neutral state coupled to the axion, but in the presence of both the magnetic and electric dipole operators.Those are not invariant under the PQ symmetry, so one has to decide how they should be introduced.We consider that they arise in the same way as the mass term, through the PQ symmetry breaking.
The fermion field is neutral under the PQ symmetry only in the derivative representation, so those effective operators can be added to Eq. ( 37) as Then, if we use again Friar's trick, Eq. ( 42), the Lagrangian interpolating between derivative and polar representations is where and S a = m sin(2µa/m), C a = m 2 (cos(2µa/m) − 1) as before.The fact that the dipole operators end up proportional to the axion field in the exponential representation is similar as in the SM, where they necessarily involve the Higgs boson field [67].
The Hamiltonian in the non-relativistic limit can be obtained by plugging the odd and even elements in Eq. ( 15).After some algebra, and noting that δµ = δ µ + O(1/m 2 ) and d = d + O(1/m 2 ), we arrive at with the same H 3 as before, Eq. ( 46).We have used Notice that the same combination of S a and C a as in the last line of Eq. ( 45) has already been dropped.Similarly, S 2 a + 2C a = −4m 2 sin 4 (µa/m) ∼ O(1/m 2 ) can be discarded.Again, we observe that the infinite towers of interactions in the polar representation, the exp(2ia/m)σ µν F µν and exp(2ia/m)σ µν Fµν terms in Eq. ( 80) for µ = 1, are automatically truncated when expanded in the non-relativistic limit.This fact would have been totally missed if we had truncated the series already in Eq. ( 80).There is another interesting aspect of this truncation.Setting µ = 1 in Eq. ( 80), a direct coupling of the axion to the electric dipole operator a ψσ µν γ 5 ψF µν is present in the polar representation, but not in the derivative one, and with a coefficient proportional to the magnetic moment of ψ.We now see that this coupling disappears in the non-relativistic limit, making both representations compatible.This information will play an important role in analyzing nucleon EDMs in Sec. 4.
At this point, we start to perform some Schiff transformations, specifically, S 1 as given in Eq. ( 50), S 2 in Eq. ( 53), S 3 in Eq. ( 54), and finally, to remove the operator involving δ µ Ė + d Ḃ.The S 2 and S 3 transformations reorganize the terms in ä occurring in H 3 , exactly as in Eq. ( 55).For S 1 , an additional term appears (compare with Eq. ( 52)) This new term combines with the {S a , [γ • p, γ • (δ µ B − dE)]} operator of Eq. ( 84) to make it µ independent.The other terms combine with those in H 3 as in Sec.3.1, and the final Hamiltonian no longer depends on µ at all in the non-relativistic limit: One can recognize in the first and fourth lines the terms of H NR EM in the e → 0 limit, Eq. ( 34), the terms in the second and third as those of the neutral fermion, Eq. ( 59), so the only new feature is the operator in the last line.It encodes higher order effects induced by the Schiff transformation of the ȧ term, and can be worked out using where we have set ∇ • B = 0 and ∇ × E = 0.These couplings are rather similar to the Darwin and spin-orbit couplings, and for a neutral fermion, should not be directly accessible.Further, whenever d is already induced by the axionic background, these couplings represent a negligible second order effect.So, the important conclusion of this calculation is that for a neutral state, there is no coupling to the time-derivative of the axion background, but the leading non-axionic EDM term is physical.

Axionic EDM observability and estimates
Whether in its axioelectric or axionic EDM form, the O(m −2 ) operators lead to oscillating EDM for charged states, i.e., to charged leptons and quarks, and thereby presumably to nucleons.Naively, if the analogy with the magnetic moments is valid, the expected size of these oscillating EDMs should be larger than those due to the axion-gluon or axion-photon local anomalous interactions, aG µν Gµν or aF µν F µν , respectively.In some sense, the O(m −2 ) operators are the equivalent of the Dirac prediction g = 2 (see Eqs. (77,78)), while the loop-level anomalous contributions play the role of the anomalous magnetic moment.At the same time, there are also reasons to think that this analogy does not hold because the O(m −2 ) operators and the anomalous contributions to a given particle EDM appear different in their scaling with the axion mass, in their connection with Schiff's screening, and in the case of the nucleons, in their hadronization.So, to confirm the naive expectation, it is necessary to delve into the details of how these operators translate into observables, from which oscillating EDMs could in principle be accessed.We will not do a systematic analysis of all the possible experimental settings, and to avoid complications arising from nuclear and atomic effects, we will focuss on the simplest system consisting of a single particle precessing in external EM fields, first in the leptonic case, and then in the more intricate nucleon case.This treatment will prove sufficient to establish once more the equivalence between the axioelectric and axionic EDM operators, and to show and characterize the situations in which they do lead to much larger EDMs than expected on the basis of the anomalous axion couplings alone.

Leptons
Imagine a charged lepton, not bound in an atom, travelling in a region where ∇a is negligible and only some electric fields are present.Then, its spin precession is dictated by the Hamiltonian in Eq. ( 72) as, to leading order in 1/m, In writing this equation, we implicitly use the fact that the spin operator is the same for all α (but for the gauge variance due to P, which cancels only when acting on the fermion wavefunction).This makes sense because, compared to the discussion in Sec.2.3.2, the Schiff transformation of Eq. (67a) relates two equivalent forms of the axion coupling.It does not affect the mass term, which anchors the rest-frame in which S is defined.In other words, at the level of the Lagrangian, the original Schiff transformation of Eq. ( 20) is related to the purely chiral rotation exp(iαγ 5 ), as was made clear in Eqs.(77,78), while that for the axion coupling, Eq. (67a), is rather related to the Goldstone boson reparametrization exp(iαγ 5 a/m).
Reparametrization invariance.Let us first prove that observables do not depend on α.In the ϕ = 0 gauge, the γ • p term does not contribute and E = − Ȧ, so the equation for S reduces to Since this equation is of the ḟ = f × g type, for which a solution involves the time integral of g, the two operators give the same contribution and α drops out (under appropriate boundary and gauge conditions).We recover the equivalence of the two representations for the O(m −2 ) axion couplings, as depicted on the top line of Table 1.
If we take instead a gauge where A = 0, the γ • p produces the required ∇ϕ term of E, as depicted in the bottom line of Table 1.Let us check this explicitly, assuming ∇a = 0. We start with the EDM form of the operator in the equation of motion of S, and first replace eE the Hamiltonian with a = 0: Acting with H NR 0 on the external spinors gives the energy difference, which comes entirely from the difference in electric potential felt by the initial and final states.This must match the energy brought in by the axion, i.e. m a , so that aγ 5 γ • [p, H NR 0 ] collapses to ȧγ 5 γ • p, as it should 7 .In full details, the correspondance is verified by writing explicitly the external spinors.According to the Ehrenfest theorem, using the equation of motion of the external states and integrating by part over time, i.e., imposing the conservation of energy.We do not include ⟨ψ| Ṡ |ψ⟩ in the evolution of ⟨ψ| S |ψ⟩ and used [S, H NR 0 ] = i Ṡ = 0 because as an operator, S does not vary in time.This is consistent since in this picture, H NR 0 = γ 0 (m + ...) and S = σ/2 = −γ 0 γ 5 γ/2, so clearly [S, H NR 0 ] = 0. We thus recover the equivalence between the axionic EDM and axioelectric operators of Eq. (74).All in all, the situation is totally analogous to that for the axioelectric effect discussed in Sec.3.2.1:observables are gauge invariant and independent of the choice of operators in H NR (α, β).Numerical estimates.Altogether, the α parameters drops out from the equation of motion of S which can be written as Naively, this could represent a significant effect.For an electron, taking the coherent classical axion background a(t) = a 0 cos(m a t) with m a a 0 = √ 2ρ DM and ρ DM = 0.4 GeV/cm independently of m a and Λ [68,69].Alternatively, this same estimate can be expressed in terms of the axion-electron coupling (which corresponds to taking a → g ae a in Eq. ( 72)), The electron EDM thus appears very promising to set competitive bounds on the axion (or ALP) couplings, especially when compared to the current limit d exp e < 1.1 × 10 −29 e cm [70] for a fixed d e σ • E coupling.Yet, the two cannot be immediately related, not least because d exp e is extracted from atoms, i.e., bound electrons, but also because even for freely precessing leptons, there are conditions hidden in Eq. ( 95).Also, trivially, the axionic EDM integrates to zero if the observation time T runs over several oscillation cycles, and the electric field |E| = E 0 is constant.
An immediate question looking at Eq. ( 95) is what happens if the axion is sufficiently slowly varying compared to T , i.e., when T is small enough compared to 1/m a .In effect, the axion field is constant, and it may seem Ṡ ∼ aS × E survives as m a → 0 since it is linear in a.This is not true though, because if the axion field is constant, then S is no longer the right spin operator.With a(t) = a 0 , the mass term becomes complex (see Eq. ( 39)), and a chiral rotation of the wavefunction becomes necessary to identify the true spin operator.This additional change of basis, required only if the axion field is constant over the observation time (or other relevant time scale), is what distinguish the situation with α = 0 from α = 1 in Eq. (91).More generally, if we write a(t) = a 0 + ȧt + ..., then the a 0 term disappears when the fermion mass is made real.When T ≲ 1/m a , the observable change in the spin orientation ends up linear in m a , with a |E| ≈ E 0 a 0 m a T .This makes the precession inobservable for very small axion masses, reproduced the expected decoupling of the axion in the m a → 0 limit, but makes the exploitation of d exp e for constraining axion interactions impossible.
More promising are situations in which the electric field is oscillating in time at some frequency ω ≈ m a .Then, no change of basis is needed, a |E| = E 0 a 0 cos(m a t) cos(ωt) is not linear in m a , and it integrates to a non-zero value over some long enough observation time.Clearly, one can no longer take the m a → 0 limit here, it is ill-defined for this situation.Also, the matching with the axioelectric form is trivial since to |E| = E 0 cos(ωt) corresponds |A| = (E 0 /ω) sin ωt, so ȧ |A| matches a |E|, up to fixed boundary terms.This shows how oscillating electric fields really takes the full advantage of the fact that the derivative of the axion field, ȧ, is coupled to the vector potential, and not to the electric field8 .Thus, provided the axion mass is not too small so that T covers more than a fraction of an oscillation, it is in principle possible to access directly to the axion-induced EDM, as predicted by Eq. ( 96).This could be particularly interesting with electric microwaves, covering the µeV range of axion masses favored by the misalignment mechanism.Note also that this reasoning remains valid in a CASPER-like situation [16] in which the spin of the charged lepton is precessing at a Larmor frequency ω L ≈ m a in a magnetic field.In that case, a constant electric field is seen in the rotating frame as oscillating at that ω L frequency and Eq.(96) applies.
Let us stress though that Eq. ( 96) does not apply to CASPER itself as currently designed, since we are dealing with free charged leptons here.Further, we do expect very significant suppressions for charged fermions bound into a neutral atomic system, as will be shown in the next section for the neutron.Further work is needed to cover these systems, to see whether some sensitivity can be retained for some range of axion masses.In that respect, the equivalence between the axionic EDM and axioelectric forms of the operator may turn out to be useful.We already know that they both have precisely the same capability to kick bound electrons out when the axion brings enough energy, and we have seen above that they act in the same way on free lepton spins, so there is no reason to think they could act differently on bound electrons.As a tool, this equivalence could thus help in obtaining realistic numerical estimates for atomic systems.This is left for future work.

Nucleons
Let us now turn to the nucleons, for which the situation is much more complicated because of hadronic effects, and because the quarks are not expected to be non-relativistic inside a nucleon.We concentrate on the connection between the quark and nucleon levels here, with an idealized nucleon precession experiment in mind, and leave the discussions about the nuclear or atomic levels to future work.To proceed, let us characterize the various contributions to the nucleon EDMs in terms of effective Lagrangian couplings.
1-Quark constant EDMs.The simplest mechanism to generate a nucleon EDM occurs when quarks develop constant EDMs, like in the presence of some new CP violating sources.The quark electric moment operators then translate naturally into the corresponding nucleon operators with q = u, d.Given current lattice estimates [59], hadronization appears essentially transparent on these local operators, and the naive SU (6) estimate in Eq. ( 36) is actually quite good.In case the short-distance quark EDMs scale as 1/m q , analogy with the magnetic moment would suggest a similarly transparent hadronization, but with running quark masses replaced by constituent masses.The important information is that even if the neutron is neutral, it is not "EDM neutral", because the EDM interaction involves a combination of spins and electric charges.As a result, even very soft photons, insensitive to the quark structure, do interact with the neutron via its spin.
2-Gluonic contributions.The fundamental axion interactions with quarks or gluons do not involve the photon field at leading order.Yet, the quarks being electrically charged, nonlocal processes at the partonic level can induce local EDM operators at the nucleon level.The most well-known such non-local EDM contribution comes from the θ term of QCD.Indeed, in the presence of an axion field, one expects a coupling The constant θ term is cancelled by the axion field falling to its true minimum, but this leaves a aG µν Gµν coupling.In the presence of a dark matter axion background, a(t, x) = a 0 cos(m a t− k • x + ϕ), and from Eq. ( 35), this term then induces an EDM for the nucleons [68,69].Using the matrix element estimates quoted in Ref. [57]: (we use d N (t) to denote the total nucleon EDM, and d N the coefficient of the operator in Eq. ( 100)).Note though that recent lattice estimates reduce this matrix element by about a factor of two [71,72].Anyway, this prediction has motivated dedicated experimental searches [16,73], with specific strategies designed to tackle the oscillatory nature of the EDM.
At this stage, we should point out though that strictly speaking, the matrix elements was extracted for a constant θ, by extrapolating the form-factor for N → N γ(q) to q 2 → 0. The idea is that provided the axion background is not varying too quickly, QCD has time to account for the presence of the action field as a kind of effective axionic θ term.Yet, notice that in the general case, the axion field is also injecting some energy, albeit a small amount, and d N in Eq. ( 100) is actually a form-factor that depends on both the photon and axion momenta.When applying this estimate to the axion EDM coupling, one implicitly makes the assumption that the axion is very soft and that the limit m a → 0 is smooth.We will come back to this point below.
3-Ward Identity.The pseudoscalar and/or derivative couplings of the axion to the quarks also contribute to the nucleon operator in Eq. (100) via similar non-local processes (see Fig. 1).To estimate their relative size compared to the gluonic contribution, a crucial piece of information comes from the Ward identity of the anomalous PQ symmetry (we assume for now that the axion couples only to a single quark ψ q ): This means that the Lagrangian interpolating between the derivative and polar representations should read where one can recognize Eq. ( 39) when α = 0, and Eq. ( 37) when α = 1, plus the anomalous term.Thus, the gluon-induced and quark-induced contributions cannot truly be disentangled and must be treated together.This fact is actually often used to estimate theoretically the nucleon EDM for constant θ, as it can be easier to deal with a phase for the quark masses than with the anomalous G G coupling [74,75].Note, though, that if the aG G coupling also receives contributions from other heavy states, either SM quarks or new heavy fermions like in the KSVZ scenario, then the quark Lagrangian should rather read where we have put back the coupling g q and g g to distinguish axion-quark and axion-gluon couplings.Notice that in this case, it is always possible to chose α (or β) such that one of the coupling disappears, i.e., without the axial, the pseudoscalar, or the anomalous coupling: Also, notice that in the quark non-relativistic limit, g g never contributes to the quark axion wind or the quark axionic EDM operator.
4-Sutherland-Veltman theorem.What the Ward identity Eq. ( 102) shows is that in the soft limit, i.e., ∂ µ a → 0, the pseudoscalar axion-quark coupling is strictly equivalent to the anomalous axion-gluon coupling, and thus that Actually, this is essentially the Sutherland-Veltman theorem [76,77], well-known in the context of π 0 → γγ.Here, it proves that the non-local contributions of a ψq γ 5 ψ q to the nucleon EDM must reproduce that of aG µν Gµν quoted in Eq. (101) in the soft limit.This serves as a baseline, and the question now is whether the EDM can be enhanced compared to Eq. (101) even slightly away from that limit.Note, for completeness, that the Ward identity also relates the matrix elements of the axial current and the anomalous term in the chiral limit.Indeed, for a massless quark, the axion decouples entirely from Eq. ( 103), as can be seen taking m q = 0 and α = 0. Thus, the matrix elements of ψq γ µ γ 5 ψ q ∂ µ a and aG µν Gµν must match when m q → 0 since the two must cancel each other [61].
5-On the m a → 0 limit.The soft limit, ∂ µ a → 0, and the m a → 0 limit are not entirely equivalent.Naively, if a(t, x) = a 0 cos(m a t − k • x + ϕ) for some constant ϕ, ∂ µ a → 0 is equivalent to ȧ = 0 if k is negligible, and ȧ(t) = 0 can be attained for all time only with m a = 0. Now, simply setting m a = 0, the axion field becomes constant and some of its couplings to quarks and gluon in Eq. (103) survive, and so is the nucleon EDM operator in Eq. (100).At the same time, theoretically, m a → 0 requires Λ → ∞ since m a comes from the gluon coupling aG µν Gµν .But if we send Λ to infinity, the axion entirely decouples, there will be no axionic EDM at all, and the strong CP puzzle is back.Actually, there is probably a threshold for m a below which the axion would take too much time to realign to compensate for some preexisting θ term, in which case all the axion contributions would be overshadowed by the large constant EDM due to θ.This goes beyond what is discussed here, and we still assume that term is absent.Yet, if we expand the coherent background axion field as a(x µ ) = a 0 + x µ ∂ µ a + ..., how to treat the constant term a 0 needs caution.In effect, QCD with such an external field looks very much like QCD with a non-zero θ term, which is CP-violating and quite different from QCD with an axion and θ = 0, which is CPconserving (see e.g.Ref. [78] for a comparison).In particular, the a 0 term generates complex phases for the quark condensate, and changes how to treat the chiral symmetry breaking terms.Specifically, when a = a 0 is constant, Dashen theorem must be called in, and after the necessary realignment of the chiral vacuum [74], This drastically changes the character of the axion-quark coupling, and suppresses it significantly (it now vanishes if any of the quark masses vanishes, as it should).This isospin singlet quark current accounts for the isospin singlet anomalous term a 0 G µν Gµν .This explains how Sutherland-Veltman theorem sets in: the L q,g (α = 0) term contains both a suppressed term collapsing to the anomalous one, and a term matching the axial interaction (i.e., proportional to ∂ µ a), so that any observable calculated from L q,g (α = 0) or L q,g (α = 1) are equal.
6-Axion couplings to nucleons.The three axion couplings in L q,g (α) of Eq. ( 103) generate the nucleon EDM operator, but also simpler axion couplings to nucleons.With them, the nucleon effective Lagrangian becomes: for some prefactor g N a priori of O(1), and D µ = ∂ µ − ieQ N A µ with Q N the nucleon electric charge.From the previous points, all the L q,g (α) couplings appear to contribute to all the L N (α ′ ) couplings.For instance, the axion-gluon coupling contributes to both d N and g N , since aG µν Gµν does not need photons to generate a CP-violating coupling.Similarly, the axion-quark couplings in L q,g (α) naturally induce axion-nucleon couplings, but also the local anomalous term d N , as is obvious starting from L q,g (α = 0) or invoking Sutherland-Veltman theorem, Eq. ( 106).Yet, there are subtelties at play here.First, even if the the two axion couplings to nucleons are necessarily present and related under the reparametrization ψ N → exp(ig N αγ 5 a/Λ)ψ N , α ′ is not necessarily equal to α.We neglect the QED anomaly here, so the Ward identity underpinning the Goldstone-boson reparametrization invariance at the nucleon level is simply the classical one (and α ′ has to cancel from observables).Second, for a constant axion background a(t, x) = a 0 , the axion couplings to nucleons vanish.This is clear for α ′ = 1, while it requires a chiral rotation ψ N → exp(iαg N γ 5 a 0 /Λ)ψ N to make the mass term real for α ′ = 0 (this can also be understood as a nucleon-level Sutherland-Veltman theorem: the pseudoscalar coupling is equivalent to the axial one, which vanishes in the ∂ µ a → 0 limit).Third, the d N coupling does not represent the whole axionic EDM of the nucleon.From an effective theory point of view, d N only represents the short-distance contribution, to which tree-level (and loop-level once pions and other light mesons are included) contributions from the leading couplings have to be added, see Figs. 1 and 2.
Figure 2: Partonic contributions to the nucleon-axion and nucleon-photon interaction generate a non-local, long-distance EDM effect for the proton, at the hadronic level and in the non-relativistic limit.The axion-gluon and axion-quark contributions are related by a Ward identity, but their combination is imposed to ensure that the axion-nucleon interactions decouple entirely in the soft limit, ∂ µ a → 0. This is justified since in that limit, a non-decoupling axion-nucleon coupling corresponds to a complex mass term, so it has to be removed by a chiral rotation of the nucleon field.
7-Nucleon magnetic moment.Deciding like in Sec.3.3 to add the magnetic dipole operator to L N (α ′ = 1), the nucleon Lagrangian becomes (see Eq. ( 81)) up to dipole terms of O(a 2 ).The magnetic dipole operator accounts for the proton and neutron magnetic moments 2 + µ p ≈ 2.8 and µ n ≈ −1.9, respectively.This Lagrangian is invariant under the Goldstone boson reparametrizations ψ N → exp(ig N α ′ γ 5 a/Λ)ψ N , again up to terms quadratic in the axion field.In the non-relativistic limit, a single axionic EDM operator arises at leading order: The contribution from µ N drops out, independently of α ′ (as was already apparent for the neutron in Eq. ( 89)).It is interesting to understand the mechanism of this cancellation (see Ref. [79]): the µ N -dependent shift of d N is compensated by the long-distance contributions arising from one-nucleon reducible diagrams with the photon emitted via the µ N ψN σ µν F µν ψ N vertex.Now, this cancellation crucially relies on the assumption that one should introduce the magnetic dipole operator to L N (α ′ = 1).In Sec.3.3, this was justified by the fact that quarks are PQ neutral in the derivative representation.Nucleons, on the contrary, do not have definite PQ charges if the PQ-breaking aG G coupling contributes to the a N N vertex.
To circumvent this, we require that only the local term d N should be present in the ∂ µ a → 0 limit.Indeed, the reparametrization then becomes a chiral rotation, i.e., a change of basis.What is to be called the magnetic moment and the EDM have to be defined in the basis in which the nucleon mass is real (see e.g.Ref. [80] for a detailed discussion).With the assumption that this limit is smooth despite the fact that a change of basis for ψ N is implied, this restores a definite PQ charge for the nucleons.The representation in Eq. ( 109) confines the contributions of the aGG coupling to the local term d N , together with some axion-quark contributions given Eqs.( 105) and (106), leaving the quark couplings to induce g N .Thus, we expect d N (g q , g g ), but g N (g q ), with the g N (g q ) contribution cancelling out when ∂ µ a = 0. Note, finally, that L N (α ′ = 1) corresponds to the usual form employed in the literature, see e.g.Ref. [69].

7-Proton EDM:
The important property of the first two operators of L N (α ′ ) is that they combine with the electromagnetic coupling to produce non-relativistic EDM operators for the proton, but not for the neutron.The phenomenology of a precessing proton is thus very similar to that discussed for charged leptons.In the non-relativistic limit for the proton, using L N (α ′ = 1), the axion EDM couplings are The covariant axioelectric operator cannot be rotated away, and does induce spin precession, see Eq. ( 94).We choose to write H NR starting from L N (α ′ = 1) instead of as in the equivalent form of Eq. ( 110) to emphasize the different nature of these two contributions.First, the g p is explicitly vanishing if ∂ µ a = 0, but not the d p term.Yet, both sum up to an EDM-like precession, and even more, both operators becomes identical if |E| = E 0 sin(ωt) since then ȧγ 5 γ • P ⊃ ȧγ 5 γ • A = aγ 5 γ • E provided m a = ω ̸ = 0. Second, these two operators encode different physics: the d p contribution is local already at the hadronic level, but the axioelectric contribution has an intrinsically non-local origin, see Fig. 2, and becomes local in the nonrelativistic limit only (as said earlier, we do not consider nuclear systems here).From the discussion in the lepton case, we do expect that the axion-induced proton EDM to increase to in the "resonant" situation in which the EM field matches the axion frequency, and assuming g p ∼ O(1).This represent an enhancement of the long-distance contribution by about two orders of magnitude compared to the local contribution tuned by d p , Eq. (101).Beware though that, obviously, the same provisions about the implicit observation time constraints as in the lepton case do apply, since the g p contribution does decouple if ∂ µ a = 0, as is manifest in Eq. (111).Note, finally, that there is some model dependence in comparing Eq. (101) to Eq. (112).For instance, hadrophobic scenarios can be designed in which the axion couplings to the u and d quarks conspire to suppress g p (see e.g.Ref. [81] and references there).Barring these possibilities though, Eq. ( 112) may represent our best window into the axion-light quark couplings, even compared to axion wind operator that are suppressed by the local galactic axion speed [14].

8-Neutron EDM:
The purely long-distance enhancement mechanism at play for the proton is not available for the neutron since it is neutral, see Fig. 2 (as said earlier, this crucially rely on how the neutron magnetic dipole operator is introduced though).Instead, with only the a ψn σ µν γ 5 F µν ψ n coupling, there will be an enhancement if the non-local quark-level matrix elements can be significant away from the m a = 0 limit, so that the enhancement identified at longdistance somehow spills over at short-distance.If that is the case, this violation would show up in d N , which should be understood to be a form-factor: d N (g q , g g ) = d N (g q , g g ; q 2 γ , q 2 a , q a • q γ ) .
While we know that d N (g q , g g ; q 2 γ , q 2 a , q a •q γ ) → Eq. ( 101) when ∂ µ a = 0, the behavior reaching that limit may not be that smooth if Eq. ( 113) does not go to zero sufficiently fast as ∂ µ a → 0. Let us imagine that the proton and the neutron are simply collections of loosely bound nonrelativistic constituent quarks.Then, the long-distance hadronic mechanism at play for the proton would have a direct counterpart as a non-local constituent quark mechanism (e.g. from the third diagram in Fig. 1).Both the proton and the neutron EDM would then be expected to reach Eq. (112) in the presence of "resonant" EM fields since, as explained in point 1 above, the neutron is not neutral for spin-dependent electric interactions.In practice, in this picture, one way to understand Eq. ( 111) would be from a term in d N scaling like q 2 a /q a • q γ = m a /ω, vanishing in the m a → 0 limit, but of O(1) when ω ≈ m a and m a is not too small.Of course, this consituent quark picture is not particularly realistic, but in our opinion, it nevertheless suggests that some level of enhancement of the neutron EDM is possible.Indeed, the real world situation should lie somewhere in between no enhancement, as expected looking at Fig. 2 with the neutron not interacting with photons, to a significant enhancement thanks to residual interplays between the axion and photon couplings to the quarks inside the neutron.Obviously, to get a definitive answer from first principle is complicated and probably requires detailed lattice simulations starting from the general Lagrangian of Eq. ( 103), away from the m 2 a = 0 and q 2 γ = 0 limit.
To close this section, we stress once more that the above discussion does not immediately apply to nuclear or atomic probes of the axion-induced proton and neutron EDMs (assuming the axion-electron coupling is absent).The non-relativistic limit appears crucial to collapse the axion couplings to an EDM-like operator for the nucleons, which can then be enhanced with suitable EM fields.Further, estimating how an oscillatory external electric field can penetrate the nuclear, atomic, and/or even the molecular system, accounting in addition for the presence of resonances, and estimating the resulting observable EDM it would induce is beyond the scope of the present work [82][83][84][85].

Summary
In this paper, the non-relativistic description of the axion interactions with fermions was systematically analyzed.We relied on rather old and well-established techniques like the Foldy-Wouthuysen transformation [46], the unitary transformations of Ref. [40], and Schiff theorem [47].Yet, as these techniques had not been fully combined and supplemented by the reparametrization invariance for the axion field, to our knowledge, none of the final non-relativistic expansions for the Hamiltonian presented here were derived before.Our results can be summarized in three points: • For a neutral fermion, we demonstrated by adapting Schiff theorem that the axioelectric operator γ 5 {γ • p, ȧ} is totally screened.As shown in the final Hamiltonian for this scenario, Eq. ( 59), there are only axion wind operators up to O(1/m 3 ), except for a very suppressed ȧ2 coupling.Since there should be no finite-size effects, and because O(1/m 3 ) relativistic corrections are of a different nature, this screening should even hold to a much higher level than the usual Schiff screening of charged fermion EDMs.Phenomenologically, this scenario is not very relevant since normal matter is essentially made of charged particles, but it provides the basis to understand the result in the charged case.
• Specifically, for a charged fermion, the final Hamiltonian is in Eq. ( 72).The covariant axioelectric operator γ 5 {γ • P, ȧ} is found equivalent to an axion-induced EDM operator aσ • E, see Eq. ( 74) and Table 1.Both operators encode the same physics, but the latter makes it manifest that this coupling disappears in the absence of EM fields, or for a neutral fermion.Phenomenologically, the usual axioelectric effect is recovered whatever the chosen form of the operator, both having the same matrix elements for observables.Besides the axioelectric effect, these operators can also induce EDMs for all charged particles.The important points are first that these EDM operators are, in some sense, tree-level.They are directly predicted by the Dirac equation itself for all charged fermions, in a way totally analogous to the magnetic moment factor of 2. Secondly, these EDMs are not constant in time, and cannot be screened since Schiff transformation would simply change the relative weight of γ 5 {γ • P, ȧ} and aσ • E, something irrelevant since they lead to the same observables.Thus, though specific search strategies have to be designed to tackle the oscillatory nature of these EDMs as well as their decoupling in the m a → 0 limit, their relatively large sizes, especially for the electron, makes them particularly promising.
• Finally, concerning the proton and the neutron, the final Hamiltonians are in Eq. ( 72) and Eq. ( 89).The main issue here is whether the axion-induced quark EDMs, which are intrinsically non-relativistic, can manifest themselves at the hadronic level.We find that this is the case for the proton, whose axion-induced EDM is significantly enhanced by long-distance effects compared to current estimates based solely on a local EDM operator induced by the axion-gluon coupling (see Fig. 2).For the neutron, if taken as point-like in a first approximation, the axion-induced EDM coupling coming from the axion-quark couplings vanishes exactly since the purely long-distance hadronic contribution is absent.Beyond leading order, some effects are likely as the neutron is not transparent to quark EDM interactions, but further work is needed to estimate these finite-size effects and establish whether they can compete with the axionic EDM coming from the axion-gluon coupling.
All these results clarify the construction of non-relativistic expansions in the presence of Goldstone bosons.Yet, to conclude, we would like to stress again that this formalism, in itself, has some limitations.For instance, our starting point was the Dirac equation for a single fermion in the presence of external background fields, electromagnetic and axionic.In our opinion, further work is urgently needed to obtain estimates for realistic experimental settings, in particular in the atomic or nuclear contexts (or even for the neutron-antineutron system [86]).Thus, extending the formalism itself, or even grounding it within a fully relativistic quantum field theory setting, would be very welcome, not least to confirm the promising phenomenological opportunities we identified for the detection of dark matter axions.

A Foldy-Wouthuysen transformation
The Foldy-Wouthuysen (FW) procedure [46] is a systematic order by order method to blockdiagonalize the Dirac Hamiltonian via a sequence of unitary transformations.Though it is wellknown and can be found in many textbooks on relativistic quantum mechanics, for completeness, we here include a brief derivation up to O(1/m 4 ).Also, compared to the literature, we stick to the usual gamma matrices instead of the original Dirac matrices.Though inessential, this permits to immediately take advantage of computer packages, in particular FeynCalc [87].
Being perturbative, the first step is to expand the impact of a specific unitary rotation ψ → ψ ′ = e iS ψ.If i∂ t |ψ⟩ = H |ψ⟩, then i∂ t |ψ ′ ⟩ = H ′ |ψ ′ ⟩ with Using the CBH formulas, where etc, and d is a differential operator acting only on e −X , the expansion of Eq. ( 115) is Proceeding further to eliminate O ′′ with iS ′′ = O ′′ /(2m) ∼ O(1/m 3 ) does not change the diagonal term anymore since E ′′′ − E ′′ ∼ O ′′2 /(2m) is already O(1/m 5 ).So, the final Hamiltonian can be read off the result after only the S and S ′ transformations, even though a total of four FW transformations are actually necessary: where all the higher order E and Ȯ dependences occur in the chain of odd operators V i+1 ≡ [V i , E] + i Vi (this remains true at higher orders).In all the applications here, only the terms in the first line are kept.Those are obtained by the sequence of transformations ψ → e iS ′′ e iS ′ e iS ψ with iS = O/(2m), iS ′ = O ′ /(2m) and iS ′′ = O ′′ /(2m).

B Non-relativistic electromagnetic interactions
Let us start from the Hamiltonian after the Schiff transformations S 1 of Eq. ( 20) with α = md/e and S 2 of Eq. ( 24) with β = 1, keeping terms at most linear in a or d: Strictly speaking, there are not enough constraints to point us towards a specific form for the Hamiltonian.To proceed, we therefore add the requirement that the pure field-dependent terms should involve only the electromagnetic invariants E 2 − B 2 and E • B. This matches the comments made in the text about higher order operators, in particular F µν F µν or F µν F µν , that could be added to the initial Hamiltonian and would immediately contribute to these terms.With these requirements, we obtain ).This form is rather suggestive, with the 1+a factor occurring for both the Zeeman term and the γ 0 {P 2 , σ • B} operator once the B 2 term is properly tuned to force the appearance of the E 2 − B 2 invariant.Remember though that some redundancies remain in this Hamiltonian, as encoded in Eq. ( 127).
) Thus, because of the 2m term, χ is essentially determined by φ.It corresponds to a small O(v/c) component relative to the large φ component.Plugging χ ≈ σ • Pφ/2m back into the equation of φ permits to reduce the Dirac equation to a Pauli equation for φ,

with V 1 ≡
[O, E] + i Ȯ an odd operator.The leading non-block diagonal term has disappeared, and non-block diagonal terms in O ′ start now at O(1/m).Those can be removed at that order by performing a second FW transformation with iS ′ = O ′ /(2m) ∼ O(1/m 2 ).Keeping terms up to O(1/m 4 ) only, and using the above formulas, we arrive at H ′′ = γ 0 (m + O ′′ ) + E ′′ with