Improved indirect limits on charm and bottom quark EDMs

We derive indirect limits on the charm and bottom quark electric dipole moments (EDMs) from paramagnetic AMO and neutron EDM experiments. The charm and bottom quark EDMs generate $CP$-odd photon-gluon operators and light quark EDMs at the $c$- and $b$-quark mass thresholds. These $CP$-odd operators induce the $CP$-odd semi-leptonic operator $C_S$ and the neutron EDM below the QCD scale that are probed by the paramagnetic and neutron EDM experiments, respectively. The bound from $C_S$ is $\vert d_c \vert<1.3\times 10^{-20}\,e\,\mathrm{cm}$ for the charm quark and $\vert d_b \vert<7.6\times 10^{-19}\,e\,\mathrm{cm}$ for the bottom quark, with its uncertainty estimated as 10%. The neutron EDM provides a stronger bound, $\vert d_c \vert<6\times 10^{-22}\,e\,\mathrm{cm}$ and $\vert d_b \vert<2\times 10^{-20}\,e\,\mathrm{cm}$, though with a larger hadronic uncertainty.


Introduction
Searches for electric dipole moments (EDMs) play a crucial role in probing physics beyond the Standard Model (SM). Given the smallness of the SM contribution from the Cabbibo-Kobayashi-Maskawa phase, EDMs of elementary particles are ideal observables to look for new physics that violates CP symmetry (see [1] for a review).
EDMs of elementary particles may be classified into three categories: lepton EDMs, light quark EDMs and heavy quark EDMs, where "light" and "heavy" are in comparison with the QCD scale. Among the lepton EDMs, the electron EDM is tested by paramagnetic EDM experiments such as the ACME experiments [2], while the muon and tau EDMs can be probed indirectly through paramagnetic and diamagnetic EDM experiments [3][4][5]. Motivated by the muon g − 2 anomaly [6,7], the muon EDM has been and will be measured by the storage ring experiments as well [8][9][10][11][12]. Hadronic EDMs can be induced by the EDMs of light quarks, as well as by other operators such as chromo-EDMs (CEDMs) and three-gluon CP -odd operators. The light quark EDMs, in particular the up and down quark EDMs, contribute to the neutron EDM and are constrained by the neutron EDM experiments [13,14]. Discussion of the strange quark contribution was typically focused on the color EDM contributions [15,16], while no reliable estimate of the neutron EDM induced by the strange quark exists at this point. Among the heavy quark EDMs, the top quark EDM is arguably most deeply investigated. In particular, due to its large Yukawa coupling, the top quark EDM generates light fermion EDMs through diagrams with an intermediate Higgs, which allows one to indirectly constrain the top quark EDM [17,18]. The light quark CEDMs directly contribute to the neutron EDM, while the heavy quark CEDMs generate, e.g., the three-gluon Weinberg operator at the threshold, which induce the neutron EDM after the confinement [19][20][21][22][23][24].
Direct measurements of the heavy quark EDMs are difficult due to their short lifetimes. The current strongest direct bound is from e + e − → qq at LEP and is only of the order of 10 −17 e cm [25].
Recently, an LHC based experiment is proposed that aims at directly measuring charm baryon dipole moments [26][27][28][29][30], potentially improving the direct bounds on the heavy quark EDMs. Given this situation, in this paper, we address indirect limits on the heavy quark EDMs. The heavy quark EDMs induce several CP -odd operators below the heavy quark mass scale. These operators generate observable signals in neutron/atomic/molecule EDM experiments, which allows us to impose indirect limits on the heavy quark EDMs. Thus, our goal of this paper is to understand the current status of indirect limits on the heavy quark EDMs.
Indirect limits on the charm and bottom quark EDMs were previously considered in [31] and the CEDM case was also analyzed in [32]. There the authors derived limits based on the constraints on the heavy quark CEDMs. Indeed, the charm and bottom quark EDMs well above the heavy quark mass scales, say 1 TeV, generate the CEDM operators through the renormalization group (RG) running at a lower energy scale. These CEDMs in turn source the three-gluon Weinberg operator after integrating out the charm and bottom quark. The Weinberg operator then generates the neutron EDM, and this allows one to derive limits on the charm and bottom quark EDMs. We may phrase it a re-interpretation of bounds on the heavy quark CEDMs at the heavy quark mass threshold, with an assumption that there is no cancellation between the EDM and CEDM contributions (see also discussion at the end of Sec. 3). In contrast, in this paper, we study CP -odd operators generated from the EDM operators at the heavy quark mass threshold. Therefore our consideration directly applies to the heavy quark EDMs at the quark mass threshold and is independent of [31]. Even though the CP -odd operators induced by the heavy quark EDMs are formally higher dimensional, suppressions from the charm and bottom quark masses are not severe as there is only a little hierarchy between the quark masses and the QCD scale.
In order to make our limits robust, we derive indirect limits on charm and bottom quark EDMs based on multiple observables, paramagnetic and neutron EDMs. Indirect limits always have a potential of having a cancellation among different operators. For instance, paramagnetic EDM experiments are sensitive to only a particular linear combination of the electron EDM d e and the CP -odd semi-leptonic operator C S (see Sec. 3). Therefore a new physics contribution to d e and C S can be such that the linear contribution almost vanishes even though each term is sizable. Deriving constraints from two completely different observables, paramagnetic and neutron EDMs, makes the probability of having such a cancellation less likely. We also minimize the QCD uncertainty as much as possible. For this purpose the paramagnetic EDM plays an essential role as it is sensitive to the heavy quark EDM through the semi-leptonic operator C S . The estimation of C S is not polluted by hadronic uncertainties, and its uncertainty is theoretically well under control.
The structure of this paper is as follows. In Sec. 2, we derive CP -odd operators generated by a heavy quark EDM after integrating out the heavy quark. In particular, we see that the heavy quark EDM generates CP -odd photon-gluon operators and light quark EDMs. We restrict ourselves above the QCD scale in this section, and thus the operators are written in terms of light quarks and gluons. In Sec. 3, we calculate the semi-leptonic CP -odd operator C S induced from the CP -odd photon-gluon operator below the QCD scale. This allows us to derive a clean constraint on the heavy quark EDMs from the paramagnetic EDM experiments. Sec. 4 is devoted to an estimation of the neutron EDM induced by the CP -odd photon-gluon operators and the light quark EDMs. There we derive limits on the heavy quark EDMs from the neutron EDM experiments, which are stronger than the paramagnetic EDM experiments but have a larger hadronic uncertainties. The conventions of this paper is summarized in App. A, while some technical details of our computation can be found in App. B. The effective action to linear order in d Q after integrating out the heavy quark. The thick line on the left hand side indicates the full heavy quark propagator with the covariant derivative D µ , and the cross dot indicates the EDM operator insertion. The full propagator is expanded with respect to the field strengths, which results in the CP -odd photon-gluon operators as shown on the right hand side.

