Sum rules for CP-violating operators of Weinberg type

We estimate the size of the hadronic matrix elements of CP-violating three-gluon and four-gluon Weinberg operators using sum-rule techniques. In the three-gluon case, we are able to reproduce the expressions given in earlier works, while the four-gluon results obtained in this article are new. Our paper therefore represents the first systematic study of contributions to the electric dipole moment of the neutron due to CP-violating dimension-six and dimension-eight operators. We provide many details on both the derivation of the sum rules as well as the analysis of the uncertainties that plague our final predictions.


Introduction
Searches for electric dipole moments (EDMs) place stringent constraints on any beyond the standard model (BSM) scenario with additional sources of CP violation (see  for reviews and recent discussions). At present the strongest limits are set by measurements of the electron spin precession in thorium monoxid [22,23], the EDM of the neutron (nEDM) [24,25] and the mercury atom [26,27]. While the thorium monoxid measurements can be interpreted as a probe of the electron EDM with small theoretical uncertainties [28,29], nucleon, nuclear and diamagnetic EDMs receive contributions from several effective operators that are plagued by theoretical uncertainties of different sizes. For instance, the EDM contributions from down and up quarks to the nEDM have been calculated with an accuracy of O(5%) using lattice QCD (LQCD) [30][31][32], while sumrule calculations [33][34][35] allow to determine the nEDM contributions from the down-quark and up-quark chromomagnetic EDMs (CEDMs) with uncertainties of O(50%). To date only estimates of the hadronic matrix element of the leading operator of Weinberg type exist. These rely on either naive dimensional analysis (NDA) [36], the vacuum insertion approximation (VIA) [37] or sum rules [38]. The resulting uncertainties are hard to quantify, but commonly said to be of O(100%). LQCD computations of the contributions of the CEDMs and the leading Weinberg operator have gained significant momentum in recent years [39][40][41][42][43][44][45][46], and considering the ongoing efforts by several LQCD groups, calculations with uncertainties similar to those of the sum-rule estimates may be achievable within the next five years [47,48]. To fully exploit the expected increase in sensitivity of future EDM searches (see for instance [49][50][51] for discussions), improved calculations of the hadronic matrix elements of CEDMs and Weinberg-type operators are direly needed.
The goal of this article is to determine the hadronic matrix elements of the following effective operators of Weinberg type [36,[52][53][54][55][56]] O k imply that the correlator is evaluated in the presence of a constant external electromagnetic (EM) source and one of the operators introduced in (1.1). The basic idea is to calculate (2.1) using two different approaches and to match the results to obtain an analytic expression for the nEDM in terms of hadronic quantities. In the first approach, one defines a phenomenological form Π phen of the correlator, which incorporates the wave function of the neutron, its EDM and other parameters. The second approach relies instead on an OPE of the correlator leading to the object Π OPE that depends on the expectation values of effective operators, such as the three-gluon and four-gluon interactions introduced in (1.1). Matching the expressions for Π phen and Π OPE then yields the contribution of the effective operators of interest to the nEDM. To improve the accuracy of the sum rules, the correlators are, however, not matched themselves but their Borel transforms are considered. As we will explain in Section 3, such a procedure removes higher-order polynomial terms and suppresses excited states.
3 Phenomenological side of the sum rules 3

