Electric dipole moments from colour-octet scalars

We present the contributions to electric dipole moments (EDMs) induced by the Yukawa couplings of an additional electroweak doublet of colour-octet scalars. The full set of one-loop diagrams and the enhanced higher-order effects from Barr-Zee diagrams are computed for the quark (chromo-)EDM, along with the two-loop contributions to the Weinberg operator. Using the stringent experimental upper limits on the neutron EDM, constraints on the parameter space of this model are derived.


Introduction
The Standard Model (SM) has accurately anticipated a wide range of phenomena and satisfactorily described almost all experimental outcomes, resulting in the best description of nature that we have up to date. Nevertheless, it has clear shortcomings which require invoking new physics (NP). One of them is the observed baryon asymmetry of the universe, which differs by several orders of magnitude with respect to the SM prediction. To generate this asymmetry, the Sakharov conditions [1] require additional sources of CP violation (CPV) beyond the SM that, in turn, can be constrained by the stringent experimental limits on the permanent electric dipole moment (EDM) of particles. Due to the small size of EDMs as predicted by the SM 1 , any signal of a nonzero EDM in current or planned experiments would be an indisputable sign of NP. In this sense, these observables provide a background-free search for new sources of CPV beyond the SM.
One of such sources can be found in extensions of the SM with additional scalar particles that transform as doublets under SU (2) L and as octets under SU (3) C . These scalars were first proposed by Manohar and Wise (MW) [4], the original motivation being that they are one of the few scalar representations of the SM gauge group that can implement Minimal Flavour Violation (MFV) [5,6]. In addition, these scalars emerge naturally with a mass of few TeVs from SU (4), SU (5) or SO(10) unification theories at high energy scales [7][8][9][10][11][12][13].
In this work, we study the phenomenology of EDMs arising from the Yukawa couplings of these scalars, complementing the extensive literature of phenomenological studies within this model .
Only few publications have studied the effect of colour-octet scalars on the EDM of particles [28,[36][37][38]. The restrictions on their Yukawa couplings from heavy quark EDMs were studied in Ref. [28] and numerically updated with more restrictive bounds in Ref. [32]. However, to derive robust limits on the model accounting for cancellation effects, the direct contributions to the neutron EDM through gluonic or light-quark operators need to be considered. The relevant contributions include two-loop diagrams which, to our knowledge, have not been computed in the literature within the MW model. Due to the colour structure of the new scalars, the light-quark EDMs are greatly enhanced when compared to the contribution from colourless scalars appearing in the two-Higgs-doublet model (THDM), as shown in Figure 1. This feature makes the hadronic EDMs powerful observables to assess the viability of the MW theory.
This paper is organised as follows. First, in Section 2, we describe the relevant operators in the low-energy effective field theory that will contribute to the EDM observables. There, we also describe how to match the contributions of the new degrees of freedom to the operators at the hadronic scale. The contributions of these operators to the hadronic and nuclear EDMs, together with their current experimental limits, are presented in Section 3. In Section 4 we briefly introduce the MW model. The main results of our paper are shown in Sections 5 and 6. First, in Section 5, we provide analytical expressions for the oneand two-loop contributions of coloured scalars to the chromo-EDM (CEDM) and EDM of the quarks, as well as to the Weinberg operator. Furthermore, we present a detailed  Figure 1. Opposed to the THDM with colourless scalars (left), the leading contribution to the light quark EDM appears in the MW model through gluon exchange (right), enhancing the EDM by a factor (C F α s )/α ∼ 50, with the couplings evaluated at the hadronic scale.
comparison between the contributions to the quark EDM for each flavour. In Section 6 we provide constraints on the CP-violating couplings of the colour-octet scalars from the neutron EDM, whose competitiveness with other observables constraining the same planes of the parameter space is demonstrated. Finally, in Section 7 we summarise our results.

Effective theory framework
In this section, we provide the effective theory framework that describes the relevant operators contributing to the nuclear and atomic EDMs and their evolution to low energies through the renormalisation group. Below the electroweak scale µ < Λ EW the relevant flavour-diagonal CP-violating effective Lagrangian is given by 2 where the sum of q runs over all quark flavours. The effective operators are defined as Here, F µν and G a µν with a = 1, ..., 8 are the electromagnetic and gluon field strength tensors, g s is the strong coupling constant (α s ≡ g 2 s /4π), and σ µν = i 2 [γ µ , γ ν ]. The matrix T a represents the generators of the SU (3) C group with normalisation Tr(T a T b ) = δ ab /2, and the tensor f abc the structure constant. The charge of up-and down-type quarks is Q q = (2/3, −1/3). The expression for the covariant derivative will be relevant to define the anomalous dimension matrix later in Eq. (2.5). In this work, it is defined as µ T a , where A µ and G a µ are photon and gluon fields, respectively. We do not consider the contributions from four-quark operators since their moderate effect is well below the two-loop contributions as shown in Section 5.3. The quark EDM d q (µ), chromo-EDM d q (µ), and the usually defined coefficient w(µ) of the Weinberg operator are related to the Wilson coefficients by To compare the theory prediction with the low-energy EDM observables, these parameters need to be evaluated at the hadronic scale, µ had ∼ 1 GeV. For that purpose, we employ the renormalisation group equations (RGEs), , and γ(µ) is the anomalous dimension matrix. At leading order in α s it reads [38][39][40][41][42][43][44][45] and n f denotes the number of active flavours. Solving Eq. (2.4), we obtain the scale dependence of the Wilson coefficients − → C (µ) for a theory with constant number of active flavours. Starting at the NP scale µ = Λ NP , close to the top quark mass, 3 where the fundamental theory is matched to the effective one, the Wilson coefficients are evolved down to the bottom-quark mass scale with n f = 5. At this point, the bottom quark is integrated out, generating a threshold contribution of the bottom CEDM to the Weinberg operator as [40,46,47] Here µ + b and µ − b refer to the scale µ b ∼ m b in the theories with n f = 5 and n f = 4, respectively. Analogously, also the charm CEDM induces a threshold correction to the Weinberg operator, although it is numerically irrelevant for our study on the MW model. The final running with n f = 4 and n f = 3 brings the Wilson coefficients down to the hadronic scale µ had ∼ 1 GeV.
It is interesting to note that also the EDM of heavy quarks contributes to the hadronic EDMs [32]. This occurs through photon-loop corrections to the corresponding quark CEDM which, in turn, provides threshold corrections to the Weinberg operator, as mentioned above. If family-specific couplings exist that appear only in heavy quark EDMs, e.g. in models of 3 The masses of the new scalars are in fact constrained to be above 1 TeV [33,35]. By simultaneously integrating out the new scalars and the top quark, we ignored the running in the range [1 TeV, mt]. Conservatively, we estimate this effect to be at most of 30%, which is smaller than the leading systematic error from the hadronic matrix elements in Eq. (3.1). scalar Leptoquarks, this effect, suppressed in comparison to gluon-loop corrections, provides a unique window to access these family-specific couplings [32]. In the model considered in this work, we have checked that the leading constraints to the Yukawa couplings appear from the light quark (C)EDMs, as we will see in Section 6, and therefore we neglected O(α) corrections in Eq. (2.5) for the sake of simplicity.

Matching onto nuclear and atomic EDMs
The contributions of the operators in Eq. (2.2) at the hadronic scale to the electric dipole moments of the neutron d n , proton d p , and mercury d Hg are computed with non-perturbative techniques of strong interactions at low energies. State-of-the-art coefficients for the CEDMs and Weinberg operator have been obtained in the literature with QCD sum rules [48][49][50][51] and the quark model [52], while the contributions of the quark EDMs have been computed in lattice QCD [53][54][55][56][57], with significantly smaller errors. Assuming a Peccei-Quinn mechanism, these read [58] where the contributions to d Hg from pion-nucleon interactions that are compatible with zero within 1σ have been left out. In these expressions, the tensor charges describing the light quark EDM contributions read g u T = −0.213 ± 0.012 and g d T = 0.820 ± 0.029. The current experimental limits of d n , d p , and d Hg are collected in Table 1 [59], and mercury d Hg [60], and the future experimental sensitivity, including the proton d p in e cm units.
The current bound on the Mercury EDM [60] is four orders of magnitude stronger than that of the neutron, which has been recently reported by the nEDM collaboration [59]. However, this difference is compensated by the suppression factor in front of d n,p in its contribution to d Hg , as shown in Eq. (3.1). As a result, we obtain almost identical constraints on the model parameters by using the d n or d Hg experimental limits (within less than a 10% difference), and in the numerical analysis of Section 6 we will use only the direct limit on d n .

The scalar colour-octet model
The inclusion of a scalar field transforming as (8, 2) 1/2 under the SM gauge group SU (3) C × SU (2) L × U (1) Y allows the construction of more renormalisable invariant interactions [4]. In general, the Lagrangian describing colour-octet scalar interactions can be written as where L SM , L kin , L Y , and L S represent the SM Lagrangian, the colour-octet scalar kinetic term, the interaction with SM fermions, and the scalar sector, respectively. The kinetic term generates interactions with the SM gauge particles through the covariant derivative The scalar sector L S encodes the interactions between the octet scalars which are given by [4] GeV the vacuum expectation value (VEV). Here, i and j are SU (2) L indices, and all traces are in colour space. It can be seen from Eq. (4.3) that λ 3 , λ 4 and λ 5 can have complex phases. However, through an appropriate phase rotation of S, λ 3 can be chosen to be real. The phenomenological study of the present work will be limited to EDM observables arising from the CPV in the Yukawa couplings. The effect of the CP-violating parts of the potential will be considered in a future global-fit analysis using several observables at the same time.
The VEV, v, causes a mass splitting between the octet scalars. Decomposing the neutral octet scalars into two real scalars, one obtains the following relation between the physical masses where m S ± is the mass of the charged scalar, m S 0 R the mass of the neutral CP-even scalar and m S 0 I the mass of the CP-odd scalar. Assuming MFV, the new Yukawa interaction can be parametrised by two complex numbers inducing new CP-violating sources beyond the SM that contribute to the (C)EDM of quarks. In Eq. (4.6), Q L represents the left-handed quark doublet, and u R and d R correspond to the right-handed up-and down-quarks singlets, respectively. The new scalar fields are written as S = S a T a with S a = (S a,+ , S a,0 ) T . In Eq. (4.6), the shorthand notation S = i σ 2 S is employed, where σ 2 is the usual Pauli matrix.

Contributions to EDMs in the colour-octet scalar model
In this section, we analyse the different contributions to the neutron EDM from the colouroctet scalars. Namely, we will derive the expressions for the quark (C)EDM at one-loop level and the enhanced contributions at two-loop level, together with the leading twoloop contributions to the Weinberg operator. The constraints on the model using these expressions are discussed in Section 6.

One-Loop contributions
At one-loop level, the (C)EDM of a quark q receives contributions from neutral and charged scalars, S 0 I,R and S ± , as shown in Figure 3 ( Figure 2). These contributions can be computed using standard techniques, and are finite since the Lagrangian does not contain any treelevel (C)EDM. Our results for the CEDM of a quark q reads are colour factors emerging from the color structures appearing in the Feynman diagrams The loop functions 4 Eqs. (5.2) and (5.7) can be translated to the basis of Ref. [28] through the following replacements dq → − dq and dq → − dq, as well as ηU → ηU e i(α U +π) and ηD → ηD e −iα D . 5 Sum over repeated (colour) indices is understood here and in the following.
F n,m (r) and G n,m (r 1 , r 2 ) are defined in the notation of Ref. [28] as  While both neutral and charged scalars couple to gluons, only charged scalars couple to photons, making the CEDM receive an extra contribution with respect to the EDM, depicted in diagram (b) of Figure 2. The contributions to the EDM of a quark q, shown in Figure 3, share the same loop functions with the CEDM diagrams. Therefore, the results for the EDM can be given in terms of the CEDM expressions (Eqs. (5.2)) as and The colour factor C F emerges from the combination of two colour matrices (T a T a ) ij = C F δ ij , each of them provided by one of the two Yukawa couplings in the diagrams of Figure 3.
Correcting by colour factors, our results are in good agreement with the literature of the colourless THDM [61,62]. Additionally, for the diagrams that do not appear in the THDM but emerge through the introduction of colour-octet scalars (see Figure 2 (b) and (d)) we found agreement with the previous calculation in Ref. [28]. However, in this reference, the loop function F 2,0 in Eqs. (5.2) is replaced by F 0,0 , which differs from our results and those of Refs. [61,62].

Two-loop contributions
In the previous section, we have studied all one-loop contributions to the (C)EDM of the quarks. Since light quark (C)EDMs are heavily suppressed at one-loop level by powers of the quark masses, the leading contributions to the neutron EDM will happen at two-loop level. In the following, we derive these two-loop contributions to the quark (C)EDM and the Weinberg operator.

Barr-Zee diagrams
Although being suppressed by additional coupling constants and loop factors, the Barr-Zee type diagrams, shown in Figure 4, benefit from the enhancement of the top-quark Yukawa coupling in flavour models. Among the contributions to the CEDM, the one depicted by diagram (b) (Figure 4) largely dominates, since it is enhanced by the strong coupling constant (from the internal gluon propagator) and it is not suppressed by any mass of the electroweak bosons. For this reason, this will be the only diagram that will be included in the expressions of the CEDM. Contributing to the quark EDM, only the diagram (a) (Figure 4) is numerically relevant.  The Barr-Zee contribution to the CEDM is given by where we have defined the loop functions 6 , (5.10) , (5.11) and the colour factor N qg The first two traces correspond to the current flow of the top quark in the loop, clockwise and counterclockwise. Similar to the one-loop case, the EDM depends on the same loop functions as the CEDM and they are related by with Q t = 2 3 being the charge of the top quark. For details on the computation of these diagrams see Appendix B. The Weinberg operator can also play an important role since it contributes directly to the neutron EDM and does not suffer from light quark mass suppressions. Furthermore, it also has an impact on the light quarks' (C)EDM due to the operator mixing in the RGEs. The first-order contribution appears at the two-loop level via the exchange of neutral 6 F (1) and F (1) are related to f and g from Ref. [62], through F (1) ( √ r) = −f (r) and F (1) ( √ r) = g(r) .

Weinberg contribution
are the corresponding loop functions. In turn, the diagram (c) of Figure 5 vanishes. The details of the calculation are found in Appendix A. Notice that h(r) corresponds to the wellknown Weinberg loop function of Ref. [63] but g(r) is only appearing for coloured scalars since they can couple to gluons. Looking at Eq. (5.13) we see how, as it happened for the one-loop contribution to the (C)EDMs of the quarks, there is a relative minus sign between the contribution of the CP-even and CP-odd neutral scalars. Therefore, this contribution will be suppressed by the mass splitting of the neutral scalars which, as can be seen in Ref. [35], is around two orders of magnitude smaller than the mass of the scalars. This suppression is large and, for the numerical analysis in Section 6, we will assume that all scalar masses are the same, effectively neglecting the contribution of the neutral scalars in the Weinberg operator.
The calculation of the charged scalar contributions, shown in Figure 6, proceeds differently. Having a bottom quark propagator in the loop, these diagrams will only induce a contribution to the Weinberg operator below the bottom-quark mass scale. At this scale, the NP particles and the top quark have already been integrated out and their information is encoded in the effective vertices shown in diagrams (c) and (d) of Figure 6. These oneloop diagrams generate a threshold contribution to the Weinberg operator from the bottom CEDMd b , as shown in Eq. (2.6). The main contribution tod b comes from the one-loop diagrams (c) and (d) of Figure 2, which are not suppressed by any kind of light quark mass or CKM factor.

Four-quark contributions
Four-quark interactions due to the exchange of colour-octet scalars (Figure 7) have been studied in detail in Ref. [38]. Using the results from that work, we observe that they represent a sub-leading contribution to the neutron EDM well below the two-loop contributions presented in Section 5.2. Therefore, we neglect them in the following.

Discussion
Before moving on with the phenomenological analysis to compare the model predictions against current experimental upper limits, it is worthwhile to compare the size of the Barr-Zee contributions, obtained above, to the one-loop contributions presented originally in Ref. [28]. 8 Due to the strong suppression from powers of the light quark masses, the light quark (C)EDMs are dominated by two-loop Barr-Zee diagrams. Even though these contributions include additional coupling constants and loop suppression factors, the top-quark Yukawa coupling, proportional to m t , makes the two-loop diagrams the dominant contribution for light quark (C)EDMs. With the same argument, one would expect heavy quark (C)EDMs to be dominated by one-loop diagrams, as they are not suppressed neither by light quark masses nor by loop suppression factors. This hierarchy of contributions is illustrated in

Phenomenological analysis
With all the relevant contributions of the coloured scalars to the EDM of hadrons obtained, we can study the constraints that these observables impose on the model. Currently, there is no direct limit on the EDM of the proton and it will not be used for this analysis. Furthermore, as commented previously, the implications from the neutron and mercury EDM on the MW model are extremely similar. Therefore, in order to provide clearer limits, we will only consider the direct limit on the neutron EDM in the numerical analysis. The values for the hadronic and nuclear matrix elements of Eq. (3.1) have been set to the central values. Of course, with only one observable and a total of seven parameters, |η U |, |η D |, arg(η U ), arg(η D ), m S R , m S I , and m S ± , we do not provide a complete phenomenological analysis here, but just a brief study showing the potential of EDM observables to constrain the parameter space of the model. To this end, the interplay of the model parameters appearing in the neutron EDM is discussed, comparing its limits to those imposed by other powerful observables studied in the literature. A global fit with all the relevant observables for the CP-violating MW model, although interesting, is beyond the scope of this paper.
To reduce the total number of free parameters and provide some sensible plots, we will first assume that the mass of the scalars is degenerate, m S R = m S I = m S ± = m S . This assumption is completely reasonable given the constraints on the mass splitting found in Ref. [35]. In addition, we will show here the results for |η U | = 1, which is below the maximum allowed value found in Ref. [35] for masses of the coloured scalars of a few TeVs. If |η U | = 0, the contributions with the top quark Yukawa coupling studied here vanish, and only the suppressed diagrams with bottom-quark propagators running in the loops contribute.

Neutron EDM predictions
An intuitive way to see the effect of the neutron EDM bound on the MW model is to compare its prediction, as a function of the model parameters, with the current experimental limit. Besides giving clues on the allowed size of the model parameters (studied in more depth in Section 6.2), this also allows to evaluate the effect of future neutron EDM limits on this model. In Figure 9 we show the size of the neutron EDM as a function of arg(η U ), fixing |η U | = 1, |η D | = 10 and arg(η D ) = 0, for different values of the mass of the coloured scalars.
Here we see how a strong constraint for arg(η U ) can be obtained even for masses of the scalars around 1.5 TeV, using the current experimental limits for the neutron EDM. It is also interesting to look at the possible constraints on the scalar masses, obtained by fixing the other parameters. This can be seen in Figure 10 where we have fixed arg(η U ) = π/2 (which gives the strongest contribution to the neutron EDM), |η U | = 1 and arg(η D ) = 0, and we have varied |η D |. Here we can see how for reasonable values of |η D |, the mass of the scalars could be constrained to be higher than 3 TeV, far beyond the current experimental limit from direct searches. Figure 9. Electric dipole moment of the neutron as a function of the complex phase of η U . The shaded region is excluded by the current experimental limit.

Constraints on the model parameters
To assess the restrictive power of the neutron EDM bounds, we can compare them to the most restrictive observables on the same planes of the parameter space. Following the global-fit analysis of Ref. [26], we chose the observable B(B → X s γ) as a benchmark to study the EDM restrictions. This comparison is done in Fig. 11, where the region of the parameter space allowed by each observable is shown. In the following we describe these results, pointing out the main patterns in this figure.
We have fixed the phases arg(η U ) = 0 (arg(η D ) = 0) in the top (bottom) panels in order to study the effect of CP violation as coming from the down-type (up-type) Yukawa couplings. First, looking at the |η D | − arg(η U ) plane (top-left panel) one can see that the constraints from the neutron EDM are stronger than those of B(B → X s γ). The only exceptions lie in the vicinity of the values arg(η U ) = 0, ± π, where d n vanishes and it cannot impose any restriction on the model parameters. Fortunately, an excellent experimental precision on B(B → X s γ), which is sensitive to both CP-violating and CP-conserving interactions, allows to restrict those directions even for these limiting values of the phases. This feature shows the power of combining the stringent experimental limits on EDMs with the complementary information from flavour observables. Nonetheless, as in any interaction beyond the SM, when the absolute value of the new coupling is small enough, restrictions on other model parameters cannot be found. In our case and with the current experimental precision of the neutron EDM, this feature appears as an horizontal band at |η D | 1, spanning along all the domain of arg(η U ) in the top-left panel.  Figure 11. Constraints on the parameter space of the MW model from the limits on the neutron EDM (blue) compared to those from B(B → X s γ) (red), calculated in Ref. [26]. The coloured areas represent the allowed regions of parameter space. In the top panels one can see the large impact of the new EDM bounds on some regions of the parameter space. However, its restrictive power is diminished either when the absolute value of the Yukawa couplings is small enough, when the phases are arg(η U,D ) ≈ 0, ±π, or in the presence of cancellation effects as in m S ≈ 1 TeV (see text for details).
Regarding the |η D | − m S plane, the constraints from d n are much stronger than those from B(B → X s γ) when arg(η U ) = π/2 and arg(η D ) = 0 (top-right panel in Fig. 11). However, exchanging the values of the phases, arg(η U ) = 0 and arg(η D ) = π/2 (bottom-right panel), an unconstrained direction appears for masses around m S ∼ 0.9 TeV (bottom-right panel). In this specific region of the parameter space, different contributions to the neutron EDM cancel out. In particular, the light quark (C)EDMs (through Barr-Zee diagrams) and the Weinberg operator (through the threshold contribution proportional to the bottom CEDM,d b ) have similar sizes. For arg(η U ) = π/2 and arg(η D ) = 0, both contributions interfere constructively resulting in very stringent limits (top-right panel). Conversely, when the values of the phases are switched, arg(η U ) = 0 and arg(η D ) = π/2, the Weinberg contribution flips sign, and the interference becomes destructive, preventing any constrain from the neutron EDM. To illustrate this dilution of the constraints, we fixed the mass value close to where the destructive interference is produced, m S = 1 TeV, and plotted the allowed regions in the |η D | − arg(η D ) plane (bottom-left panel), keeping arg(η U ) = 0.

Summary
In this paper we have analysed the relevant contributions to the neutron EDM in the MW model. Expressions for the quark (C)EDM and Weinberg operators have been obtained, which can easily be generalised to other models with colour-octet scalars through the appropriate relations between the coupling constants. In the case of the Weinberg operator, the neutral scalar contributions turn out to be irrelevant due to the cancellation between CPodd and CP-even scalars, with the charged scalar contribution being completely dominant for this operator. In turn, only the neutral scalars produce sizable effects in the (C)EDM of light quarks through Barr-Zee type diagrams.
Using the current experimental limits on the neutron EDM, we found new stringent limits on the parameter space of the MW model when the Yukawa CP-violating phases are different from zero. Additionally, in the presence of strong cancellations between the contributions to the neutron EDM, or when the Yukawa phases are zero, we found a valuable complementarity of the neutron EDM with other flavour observables. In future works, the combination of these observables in a global-fit analysis will lead to the most stringent limits on the general CP-violating MW model. areas of theoretical physics. We see some of his footprints in our analysis of EDMs as well: he formulated the CP-odd three-gluon operator, in Eq. (2.2), and its two-loop leading contribution through the exchange of a scalar field. In his original paper [63], he provided the h(r) loop function, while the full expression was obtained by D. Dicus in Ref. [64], who referred to the computation of this diagram as straightforward.
At various points in the calculation we found, however, that equally justified choices of parameterisation or approximations can take the parametric integral off the straight track towards this simple analytic expression. To our knowledge, these technical details are not found in the literature in a comprehensive summary. To facilitate the reproducibility of this analytical shape, we describe these details in the following. Following the same procedure, we arrived at the expression of g(r), which is surprisingly simple as well, and at the cancellation of the diagram in Figure 5 (c). The Dirac trace of this two-loop amplitude contains up to eight γ µ , and two γ 5 matrices. Since the final result is finite, the traces with γ 5 can be solved with the usual relation Tr(γ µ γ ν γ ρ γ σ γ 5 ) = 4iε µνρσ . The CP-violating parts of this amplitude are proportional only to the index structures with Levi-Civita tensors ε µ 1 µ 2 σρ , ε µ 2 µ 3 σρ , ε µ 1 µ 3 σρ , and ε µ 1 µ 2 µ 3 σ , where the indices σ and ρ are contracted with external momenta. To ease the calculation, it is convenient to select only one of these linearly-independent structures, as the final result is independent of this choice. The directions of the internal loop momenta were chosen as in Figure 12. With this, only three propagator denominators contain external momenta, which are small compared to the heavy mass M . Thus, we can expand these denominators in powers of (p 2 /M 2 ) with and carefully removing higher-order terms after the expansions. Once the denominator is free of external momenta p, the tensor integrals with an odd number of open indices vanish, and, for the rest, we can apply the identity k µ i k ν j → (k i · k j /D)g µν . The resulting master integrals have the shape To re-express the W functions in terms of Feynman parameters, one must use Feynman parameterisation of two denominators at a time, and the standard Wick rotation to sequentially integrate over the loop momenta k 1 and k 2 . The solution to the master integral that leads to the desired analytical shape of h(r) reads Selecting the tensor structure ε µ 1 µ 2 σρ , only the function W 11 contributes to this amplitude in the diagram of Figure 12. Finally, the three permutations of this diagram, obtained by rotating the internal scalar propagator by 120 o at a time, are obtained from the first diagram by renaming the indices and external momenta. This step is crucial to analytically cancel the divergencies of different parametric integrals. To obtain the Wilson coefficient, the fundamental amplitude must be matched onto the effective one, that reads [64], iM eff = − 2 3 f abc g s w ε µ 1 (p 1 ) ε µ 2 (p 2 ) ε µ 3 (−p 1 − p 2 ) (A.5) [(p 1 − p 2 ) µ 3 ε µ 1 µ 2 σρ + 2 (p 1 µ 2 ε µ 1 µ 3 σρ + p 2 µ 1 ε µ 2 µ 3 σρ )] p σ 1 p ρ 2 .
To develop these expressions, we used the help of the open-source packages FeynArts and FeynCalc [65][66][67].

B Barr-Zee diagrams
To simplify the calculation of the Barr-Zee diagrams, the two loops can be computed sequentially. The loop attached to the external photon (gluon) shall be obtained first. The result, in terms of Feynman integrals, can be written in the shape [68] iΓ µν g,γ = i(g µν k · q − k µ q ν )S g,γ + i µναβ k α q β S g,γ , (B.1) Only the top quark Yukawa coupling gives a sizeable contribution to this vertex. Thus, we only considered top quarks running in the inner loop, as shown in Figure 4. To arrive to this result, we use Feynman parametersiation in the shape of Eqs. (6.23) and (6.25) of the detailed guide for computations [69]. Furthermore, the photon (gluon) is assumed to be soft, i.e. k·p → 0, following the arguments of Ref. [68]. Once the expressions for the first loop are parametrised as in Eq. (B.1), this effective vertex is plugged in the second loop (Figure 13), rewriting the denominator k 2 (x−1)x+m 2 t as another propagator with momentum k. Then, the integrals over k can be identified in terms of Passarino-Veltman functions. Expanding the result in powers of (m q /M ), where M is a heavy mass and q = u, d, only the first term is numerically relevant. In this way, we obtained the loop functions F and F, in terms of the Feynman parameter x, which comes from the inner loop. To match the fundamental amplitude to the effective (C)EDM operator, it is convenient to express the Levi-Civita tensor in terms of products of gamma matrices, through the Chisholm identity.