CP -odd operators from heavy quark EDMs
In this section, we compute CP -odd operators induced by integrating out charm and bottom quarks with EDMs. The Lagrangian is given by where Q is the heavy quark with its mass m Q and its EDM d Q . The covariant derivative is defined by where G a µ is the gluon field with its field strength G a µν and coupling g s , and A µ is the photon field with its field strength F µν and coupling e. The SU(3) generator is denoted by T a while the U(1) charge by Q Q , with Q c = Q t = +2/3 and Q b = −1/3. See App. A for more details on our conventions. Throughout this paper, "heavy" means that the mass is larger than the QCD scale. Although our discussion applies equally to the top quark, it has another contribution that puts a stronger constraint on the EDM [17,18]. Therefore our main focus is on the charm and bottom quarks, i.e., Q = c, b.
Once we integrate out the heavy quark, d Q generates several CP -odd operators. For our purpose, there are two types of operators that are important: CP -odd photon-gluon operators and light quark EDMs. Below the QCD scale, the former generates the neutron EDM and the semi-leptonic CP -odd interaction C S that induces paramagnetic EDMs, while the latter contributes to the neutron EDM. In this section, we keep ourselves above the QCD scale. The operators are then written in terms of the quarks and gluons. We will deal with nonperturbative effects at the QCD scale in the subsequent sections.