.1 Hadronic representation
In this section, we derive the phenomenological form Π phen of the correlator (2.1) following the argument presented in [33,35,58,59]. An often considered approach for the phenomenological side of two-point correlators is the use of dispersion relations [60][61][62][63][64][65][66]. Since we are interested in the correlator of two nucleon currents η in an external EM field, we are, however, effectively dealing with a three-point correlation function. Dispersion relations for three-point correlators are less constraining than those of two-point correlators due to the lack of positivity constraints [59]. Therefore, we relate the correlator (2.1) to a perturbative expansion of the nucleon propagator in a non-zero EM background. We write Π N (q 2 ) = Π (0) N (q 2 ) + eΠ (1) N (q 2 ) + . . . , (3.1) where e is the electron charge magnitude that serves as the expansion parameter. The first nontrivial term in (3.1) describes the response of the nucleon states to the weak external perturbation and arises from a single insertion of the EM interactions L EM (x) = J µ (x) A µ (x) , J µ (x) = e q=d,u Q qq (x)γ µ q(x) . Here J µ denotes the EM current, A µ is the photon field and Q q is the fractional electric charge of the relevant quark. Note that the EM field is a non-dynamical, classical field in this approach. In order to evaluate the first-order contribution to (3.1), we insert a complete set of hadronic states N and N with the quantum numbers of the neutron into (3.3), i.e. we make use of the identity 1 = N |N N | twice. Working in the so-called fixed-point gauge (see Appendix A), which allows one to express the photon field through the QED field strength tensor employing A µ (y) = −1/2 y ν F µν (0), one obtains eΠ (1) N (q 2 ) = N ,N d 4 x d 4 y e iqx θ(x 0 − y 0 ) θ(y 0 ) 1 2 y ν F µν (0) × Ω| η(x) |N N | J µ (y) |N N |η(0) |Ω + . . . , (3.4) where θ(z) denotes the Heavyside step function, the subscript 0 indicates the time component and the ellipses represent the different combinations due to time ordering. The double sum in (3.4) involves three types of matrix elements of the EM current. These correspond to nucleon transitions of (i) ground state to ground state, (ii) ground state to excited states and vice versa, and (iii) excited states to excited states. Let us first focus on the ground-state contributions, i.e. the terms of the hadronic sums that involve only neutron states |n . Up to an arbitrary chiral phase χ the matrix elements involving |n can be parametrised by the coupling λ between the physical neutron and the interpolating current η as follows Here u is the neutron spinor which satisfies with / p = p µ γ µ , m n denoting the neutron mass andū(p, s) = u † (p, s)γ 0 . Notice that for our correlator (3.4) a spin summation is implicit in the sum over all hadronic states. The product of the matrix element involving the EM current and the photon field can be reduced to a set of four neutron form factors (see for instance [67]) Here q = p 2 − p 1 is the outgoing momentum carried by the photon and σ µν = i/2 [γ µ , γ ν ] where the brackets represent the usual commutator. At q 2 = 0, the form factors in (3.8) can be identified with the fractional electric charge Q n , the magnetic moment µ n , the EDM d n and the anapole moment a n of the neutron. Since the electric charge of the neutron is zero and its anapole moment, as a result of the constant EM background, vanishes as well, one has explicitly It then follows that the tensor structures in (3.8) associated with µ n and d n only differ by a factor iγ 5 , meaning that at zero-momentum transfer one can write with σ · F = σ µν F µν etc. and we have used that γ 5 σ · F = iσ ·F. Inserting (3.5), (3.7) and (3.10) into (3.4) and using (3.6), one obtains for the |n contributions to the first-order correction (3.4) of the nucleon correlation function the following expression with Here the ellipses denote contributions due to exited states and other operators that turn out to be suppressed in the course of our analysis. Up to O(χ) the Lorentz structure in (3.11) behaves under chiral transformations as This result implies that the anti-commutators σ · F, / q and σ ·F, / q are the only structures that are invariant under chiral rotations.

Phenomenological parametrisation
In calculating d n it should then be clear from the above discussion that one should study the operator σ ·F, / q as this structure is the unique choice with an unambiguous coefficient for what concerns the EDM. We thus make the following ansatz [33,35,58,59] , (3.15) for the first-order contribution to the nucleon propagator (3.1). The first term in (3.15) corresponds to the ground-state contribution. It matches the result that we have already derived in (3.11). The second and third term describe transitions of the ground state to excited states and vice versa and transitions of excited states to excited states, respectively. The corresponding form factors are called f N and f N N . They do not have definite signs due to the lack of positivity constraints of the considered correlator [59]. Applying the Borel transformation defined in Appendix B, one finds that the numerically leading contributions of the Borel transforms of the three terms in (3.15) are given by , (3.17) where we have used the values m n ' 0.94 GeV and m N 0 ' 1.44 GeV [68] for the mass of the neutron and its lightest excitation to obtain the quoted numerical prefactor. Under the assumption that | f N 0 | ' | f N 0 N 00 | ' 2 m n d n and setting M = 2⇤ QCD ' 0.6 GeV with ⇤ QCD the QCD scale, the mixed ground-state and exited-state (excited-states only) contributions therefore naively amount to relative corrections of the order of 30% (5%). In the following, we only include the ground-state contribution to (3.14) in our sum-rule calculation, and estimate the uncertainties that are associated to this simplification by a variation of the Borel mass M (cf. Section 6). The appropriate form of the phenomenological side of our sum rule can be established by realising that the contributions to the nEDM induced by CP-violating Weinberg-type operators (1.1) have a simple pictorial interpretation [37]. As illustrated in Figure 1, there are two types of graphs that one needs to consider in general. The first type of diagrams (left and middle) factorises into a propagator with a CP-violating mass insertion proportional to i 5 and into a part that couples to the external photon field. The e↵ect of Weinberg-type operators in this context is to rotate the nucleon wave function by an amount proportional to d n /µ n as in (3.10). The second type of diagrams (right) only exists if either an insertion of an operator is considered that couples several gluons to a single photon or if at least one of the external legs corresponds to an excited state [37]. The former possibility is not viable at the dimension-six level, because there is no gauge-invariant operator Applying the Borel transformation defined in Appendix B, one finds that the numerically leading contributions of the Borel transforms of the three terms in (3.15) are given by where we have used the values m n 0.94 GeV and m N 1.44 GeV [68] for the mass of the neutron and its lightest excitation to obtain the quoted numerical prefactor. Under the assumption that | f N | | f N N | λ 2 m n d n and setting M = 2Λ QCD 0.6 GeV with Λ QCD the QCD scale, the mixed ground-state and exited-state (excited-states only) contributions therefore naively amount to relative corrections of the order of 30% (5%). In the following, we only include the ground-state contribution to (3.14) in our sum-rule calculation, and estimate the uncertainties that are associated to this simplification by a variation of the Borel mass M (cf. Section 6).
The appropriate form of the phenomenological side of our sum rule can be established by realising that the contributions to the nEDM induced by CP-violating Weinberg-type operators (1.1) have a simple pictorial interpretation [37]. As illustrated in Figure 1, there are two types of graphs that one needs to consider in general. The first type of diagrams (left and middle) factorises into a propagator with a CP-violating mass insertion proportional to iγ 5 and into a part that couples to the external photon field. The effect of Weinberg-type operators in this context is to rotate the nucleon wave function by an amount proportional to d n /µ n as in (3.10). The second type of diagrams (right) only exists if either an insertion of an operator is considered that couples several gluons to a single photon or if at least one of the external legs corresponds to an excited state [37]. The former possibility is not viable at the dimension-six level, because there is no gauge-invariant operator that couples two gluons to a single photon. In the approximation that neglects the contributions of vertex diagrams and excitations, one can therefore use the following parameterisation Π phen (q 2 ) = − λ 2 m 2 n µ n 2 q 2 − m 2 n 1 + r(q 2 ) iγ 5 , (3.18) where the coefficient function r(q 2 ) has to be determined by matching the phenomenological side of the sum rule to the corresponding OPE calculation. Since from (3.10) we know that the EDM and the magnetic moment of the neutron are simply related by a chiral rotation with iγ 5 though, the following relation holds d n = µ n r(q 2 ) . In physical terms this result means that the Weinberg-type contributions to d n can be approximated by calculating the iγ 5 rotation of the nucleon wave function and relating it to the corresponding chiral rotation of µ n [37,38]. In Section 4 and Section 5 we will use (3.18) and (3.19) to extract the hadronic matrix elements of O 6 and O 8 , respectively.