CP -odd photon-gluon operators
The heavy quark EDM induces CP -odd photon-gluon operators through the diagrams in Fig. 1. The effective action after integrating out the heavy quark, to linear order in d Q , is given by where the trace is taken over the spinor, the color and the spacetime. We further expand this expression with respect to the field strength H µν , given by where tr c is the trace only over the color. This contains CP -odd photon-gluon operators as well as a CP -odd pure photon operator. The operator quadratic in the photon field is given by while the one linear in the photon field is given by Carrying out color traces explicitly, one finds that (2.6) contains δ ab , while (2.7) has d abc structure. In that sense, (2.6) would exist for any choice of the gauge group for G µν including a U(1), while (2.7) requires N ≥ 3 for SU(N ), which includes of course the color group of the Standard Model. As we will see, below the QCD scale, operator (2.6) contributes to paramagnetic EDMs through the semi-leptonic CP -odd operator C S while (2.7) contributes to the neutron EDM. The implication of the CP -odd pure photon operator is studied in detail in [5] in the context of the muon EDM. This operator puts only a subdominant constraint on the heavy quark EDM, and thus we do not discuss it any further in this paper.

Light quark EDMs
The heavy quark EDM induces the up and down quark EDMs at the three-loop level through the diagram in Fig. 2 (and its permutations), where the closed loop is the heavy quark and the cross dot indicates the heavy quark EDM insertion. In our case the heavy quark EDM operator already contains a derivative acting on the external photon, and thus we can ignore the momentum flow of the photon. The mass and the external momentum of the light quark are small compared to the heavy quark mass scale, and this allows us to expand the diagram with respect to the light quark mass and momentum. To linear order, the amplitude is schematically given by where q = u, d is the light quark field and p is the light quark momentum. Here we used the equation of motion to express quantities that explicitly depend on the light quark mass m q in terms of p. After projecting S µνρ onto the EDM part, we obtain the light quark EDM d q as where the trace is taken only over the spinor index, and {.., ..} represents anticommutator. We evaluate this expression with the dimensional regularization. Since we expanded the amplitude with respect to the external momentum and the light quark mass, S µνρ depends only on the heavy quark mass and hence is given by the vacuum integrals. In our case, among six propagators, three are massive and three are massless. In this case, we can reduce the vacuum integrals to two irreducible master integrals [33]. We use FIRE6 [34] to reduce the vacuum integrals by integration by parts, and use analytical expressions of the master integrals in [35]. We then obtain where α s = g 2 s /4π and ζ(3) 1.202 with ζ(z) the Riemann zeta function. We have checked that all the gauge-dependent part of the gluon propagator cancels in our computation. Moreover, although each master integral is divergent, all the divergences cancel and the final result is finite. If we drop the color factors, this is equivalent to one contribution to the electron EDM from the muon/tau EDM, the diagram (a) in [3]. Our result differs from [3] and we will come back to the origin of this difference in a separate publication [36]. Notice that d q is suppressed only by 1/m Q since the loop integral is dominated by the loop momentum of order m Q , while the CP -odd photon-gluon operators that we computed in Sec. 2.1 are suppressed by 1/m 3 Q . Therefore the light quark EDM is expected to be more important for heavier quarks, although there is another (and a dominant) contribution in the case of the top quark as we mentioned at the beginning of this section.

Paramagnetic EDM
In this section we study a constraint on the heavy quark EDM from paramagnetic atomic/molecule EDM experiments, in particular the ACME experiment [2]. The paramagnetic EDM experiments are sensitive only to a linear combination of two operators, the electron EDM d e and the semi-leptonic CP -odd operator C S given by where G F is the Fermi constant, e is the electron (one should not confuse it with the electromagnetic coupling), p is the proton and n is the neutron, respectively. It is convenient to define an equivalent electron EDM as [37] d (equiv) where we focus on ThO molecule. The current bound on the equivalent electron EDM is given by |d (equiv) e | < 1.1 × 10 −29 e cm [2]. We now compute C S induced by the CP -odd photon-gluon operator, in particular the two-photon two-gluon operator, below the QCD scale. The effective Lagrangian is given by Because our starting point here is explicitly isospin symmetric, it will result in the same C S coupling for neutrons and protons. The nucleon matrix element of the gluon part is given by where N is the nucleon field (either p or n) with m N its mass, and · · · denotes the traceless tensor part that is irrelevant for our purpose. Here we used that, in the chiral limit, the one-loop trace anomaly dominantly contributes to the nucleon mass, where the coefficient in the right-hand-side is related to the beta function [38]. The photons are attached to the electron and induce the operatorēiγ 5 e at one-loop [39] as shown in Fig. 3. This integral is logarithmically divergent, which is regulated by the heavy quark mass scale. Therefore, to the leading logarithmic accuracy, we obtain where α = e 2 /4π and m e is the electron mass. By requiring |d (equiv) e | < 1.1 × 10 −29 e cm, we obtain for the charm quark and for the bottom quark, respectively. The heavy quark EDM also induces the electron EDM at threeloop [3,40], but its contribution to d (equiv) e is negligible compared to C S . We now estimate the precision of our calculation. There is an uncertainty associated with Eq. (3.5), which is valid only in the chiral limit. However, the quark contribution to the nucleon mass is less than 10 % [41], and hence the uncertainty of Eq. (3.5) is also less than 10 %. Thanks to the lattice computations, these contributions are rather precisely known and the uncertainty in the non-perturbative matrix element can be further reduced by including the light quark contributions into our computation. Another uncertainty originates from our photon-loop computation, where we include only the leading logarithmic term. The uncertainty associated with this treatment may be estimated by changing the cut-off scale from m Q to 2m Q , which results in 10 % change of the result for both the charm and bottom quarks. This can be further improved by performing the full two-loop computation, with the heavy quark and photon loops at the same time. Therefore, the precision of our C S calculation is estimated to be 10 % and this can be further improved by properly including the quark contribution to the nucleon mass and performing the full two-loop computation of the photon-loop.
Our constraint is stronger than the one derived from e + e − → qq at LEP by several orders of magnitude [25]. Although our constraint is weaker than [31] at its face value, there are two caveats. First, as we mentioned in the introduciton, the constraint in [31] is a re-interpretation of the bounds on the CEDMsd Q at the heavy quark mass scale, i.e.,d Q (m Q ). Both the quark EDM and CEDM at high energy scale Λ NP contribute to the CEDM at lower energy due to the RG running, and henced Q (m Q ) is given as a linear combination of d Q (Λ NP ) andd Q (Λ NP ). With only this information, one can put a constraint on only this particular linear combination of d Q (Λ NP ) andd Q (Λ NP ). Therefore the authors assume that there is no huge cancellation between d Q (Λ NP ) andd Q (Λ NP ), which allow them to derive a constraint solely on d Q (Λ NP ) and thus on d Q (m Q ). Since our constraint directly applies to d Q (m Q ), it is an independent constraint, without any assumptions on the relative size of d Q (Λ NP ) andd Q (Λ NP ). Second, it is a complicated task to estimate the size of the neutron EDM induced by the Weinberg operator. Since [31] relies on the constraints ond Q (m Q ) from the Weinberg operator, it has a large hadronic uncertainty. In contrast, as we have just emphasized, our computation of C S is precise and its uncertainty is well under control. Given the rapid progress of the paramagnetic EDM experiments, C S will continue providing important and clean constraints on the heavy quark EDMs. In the next section, we will see that the neutron EDM provides a stronger constraint than C S , comparable to [31], though with a larger hadronic uncertainty.

Neutron EDM
In Sec. 2 we have seen that the heavy quark EDM generates the CP -odd photon-gluon operator and the light quark EDM. Below the QCD scale, these CP -odd operators in turn generate the neutron EDM d n , whose experimental upper bound is [14] 1 |d n | < 1.8 × 10 −26 e cm. (4.1) In this section we translate this upper bound to the bound on the heavy quark EDM by estimating the size of the induced neutron EDM, with effects of the nonperturbative QCD taken into account (to the extent it is possible).

Light quark EDM contribution
We first study the neutron EDM induced by the light quark EDM d n . Both the QCD sum rule and the quark model suggest that [1,42,43] QCD sum rule calculations [42,43] and more recently lattice calculations [44] give support to this simple formula within ∼ 30% accuracy. We thus obtain As we will see, there is another contribution to the neutron EDM from the CP -odd photon-gluon operator that leads to the QCD condensate power corrections in the light quark propagator. The rest of this section is devoted to estimate the size of this contribution.