Interpolating current
We parameterise the interpolating current introduced in (2.1) as follows where the real parameter β is kept arbitrary throughout our calculations. The two currents form a basis for projection onto the neutron state in the case of a CP-conserving background. The current j 1 (x) is often used in LQCD simulations to describe the neutron wave function (see for instance [69][70][71]). While j 2 (x) vanishes in the non-relativistic limit, it should be included in the interpolating field since we are dealing with light quarks. In (4.2) the symbols a, b and c are colour indices, d and u denote a down-quark and up-quark field, respectively, and C is the charge conjugation matrix, which satisfies C = C * = −C † = −C T = −C −1 , C γ T 5 C = −γ 5 and C † γ 0 = γ 0 C. Notice that in contrast to the publications [33,35,58,59], we do not need to consider the two additional currents i 1 (x) = γ 5 j 2 (x) and i 2 (x) = γ 5 j 1 (x), because in our case the only source of CP-violation is provided by the Weinberg-type operators (1.1). The vacuum |Ω appearing in correlators such as (2.1) is instead taken to be CP-conserving, which in particular means that we assume that the QCD theta term θG A µν G A µν vanishes either accidentally or dynamically due to a Peccei-Quinn mechanism [72].

Weinberg contribution to the quark propagator
In the presence of a non-trivial EM background and the dimension-six operator O 6 , the OPE of the correlator (2.1) can be formally written as where C k are so-called Wilson coe cients and hQ k i = h⌦|Q k |⌦i are vacuum matrix elements or condensates of the operator Q k . One important ingredient to evaluate (4.3) is the quark propagator on the CP-conserving background including insertions of the EM field and O 6 . In position space and suppressing colour and spinor indices the sought propagator reads where the first term is the free propagator for a massless quark and the second term describes non-perturbative interactions with background quark fields. As shown in Appendix C at leading order (LO) in the OPE these two quantities take the following form with hqqi ' (0.25 GeV) 3 [65,66,73] the quark condensate. The e↵ective operator O 6 can be perturbatively inserted into the quark propagator [38]. The corresponding Feynman diagram is shown in Figure 2. It follows that the Weinberg-induced contribution to (4.4) can be written as where the amputated two-point function S O 6 amp (z) is given by  where C k are so-called Wilson coefficients and Q k = Ω|Q k |Ω are vacuum matrix elements or condensates of the operator Q k . One important ingredient to evaluate (4.3) is the quark propagator on the CP-conserving background including insertions of the EM field and O 6 . In position space and suppressing colour and spinor indices the sought propagator reads where the first term is the free propagator for a massless quark and the second term describes non-perturbative interactions with background quark fields. As shown in Appendix C at leading order (LO) in the OPE these two quantities take the following form with qq −(0.25 GeV) 3 [65,66,73] the quark condensate. The effective operator O 6 can be perturbatively inserted into the quark propagator [38]. The corresponding Feynman diagram is shown in Figure 2. It follows that the Weinberg-induced contribution to (4.4) can be written as where the amputated two-point function S O 6 amp (z) is given by where C k are so-called Wilson coe cients and hQ k i = h⌦|Q k |⌦i are vacuum matrix elements or condensates of the operator Q k .
One important ingredient to evaluate (4.3) is the quark propagator on the CP-conserving background including insertions of the EM field and O 6 . In position space and suppressing colour and spinor indices the sought propagator reads where the first term is the free propagator for a massless quark and the second term describes non-perturbative interactions with background quark fields. As shown in Appendix C at leading order (LO) in the OPE these two quantities take the following form with hqqi ' (0.25 GeV) 3 [65,66,73] the quark condensate. The e↵ective operator O 6 can be perturbatively inserted into the quark propagator [38]. The corresponding Feynman diagram is shown in Figure 2. It follows that the Weinberg-induced contribution to (4.4) can be written as where the amputated two-point function S O 6 amp (z) is given by Here g s denotes the QCD coupling constant, t A are the colour generators of SU(3) and we have expanded the quark wave function to zeroth order using (A.7) to obtain the final result. The object D O 6 AB µ⌫ (z) entering (4.7) represents the Weinberg-induced correction of the gluon propagator. Pictorially, one has where the dotted vertex represents the insertion of the operator O 6 and the cross indicates interactions with the classic background. In order to determine the form of (4.8) we rely on standard background field techniques. We start by writing the dimension-six Weinberg operator of (1.1) in a more convenient form, namely as (see for instance [53,54,56]) where T µνρλστ denotes the following trace Tr σ µν σ ρλ σ στ γ 5 . (4.10) Notice that this tensor is anti-symmetric under µ ↔ ν, etc. as well as µν ↔ ρλ etc. By splitting the original gluon field G A µ =Ḡ A µ +Ĝ A µ into a classical fieldḠ A µ and quantum fieldĜ A µ , one can then expand the QCD field strength tensor around its classical configuration to obtain When one now expands (4.9) using (4.11), one is only interested in terms that are linear inḠ A µ and bilinear inĜ A µ . Using the anti-symmetric properties of f ABC and that of (4.10), we find that the relevant terms are Employing the result (4.12) one can now calculate the Weinberg-induced corrections (4.8) to the gluon propagator. By performing all possible contractions of the time-ordered product, we obtain the expression where D (0)AB µν (x) denotes the free gluon propagator in position space. In Feynman gauge it takes the following simple form a result that can be gleaned by inspection of (D.3). Here g µν = diag (1, −1, −1, −1) denotes the Minkowski metric and we have introduced the abbreviation d 4 p = d p 4 /(2π) 4 . In order to simplify (4.13) we use the following two relations which follow from the expansion (A.6) and the explicit form (4.14) of the free gluon propagator in momentum space, respectively. Using (4.15) yields where in the last step we have employed the Fourier integral given in (D.6). Plugging this into (4.7) we find (4.17) Here we have used that the non-perturbative quark-gluon condensate appearing in the first line simplifies as follows (see for instance [71,74]) Furthermore, the colour factor in the second line evaluates to with the Casimir operators given by C A = 3 and C F = 4/3 for SU (3). From (D.9) one sees that the Fourier transform of (4.17) reads where we have dropped colour and spinor indices. Inserting this result into (4.6) leads to This result agrees up to a factor of i with the corresponding expression reported in [38] after taking into account that the definition of O 6 used in this paper differs from the one employed in (1.1) by a factor of 1/3.