CP -odd photon-gluon operator contribution
In Sec. 2.1 we saw that the heavy quark EDM generates the CP -odd photon-gluon operator Below the QCD scale the gluons confine and condense, and this operator feeds into the neutron EDM. We use the QCD sum rule technique [45] to estimate the size of the neutron EDM induced by this operator.
The starting point of the QCD sum rule is to define the two-point correlator The interpolating function η is typically chosen as and has an overlap with the neutron one-particle state. The functions j 1 and j 2 have the same quantum numbers as the neutron and are given by where i, j, k are the color indices with ijk the totally anti-symmetric tensor, and C is the charge conjugate matrix (see App. A for its definition). Finally β is a redundant parameter that should disappear in an exact computation. In practice, one can choose β, e.g., to optimize the convergence of the operator production expansion (OPE). The QCD sum rule relies on the quark-hadron duality and evaluate this two point correlator Π(p) in the hadronic (phenomenological) side and the quark (OPE) side. On the phenomenological side, Π(p) is expressed in terms of hadronic quantities such as the neutron mass and dipole moments. On the OPE side, Π(p) is evaluated in terms of perturbative quarks and gluons with non-perturbative effects included in the form of the QCD condensates. One then obtains an estimation of the hadronic quantities by equating these two expressions, after the Borel transformation to reduce effects of excited states. In our case, we have external electromagnetic field F µν multiplying three Lorentz structures that in principle we can use to derive the sum rules: / pσ µν / p, { / p, σ µν } and σ µν . In this paper, we focus on the sum rule that follows from / pσ µν / p since it depends most strongly on the momentum and the vacuum susceptibilities are less relevant [46]. Moreover, it is computationally simple to derive the QCD sum rule based on this structure. In the following we evaluate Π(p) in the OPE and phenomenological sides, respectively. OPE side. We denote the light quark propagator as where only the color-diagonal contribution is relevant for the QCD sum rule based on / pσ µν / p in our case. The CP -odd photon-gluon operator does not distinguish the up and down quarks, and thus we collectively denote u and d as q. With this expression, we obtain where S c αβ ≡ CS T C αβ and S T αβ ≡ S βα with α and β being the spinor indices, and the trace is taken only over the spinor index. In our case the light quark propagator has two contributions where S (0) corresponds to the free quark propagator and is given by while S (CP ) is the CP -odd part which we compute from now. The source inserted into the correlator has three gluon fields, and we follow the general idea of [23] where the Weinberg operator contribution to d n was first evaluated, and where two gluons are treated perturbatively, while the third gluon field contributes to the quark-gluon condensate. With the QCD condensate background, we can write down the diagram where the crosses indicate the background fields, either the external photon or the QCD condensation of the quarks and gluons, and the cross dot is the insertion of the CP -odd photon-gluon operator (4.4).
Since the QCD vacuum does not violate the Lorentz and color symmetries, we have where σ · G = σ µν G µν , and · · · corresponds to the vacuum expectation value. We thus obtain (4.14) Its Fourier transformation has an IR divergence, which in the dimensional regularization leads where we take d IR = 4+2 IR and keep only the logarithmic terms, with Λ IR the IR cut-off scale. Due to the sensitivity to IR scale, calculations of further terms in the OPE are not possible. As a consequence, this evaluation should be viewed as an estimate, which cannot be systematically improved within QCD sum rule method. With the explicit forms of S (0) and S (CP ) , we obtain (4.16) where we perform the dimensional regularization to tame the UV divergence with d UV = 4 − 2 UV , and with µ the renormalization scale. We then perform the Borel transformation, defined as [47] B Π( We thus obtain the / pσ µν / p part as Taking the imaginary part of I(p 2 ; IR , UV ) is somewhat subtle, and we provide details in App. B.2. There we also clarify the physical meaning of the scale Λ IR ; it should be identified with the mass of the constituent quarks. This expression is to be compared with the phenomenological expression.
Phenomenological side. On the phenomenological side, the two-point correlator with the neutron EDM insertion is expressed as where m n is the neutron mass and λ n parametrizes the overlap between the interpolating function η and the neutron one-particle state. After the Borel transformation, we obtain where · · · expresses contributions from excited states which we ignore in the following.
Sum rule. The QCD sum rule of the neutron EDM is obtained by equating Eqs. (4.19) and (4.21).
We thus obtain is the neutron EDM induced by the CP -odd photon-gluon operator (to distinguish it from the one induced by the light quark EDM). We may use the Ioffe formula for the nucleon mass [48,49] to eliminate λ n . Then we obtain the QCD sum rule estimation of the neutron EDM where we used qσ · g s Gq = m 2 0 qq with m 2 0 = 0.8 GeV 2 [50].

Constraint on heavy quark EDM
The neutron EDM has two contributions, induced by the light quark EDM and the CP -odd photongluon operators, and is given by and we take β = −1 following [48,51,52]. The parameters should be evaluated at m Q for the former contribution, and at the scale close to Λ QCD for the latter contribution. We use m c = 1. Finally we take α s = 0.5, Λ IR = 300 MeV and M = 800 MeV for the QCD sum rule estimation for definiteness. For these values, the contribution from the CP -odd photon-gluon operator is larger by a factor of 16 and 10 for the charm and bottom quarks, respectively. In the bottom quark case, the explicit quark mass suppression is compensated by the running of the strong coupling, and the CP -odd photon-gluon operator is still larger than the light quark EDM even with the relative suppression factor 1/m 2 b . By requiring that |d n | < 1.8 × 10 −26 e cm, we obtain |d c | < 6 × 10 −22 e cm, (4.27) for the charm quark, and |d b | < 2 × 10 −20 e cm, (4.28) for the bottom quark, respectively. The constraint from d n is stronger than that from C S by more than an order of magnitude. However, one should note that the estimation of d n has a large hadronic uncertainty. Indeed, the final result is affected by a factor two if we use the sum rule of the nucleon kinetic term instead of the mass term for λ n . There are also uncertainties related to the choice of M and β which can again result in a factor of a few difference in the final result. Therefore our constraint here should be understood as an order-ofmagnitude estimation and the numerical factor should be taken with care. In contrast, our calculation of C S is far more precise. Its uncertainty is estimated to be ∼ 10 % and can be further reduced if needed (see the end of Sec. 3). In this sense, the constraints from d n and C S are complementary to each other; d n puts a stronger constraint on the heavy quark EDM, while the calculation of C S is cleaner and its uncertainty is well under control.
Our constraint on d c is stronger than [31] by a factor two, while the one on d b is weaker by a factor two. However, as we noted in the introduction and the end of Sec. 3, our constraint directly applies to d c (m c ) and d b (m b ), and is independent of the one in [31].

Conclusion
In this paper, we have derived indirect limits on the charm and bottom quark EDMs. The charm and bottom quark EDMs generate the CP -odd photon-gluon operators and the light quark EDMs after integrating out the charm and bottom quarks. Photon-gluon operators contribute to the semi-leptonic CP -odd operator C S (and ultimately to paramagnetic AMO EDMs) as well as to the neutron EDM at a non-perturbative level. Quark EDM dominantly contributes to the neutron and nuclear EDMs. Performing our evaluation and using the current limits, we obtain from the neutron EDM experiment, respectively. Although the constraint from the neutron EDM is stronger, it has a larger hadronic uncertainty. The uncertainty of the constraint from the paramagnetic EDM is estimated as 10 % and can be improved if needed, while the uncertainty from the neutron EDM can be a factor of a few. Our constraint is independent of the one given in [31] in the sense that our constraint directly applies to the EDM operators at the quark mass threshold. By assuming a simple scaling of d Q /e ∝ (α/π)m Q /Λ 2 Q , we may translate our constraint as a lower bound on CPodd new physics scale: Λ c > 70 GeV and Λ b > 20 GeV from the paramagnetic EDM experiment, and Λ c > 300 GeV and Λ b > 100 GeV from the neutron EDM experiment.
Our result provides an important benchmark to overcome for the LHC based measurements of the charmed baryon EDMs [26][27][28][29][30]. The idea of using the bent crystal technique for studying electromagnetic properties of baryons containing a heavy quark is very appealing. However, given the strength of the bounds derived in our work, and the necessity to satisfy independent constraints from d n and C S (hence removing a chance of accidentally large d c(b) due to cancellations), one may want to re-evaluate the main goal of the charmed baryon experiment. While it will be difficult to match the indirect sensitivity to d c(b) , the planned measurement may achieve sufficient accuracy to extract the values of the magnetic moments µ c(b) and compare it with the QCD predictions.
Acknowledgements Y.E. and M.P. are supported in part by U.S. Department of Energy Grant No. de-sc0011842. M.P. would like to thank Dr. K. Melnikov for the advice in evaluating loop contributions. The Feynman diagrams in this paper are generated by TikZ-Feynman [54].

A Convention
Here we summarize our conventions used in this paper. The field strengths are given by where f abc is the SU(3) structure constant. The SU(3) generator satisfies The dual field strengths are defined as with 0123 = +1. We take the gamma matrix as and σ µν as The charge conjugation matrix C satisfies such that ψ T C has the same Lorentz transformation property asψ. In the Weyl representation it is given by

B Technical details
In this appendix we provide some technical details that we omit in the main text.

B.1 Derivation of CP -odd photon-gluon operators
In this subsection we provide derivation of the CP -odd photon-gluon operators (2.5) after integrating out the heavy quarks (see [55] for an extensive review on this procedure). Our starting point is Eq. (2.3). We can simplify it as where (iD) 2 = iD µ iD µ and we used that the traces of odd γ's vanish in the second equality. We further expand the denominator with respect to H, and obtain up to fourth-order in the gauge fields The term linear in σ µν vanishes after taking the trace of the spinor index. In order to compute T 2 , we use the following identity [55]: where f (H, F ) is an arbitrary function of H µν and F µν . We thus obtain where we ignore the derivatives acting on F . We can now replace iD by i∂ in the last term to the order of our interest, and obtain (B.8) By integrating this we obtain where the additional mass scale M R comes from the regularization which we did not write down explicitly. Next we compute T 3 . By ignoring the derivatives acting onF µν , we obtain where the dots indicate terms that contain more than two commutators and thus correspond to higher dimensional operators that are out of our interest. The first term is easily seen to generate only O(H 4F ) operators with the help of the identity (B.6). The second term is expanded as The first two terms induce only higher order terms [55], and thus we ignore them. The third term in Eq. (B.10) already contains three field strengths and two commutators, and hence we can ignore any further commutators between D and H. By taking the angular average, we obtain By combining them, we obtain Finally we compute T 4 . To the order of our interest, we can simply replace iD by i∂ in the denominator. It is then easy to see that (B.14) Therefore we obtain where we ignored the quadratic term that is irrelevant for our purpose. If we drop the gluons, it correctly reduces to the result in [5].

B.2 Borel transformation and IR divergence
In this subsection we discuss the Borel transformation of Eq. (4.16) that contains both the UV and IR divergences. As the Borel transformation is related to the imaginary part, it is equivalent to taking the imaginary part of (B.16) The limit of this function at UV → 0, IR → 0 is not well defined. For example, if the limit is taken along the line UV = a IR with fixed a, one would get an a-dependent result: (B.17) where we only kept the double logarithmic terms. The purpose of this subsection is to understand the correct prescription of evaluating the imaginary part of this function.
In order to understand the correct prescription, it is helpful to consider a simpler example: a scalar three-body decay. Indeed, the two point correlator Π(p) contains three quark propagators, and thus the imaginary part of this function is related a three-body decay phase space integral. Therefore we consider the following Lagrangian where we omit i for notational simplicity. In the coordinate space this is given by (B.24) The first order term in m 2 χ is IR divergent and this is analogous to our propagator with the CP -odd operator insertion in Sec. 4.2. The two-loop amplitude is given by To the first order in m 2 χ we obtain where we focus on the leading logarithmic term on the right hand side. Thus the correct prescription of evaluating the imaginary part of I is Im I(p 2 ; IR , UV ) = −π log Λ 2 IR p 2 , (B.28) where we identify Λ 2 IR = m 2 χ in the present case. This agrees with the Borel transformation formula in [52]. For the example in Eq.(B.17), this prescription is equivalent to neglecting the real part of the UV logarithm.
Our discussion clarifies the physical meaning of Λ IR that appears in Eq. (4.19). This IR divergence originates from the phase space integral and is regulated by the mass of the daughter particles. In the neutron EDM case, the daughter particles are the constituent quarks. Therefore Λ IR is identified with the mass of the constituent up and down quarks that we take Λ IR = 300 MeV in the main text.