OPE including the Weinberg operator
In terms of (4.4) and the general form of the correlator (4.3) can be written as where we have performed all possible Wick contractions. Here β is the real parameter that appears in our definiton (4.1) of the interpolating current for the neutron.
We are only interested in the LO result of the OPE, in other words in terms that are linear either in the quark condensate (4.5) or the Weinberg-induced correction (4.21) to the quark propagator. The relevant contributions are given by the following expression This result can be interpreted in terms of the two Feynman diagrams depicted in Figure 3. The left graph shows the S q (x) part of (4.25) which corresponds to a one-loop diagram because the background quark fields are non-dynamical. The right diagram has instead a dynamically and perturbative gluon that closes the second loop. In the case of the S O 6 (x) correction one therefore has to deal with a two-loop contribution.

Matching and discussion
In order to derive the sum rules for the O 6 contribution to the nEDM, we match the phenomenological and the OPE correlators, i.e. we set (3.18) and (4.28) equal and determine the coe cient r(q 2 ) that appears in front of the i 5 term in ⇧ phen (q 2 ). After Borel transformation and identifying the IR cut-o↵ with the QCD scale, we obtain With this result at hand, we now discuss the appropriate choice for the mixing parameter introduced in (4.1) for our sum rule. There are two commonly used ways for fixing this parameter: (i) at a value where the leading terms of the OPE are stationary under variations of or (ii) at a point that maintains a balance between OPE convergence and contributions of excited states. Both method are not applicable in our case, because (i) the result (4.29) does not possess an extremum and (ii) contributions of excited states have been ignored in our sum rule (cf. Section 3.2). The procedure advocated in [33,35,58,59] where is chosen such that subleading IR logarithms are cancelled in the QCD ✓-term and CEDM contributions to the nEDM is also not useful, since in the OPE side (4.28) of our sum rule IR logarithms appear already at LO.
Our choice of is instead based on the observation that the function f q ( ) introduced in (4.26) appears in the numerator of Io↵e's formula [76], which for arbitrary takes the following form [71] m n = 7 2

Matching and discussion
In order to derive the sum rules for the O 6 contribution to the nEDM, we match the phenomenological and the OPE correlators, i.e. we set (3.18) and (4.28) equal and determine the coefficient r(q 2 ) that appears in front of the iγ 5 term in Π phen (q 2 ). After Borel transformation and identifying the IR cut-off with the QCD scale, we obtain With this result at hand, we now discuss the appropriate choice for the mixing parameter β introduced in (4.1) for our sum rule. There are two commonly used ways for fixing this parameter: (i) at a value where the leading terms of the OPE are stationary under variations of β or (ii) at a point that maintains a balance between OPE convergence and contributions of excited states. Both methods are not applicable in our case, because (i) the result (4.29) does not possess an extremum and (ii) contributions of excited states have been ignored in our sum rule (cf. Section 3.2). The procedure advocated in [33,35,58,59] where β is chosen such that subleading IR logarithms are cancelled in the QCD θ-term and CEDM contributions to the nEDM is also not useful, since in the OPE side (4.28) of our sum rule IR logarithms appear already at LO.
Our choice of β is instead based on the observation that the function f q (β) introduced in (4.26) appears in the numerator of Ioffe's formula [76], which for arbitrary β takes the following form [71]  the formula (4.30) thus predicts a neutron mass that is in the ballpark of the experimental measured value m n 0.94 GeV. We conclude from this that the appropriate choice for the mixing parameter in the case of our sum rule (4.28) is β = −1. In fact, this choice is the one that has been employed in essentially all CP-even sum rules [75][76][77], including the sum-rule calculations of the anomalous magnetic moment µ n of the neutron [61,62,78,79]. We believe that the Ioffe interpolating current has also been used in [38].

Weinberg contribution to the quark propagator
In this section, we derive the contribution to the nEDM from the dimension-eight Weinberg operator O 8 introduced in (1.1). Like in Section 4.2 we will treat the CP-violating four-gluon operator as a perturbative insertion into the quark propagator. The corresponding Feynman graph is shown in Figure 4. In analogy to (4.4), (4.6) and (4.7) we write and The explicit LO expressions for S (0) (x) and S q (x) can be found in (4.5). The function D O 8 AB µν (z) in (5.3) corresponds to the correction of the gluon propagator due to an insertion of O 8 , namely value m n ' 0.94 GeV. We conclude from this that the appropriate choice for the mixing parameter in the case of our sum rule (4.28) is = 1. In fact, this choice is the one that has been employed in essentially all CP-even sum rules [75][76][77], including the sum-rule calculations of the anomalous magnetic moment µ n of the neutron [61,62,78,79]. We believe that the Io↵e interpolating current has also been used in [38].

Weinberg contribution to the quark propagator
In this section, we derive the contribution to the nEDM from the dimension-eight Weinberg operator O 8 introduced in (1.1). Like in Section 4.2 we will treat the CP-violating four-gluon operator as a perturbative insertion into the quark propagator. The corresponding Feynman graph is shown in Figure 4. In analogy to (4.4), (4.6) and (4.7) we write and To determine the analytic expression corresponding to (5.4), we proceed as in Section 4.2 and write the operator O 8 in the more convenient form These traces are anti-symmetric under the exchanges µ $ ⌫ etc. but symmetric under the exchanges µ⌫ $ ⇢ etc. The e↵ective operator O 8 is then expanded in terms of partial derivatives and gluon fields using (4.11). Picking out the terms that are bilinear in both the background fieldḠ A µ and the quantum fieldĜ A µ , we obtain To determine the analytic expression corresponding to (5.4), we proceed as in Section 4.2 and write the operator O 8 in the more convenient form These traces are anti-symmetric under the exchanges µ ↔ ν etc. but symmetric under the exchanges µν ↔ λρ etc. The effective operator O 8 is then expanded in terms of partial derivatives and gluon fields using (4.11). Picking out the terms that are bilinear in both the background fieldḠ A µ and the quantum fieldĜ A µ , we obtain  where the ellipses in the bracket represent the other three terms that are quadratic inḠ A µ andĜ A µ . With the expression (5.7) at hand it is a matter of simple algebra to calculate (5.4). Using the relations in (4.15) we find Here we have used the shorthand notationG A · G B =G A µ⌫ G B µ⌫ and all field strength tensors and duals are evaluated at 0. The ellipses represent three additional terms that have a structure that is similar to that of the contribution proportional to p p ⇢ .
Inserting (5.8) into (5.3) it then turns out that to obtain the amputated two-point function S O 8 amp (z), one has to calculate objects of the form where X CD is a symmetric tensor in colour space. To achieve this we expand the quark current in terms of the set of basis matrices using the Fierz identity (see for instance [80,81]) In the case of the structure M AB ⇢ , we obtain where a sum over the five di↵erent Lorentz structures in (5.10) is implicit.
-14 - where the ellipses in the bracket represent the other three terms that are quadratic inḠ A µ andĜ A µ . With the expression (5.7) at hand it is a matter of simple algebra to calculate (5.4). Using the relations in (4.15) we find Here we have used the shorthand notationG A · G B =G A µν G B µν and all field strength tensors and duals are evaluated at 0. The ellipses represent three additional terms that have a structure that is similar to that of the contribution proportional to p λ p ρ .
Inserting (5.8) into (5.3) it then turns out that to obtain the amputated two-point function S O 8 amp (z), one has to calculate objects of the form where X CD is a symmetric tensor in colour space. To achieve this we expand the quark current in terms of the set of basis matrices Γ n = 1, γ 5 , γ α , iγ 5 γ α , σ αβ , (5.10) using the Fierz identity (see for instance [80,81]) In the case of the structure M AB λρ , we obtain where a sum over the five different Lorentz structures in (5.10) is implicit.
Recalling from (4.17) that we are only interested in the pieces of S O 8 amp (z) that are proportional to iγ 5 it is readily seen that these contributions all arise from the term γ µ γ 5 γ ν = −g µν γ 5 + 1 2 µνλρ σ λρ , (5.13) in the above expression for M AB λρ . The colour factors appearing in (5.12) can be decomposed in the following way with

(5.16)
Here we have used the identity that holds for any symmetric X AB and introduced the following shorthand notation for the two types of condensates appearing in (5.16). A calculation similar to the one detailed above leads to To determine the expression (5.3), we finally specify the possible colour structures c ABCD that appear in the definition (1.1) of O 8 . There are in fact three independent CP-violating four-gluon operators [52,55,56] and the corresponding colour structures can be chosen to be (5.20) Notice that these coefficients are symmetric under the simultaneous exchange of A ↔ B and C ↔ D and the pairwise exchange AB ↔ CD. Using (5.16), (5.19) and the properties of c ABCD , we find that (5.3) can be written as follows with the new colour structures In fact, using the identities (see for instance [82]) A comparison of (5.21) with (4.17) and (4.21) then implies that the Weinberg-induced correction (5.2) to the quark propagator in position space can be written as with the coefficients c 1 and c 2 given in (5.24). Notice that the dimension-seven condensates appearing in (5.25) are the only non-zero matrix elements that can be constructed out of two quark fields, a QCD field strength tensor and its dual [83,84]. This finding provides a sanity check of the calculations leading to S O 8 (x).

OPE correlator, matching and discussion
To determine the OPE correlator (4.3) corresponding to the dimension-eight Weinberg-type operator O 8 , we also need values for the two condensates in (5.25). The only estimates that exist at present are based on the instanton liquid model [66,[85][86][87][88]. One obtains [84] qiγ 5 qg 2 sG · G = 64 5ρ 4 qq , qt A iγ 5 qg 2 s d ABCG B · G C = 32 15ρ 4 qq . In the diluted instanton gas model, the quark condensate is given by where m q Λ QCD 0.3 GeV denotes the constituent quark mass andρ 1/(0.6 GeV) is the average instanton size. Notice that for the quoted values of m q andρ one finds qq −(0.25 GeV) 3 in agreement with the standard value for the quark condensate [65,66,73].
Noticing that after employing the relations (5.26) the structure of (5.25) and (4.21) are precisely the same, the derivation of Π OPE (q 2 ) and the matching of the phenomenological and the OPE correlators for O 8 proceeds as in Sections 4.3 and 4.4, respectively. In particular, for the coefficient r(q 2 ) that multiplies the iγ 5 term in (3.18), we obtain with the functions f q (β) and f O (β) defined in (4.26). Like in the case of O 6 , we will employ β = −1 in our numerical analysis of the O 8 matrix elements, since this is the appropriate choice for our sum-rule calculations (cf. the discussion at the end of Section 4.4).

.1 Dimension-six contribution
Using β = −1 and inserting (4.29) into (3.19), we obtain the following expression for the contribution of the dimension-six Weinberg operator to the nEDM which differs by the analytic result given in [38] by a sign.
In our numerical analysis we use , g s = 4πα s = 2.13 ± 0.03 , where the input values and errors of µ n , m p , g s and m 2 0 are taken from [65,66,68,73] and the strong coupling constant corresponds to a LO α s evaluated at a renormalisation scale of 1 GeV. We note that our choice M ∈ √ 2 [1, 2] Λ QCD [0.42, 0.85] GeV covers the full range of Borel masses that has been considered in the related sum-rule calculations of the QCD θ-term and CEDM contributions to the nEDM [33,35,58,59].
With the input given in (6.2) we find from (6.1) the following numerical result where the individual uncertainties have been added in quadrature to obtain the final relative error of 50%. The dominant source of uncertainty in our prediction for (d n /e) O 6 arises from the variation of the scale ratio M/Λ QCD and amounts to almost 90% of the total error given above. We add that the quoted total uncertainty in (6.1) is larger than the naive expectation of the size of the sum-rule contributions due to excited states (cf. the discussion at the end of Section 3.2) and that the sum-rule predictions for the down-quark and up-quark CEDMs [33,35] are also accurate to about 50%. Notice that the central value of our prediction (6.1) differs by a factor of roughly 1/3 from the numerical result presented in [38]. Here a factor of 1/3 is accounted for by the different normalisation of the effective operator O 6 , while the flipped overall sign in (6.1) is compensated by the fact that in the latter article the incorrect relation µ n = 1.91 e/(2m p ) has been used to obtain a numerical result.

Dimension-eight contributions
Inserting ( where the numbers in the curly bracket correspond to the three different colour structures in (5.20).
In order to obtain a numerical result for the O 8 contribution to the nEDM, we use the input given in (6.2) and (6.5). Adding individual uncertainties in quadrature we find Here the dominant source of uncertainty stems again from the variation of M/Λ QCD and amounts close to 60% of the quoted total error of 80%.

Conclusions
In this work, we have calculated the hadronic matrix elements of dimension-six and dimensioneight operators of Weinberg type (1.1) using QCD sum-rule techniques. Calculations along the same line of the dimension-four and dimension-five contributions to the nEDM, i.e. the QCD θterm and CEDMs, have been performed in [33-35, 58, 59]. A sum-rule estimate of the dimensionsix Weinberg operator O 6 also exists [38], but this article does not provide details on the actual computation, which motivated us to carry out an independent evaluation. Our determination of the hadronic matrix elements of the dimension-eight term O 8 is instead new, and provides the first systematic study of contributions to d n due to CP-violating four-gluon operators. The main results of our article, i.e. the numerical expressions (6.3) and (6.6), will be used in a companion paper [57] to derive model-independent bounds on CP-violating Higgs-gluon interactions in BSM scenarios with vanishing light-quark Yukawa couplings. Our sum-rule estimates are based on the observation [37,38] that the Weinberg-type contributions to the nEDM can be obtained by calculating the iγ 5 rotation of the nucleon wave function induced by (1.1) and relating it to the corresponding rotation of the neutron anomalous magnetic moment µ n . In this approximation only diagrams are included that factorise into a propagator with a CP-violating mass insertion and into a part that couples to the external photon field, while nonfactorisable vertex corrections are neglected (see Figure 1). In addition, we neglect contributions from excited neutron-like states N in our estimates. These simplifications lead to uncertainties in our predictions that we estimate to be of the order of 35% using Borel techniques. The OPE computation of the dimension-six and dimension-eight sum rules is described in detail, and includes a discussion of the matching and the appropriate choice of the neutron-interpolating current. The final analytic expressions for the O 6 and O 8 contributions to d n are reported in (6.1) and (6.4), respectively, and our result for (d n ) O 6 is found to agree with that given in [38]. The hadronic matrix elements of the Weinberg-type operators turn out to be logarithmically sensitive to the ratio M/Λ QCD of the Borel mass and the QCD scale. This IR sensitivity provides the dominant theoretical uncertainty of our predictions. By varying M/Λ QCD in the range √ 2 [1,2] we find uncertainties close to 45%, which exceeds the size of the expected effects from vertex diagrams and excited states. Adding individual errors in quadrature the total uncertainties of our numerical predictions for (d n /e) O 6 and (d n /e) O 8 are 50% and 80%. Sum-rule studies of the θ-term and CEDM contributions to the nEDM [33-35, 58, 59] have also found uncertainties of a similar magnitude.
While our sum-rule estimates of the d n contributions due to the dimension-six and dimensioneight operators of Weinberg type (1.1) have sizeable uncertainties, we believe that they are more robust than other existing determinations that are based on NDA [36] or the VIA [37]. In particular, in the sum-rule approach there is no sign ambiguity between the prediction for d n and the hadronic matrix elements of O 6 and O 8 -see for instance [6,13,15,16,19] for EDM studies that allow for both signs of the O 6 contribution. To find out whether our sum-rule estimates are reliable would require first-principle calculations of the nEDM, which are in principle possible using existing LQCD methodology. While such calculations have gained significant momentum [39][40][41][42][43][44][45][46], LQCD simulations involving Weinberg-type operators are challenging [47,48], and it remains to be seen which accuracy such computations can achieve in the near future. Till then phenomenological studies of hadronic EDMs have to rely on the predictions (6.3) and (6.6) despite all their limitations.

Acknowledgments
The Feynman diagrams shown in this article have been drawn with the help of the L A T E X package feynMF [95]. Parts of the Dirac and colour algebra calculations have been cross-checked against FeynCalc [96] and TRACER [97].

A Fixed-point gauge
In the case of QCD the fixed-point or Fock-Schwinger gauge [98,99] can be defined without loss of generality for gauge-invariant quantities by where it is sufficient to restrict this choice of gauge to classical gluon fields G A µ . For the gauge choice (A.1) it is easy to show that one can express the gluon field through the QCD field strength tensor Here g s denotes the strong coupling constant. To derive the sought relation, one notices first that where we have employed the gauge condition (A.1) twice and used the anti-symmetry of G A µν to obtain the final result. Setting x ν = αy ν with an arbitrary parameter α, one can then write Now if one integrates both sides of the above relation over α ∈ [0, 1] and assumes that G A µ (x) is non-singular at x = 0, one finds Using a similar assumption for the QCD field strength, one can Taylor expand G A νµ (αy) around αy = 0 and perform the integration on the left-hand side of (A.4). It follows that where we have switched back from the variable y to the variable x. The latter expression can be further simplified by noting that as a result of (A.1), the partial derivatives in (A.5) can be promoted to covariant derivatives D µ = ∂ µ − ig s G A µ t A . In consequence, one has in the fixed-point gauge of QCD. The same result also holds in the case of EM with G A µ replaced by A µ and G A µν replaced by F µν . Finally, notice that due to (A.1) the Taylor expansion of the classical quark field q can also be formulated in terms of covariant rather than partial derivatives. One has

B Borel transforms
We define the Borel transformation of a function F(Q 2 ) with Q 2 = −q 2 in the following way The auxiliary parameter M is called the Borel mass and occurs in the final result of most sum-rule calculations. A collection of useful Borel transforms can be found for instance in [63,64,74] and reads Here k ∈ N + and Γ(z) denotes the Euler gamma function. From (B.2) it is clear that all polynomial contributions to sum-rule correlation functions vanish after Borel transformation.

C OPE for the quark propagator
The free position-space propagator S (0) (x) of a massless quark is easily derived by taking the Fourier transform of the well-known momentum-space representation where colour indices are implicit and we have applied a Wick rotation to Euclidean space using x 2 → −x 2 E and p 2 → −p 2 E . The four-dimensional integration measure can be written as d 4 p E = p 3 dp dΩ 4 , where the magnitude of the four-dimensional Euclidean momentum vector has been denoted byp = |p E |. The differential solid angle is given by dΩ 4 = 2π 0 dφ π 0 dθ 1 dθ 2 sin θ 1 sin 2 θ 2 . (C. 2) It follows that dΩ 4 = dΩ 3 sin 2 θ 2 and (C.1) hence takes the form

(C.3)
Herex = |x E |, J 1 (z) denotes the Bessel function of first kind with index 1 and in the final step we have rotated back from Euclidean to Minkowski space noting that / x E → / x.
In order to determine the non-perturbative contributions S q (x) to the quark propagator, one needs to evaluate the correlator Ω|T q i a (x)q where the fields of the condensate qq are evaluated at x = 0. Notice that the minus sign in the final result comes from the exchange of the fermion fields and the numerical prefactor can be determined by contracting the expression in the middle and on the right with δ ab δ i j . Ignoring colour and spinor indices, the expansion (C.4) thus leads to the expression for S q (x) as given in (4.5).

D Fourier transforms
We define the Fourier transform of a function F(p) by where we have performed the integration in d = 4 + 2 IR space-time dimensions with IR > 0 to regulate the IR divergences that appear in some of the Fourier integrals that we have encountered in Sections 4.3 and 5.1. The symbol γ E denotes the Euler's constant and µ IR is a mass scale needed to restore the correct dimensionality of (D.1).
In the case that F(p) is polynomial in 1/p 2 , a simple calculation along the lines of the computation performed in Section C leads to F        1 where k ∈ N + . The Fourier transforms of type (D.2) relevant for our article are Tensor Fourier integrals with a polynomial denominator of the form (p 2 ) k can be obtained from (D.2) by taking derivatives This procedure allows one to derive for example If F(x) is polynomial in 1/x 2 , it is straightforward to evaluate (D.7). For k ∈ N + we obtain F        1 The Fourier integrals of the form (D.8) that occur in our article are We also encounter in our work Fourier transforms that are both IR and UV divergent. They are of the type F 1/ (x 2 ) l F 1/ (p 2 ) k with k, l ∈ N + . Using the result given in (D.2) and (D.8) these double integrals are readily computed. We find F        1 (D.11) The only Fourier integral of the form (D.11) that is necessary to compute the two-loop contributions to the OPE correlation functions considered in this paper is where the ellipses represent terms that vanish after Borel transformation (cf. Appendix B), meaning that these contributions do not enter the analytic expressions (6.1) and (6.4).